[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: what's the difference between PetscViewerASCIIOpen() and PetscViewerBinaryOpen()?

thank you for your reply. Do you have any method to generate an ascii file of the huge sparse matrix? thanks 


I have been using HYPRE calls (from PETSc) to accomplish this.
Each processor will write his own  portion of the global matrix
in his own file.
The enclosed file should explain how.

Regards,
Aldo

--
Dr. Aldo Bonfiglioli
Dip.to di Ingegneria e Fisica dell'Ambiente (DIFA)
Universita' della Basilicata
V.le dell'Ateneo lucano, 10 85100 Potenza ITALY
tel:+39.0971.205203 fax:+39.0971.205160

      SUBROUTINE PrintMat(A,RHS,X,cols,values,maxcols,ITER,JOB)
C
      IMPLICIT NONE
C
CC#define HYPER_PARCSR 5555
C
C
C
C
C     $Id: printmat.F,v 1.3 2006/11/10 08:31:31 abonfi Exp abonfi $
C
#include "include/finclude/petsc.h"
#include "include/finclude/petscvec.h"
#include "include/finclude/petscmat.h"
#include "include/finclude/petscksp.h"
#include "include/finclude/petscpc.h"
#include "include/finclude/petscis.h"
#include "include/finclude/petscviewer.h"
C#include "HYPRE.h"
C#include "IJ_mv.h"

C
      Mat A
      Vec RHS,X
      INTEGER ITER,JOB
      INTEGER rstart,rend,ncols,i,k
      INTEGER HYPRE_PARCSR
      parameter(HYPRE_PARCSR=5555)
      INTEGER maxcols
      integer cols(maxcols)
      double precision values(maxcols)
      double  precision info(MAT_INFO_SIZE)
      INTEGER nz_a
C
C
C
C
C     .. Parameters ..
C     ..
C     .. Scalar Arguments ..
C     ..
C     .. Arrays in Common ..
C     ..
C     .. Local Scalars ..
      INTEGER IFAIL,MY_PE
C     ..
C     .. Local Arrays ..
      PetscViewer MyOpenMindedViewer
      INTEGER*8 ij
      CHARACTER* 6 matfile,rhsfile,solfile
      PetscScalar x_array(1)
      PetscOffset i_x
      PetscInt n
!     PetscErrorCode IFAIL

C     ..
C     .. External Functions ..
C     ..
C     .. External Subroutines ..
C     ..
C
C     .. Intrinsic Functions ..
C     ..
C     .. Common blocks ..
      COMMON /MPICOM/MY_PE
C     ..
C     .. Equivalences ..
C     ..
C     .. Data statements ..
      DATA matfile,rhsfile,solfile/"matXXX","rhsXXX","solXXX"/
C     ..
C     Executable statements 
C
      IF(JOB.EQ.0)THEN
         call MatGetInfo(A,MAT_LOCAL,info,IFAIL)
         nz_a = info(MAT_INFO_NZ_ALLOCATED)
         maxcols = info(MAT_INFO_COLUMNS_LOCAL)
         CALL VecGetLocalSize(rhs,n,IFAIL)
         maxcols = max(maxcols,n)
         write(6,*)'PE # ',MY_PE,nz_a,maxcols,n
         RETURN
      ENDIF
C
      IF( ITER. NE. 3 )RETURN
C
      WRITE(matfile(4:6),FMT="(I3.3)")0
      WRITE(rhsfile(4:6),FMT="(I3.3)")0
      WRITE(solfile(4:6),FMT="(I3.3)")0
      CALL MatGetOwnershipRange(A,rstart,rend,IFAIL)
      CALL HYPRE_IJMatrixCreate(PETSC_COMM_WORLD,rstart,rend-1,rstart,
     &rend-1,ij,IFAIL)
!     write(6,*)'HYPRE_IJMatrixCreate has returned ifail = ',ifail
      CALL HYPRE_IJMatrixSetObjectType(ij,HYPRE_PARCSR,IFAIL) 
!     write(6,*)'HYPRE_IJMatrixSetObjectType has returned ifail = ',
!    &ifail
      CALL  HYPRE_IJMatrixInitialize(ij,IFAIL)
!     write(6,*)'HYPRE_IJMatrixInitialize has returned ifail = ',
!    &ifail
      do i= rstart, rend-1
        CALL MatGetRow(A,i,ncols,cols,values,IFAIL)
        if(ncols.GT.maxcols)STOP 'Must increase maxcols'
        CALL HYPRE_IJMatrixSetValues(ij,1,ncols,i,cols,values,IFAIL)
        CALL MatRestoreRow(A,i,ncols,cols,values,IFAIL)
      enddo
      CALL HYPRE_IJMatrixAssemble(ij,IFAIL)
      CALL HYPRE_IJMatrixPrint(ij,matfile,IFAIL)
!     write(6,*)'HYPRE_IJMatrixPrint has returned ifail = ',
!    &ifail
      CALL HYPRE_IJMatrixDestroy(ij,IFAIL)
!     write(6,*)'HYPRE_IJMatrixDestroy has returned ifail = ',
!    &ifail
C
C     now dump rhs
C
      CALL VecGetOwnerShipRange(rhs,rstart,rend,IFAIL)
      CALL HYPRE_IJVectorCreate(PETSC_COMM_WORLD,rstart,rend-1,
     &ij,IFAIL)
!     write(6,*)'HYPRE_IJVectorCreate has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorSetObjectType(ij,HYPRE_PARCSR,IFAIL)
!     write(6,*)'HYPRE_IJSetObjectType has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorInitialize(ij,IFAIL)
!     write(6,*)'HYPRE_IJVectorInitialize has returned ifail = ',
!    &ifail
      CALL VecGetLocalSize(rhs,n,IFAIL)
      if(n.GT.maxcols)STOP 'Must increase maxcols'
      CALL VecGetArray(rhs,x_array,i_x,IFAIL)
      k = 0 
      do i = rstart,rend-1
         k = k + 1
         cols(k) = i
      enddo
!     write(6,*)rend-rstart,n
      CALL HYPRE_IJVectorSetValues(ij,n,cols,x_array(i_x + 1),IFAIL)
!     write(6,*)'HYPRE_IJVectorSetValues has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorAssemble(ij,IFAIL)
      CALL VecRestoreArray(rhs,x_array,i_x,IFAIL)
      CALL HYPRE_IJVectorPrint(ij,rhsfile,IFAIL)
!     write(6,*)'HYPRE_IJVectorPrint has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorDestroy(ij,IFAIL)
c
c     dump solution
c
      CALL VecGetOwnerShipRange(x,rstart,rend,IFAIL)
      CALL HYPRE_IJVectorCreate(PETSC_COMM_WORLD,rstart,rend-1,
     &ij,IFAIL)
!     write(6,*)'HYPRE_IJVectorCreate has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorSetObjectType(ij,HYPRE_PARCSR,IFAIL)
!     write(6,*)'HYPRE_IJSetObjectType has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorInitialize(ij,IFAIL)
!     write(6,*)'HYPRE_IJVectorInitialize has returned ifail = ',
!    &ifail
      CALL VecGetLocalSize(x,n,IFAIL)
      if(n.GT.maxcols)STOP 'Must increase maxcols'
      CALL VecGetArray(x,x_array,i_x,IFAIL)
      k = 0 
      do i = rstart,rend-1
         k = k + 1
         cols(k) = i
      enddo
      CALL HYPRE_IJVectorSetValues(ij,n,cols,x_array(i_x + 1),IFAIL)
!     write(6,*)'HYPRE_IJVectorSetValues has returned ifail = ',
!    &ifail
      CALL VecRestoreArray(x,x_array,i_x,IFAIL)
      CALL HYPRE_IJVectorAssemble(ij,IFAIL)
      CALL HYPRE_IJVectorPrint(ij,solfile,IFAIL)
!     write(6,*)'HYPRE_IJVectorPrint has returned ifail = ',
!    &ifail
      CALL HYPRE_IJVectorDestroy(ij,IFAIL)
      RETURN
      END