what's the difference between PetscViewerASCIIOpen() and PetscViewerBinaryOpen()?
Yujie
recrusader at gmail.com
Wed Jan 23 20:39:35 CST 2008
thank you very much, Aldo.
On 1/23/08, Aldo Bonfiglioli <aldo.bonfiglioli at unibas.it> wrote:
>
> > 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080123/e0e0d7c3/attachment.htm>
More information about the petsc-users
mailing list