Commits

Lisandro Dalcin committed 4e0b4df

Custom implementation of Mat{View|Load}

Comments (0)

Files changed (4)

demo/InputOutput.c

 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
-  PetscErrorCode  ierr;
+  MPI_Comm       comm;
+  IGA            iga,iga1;
+  PetscViewer    viewer;
+  PetscErrorCode ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  PetscInt  i;
-  PetscInt  dim = 3;
-  PetscInt  dof = 1;
-  PetscBool b[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};
-  PetscInt  N[3] = {16,16,16}; 
-  PetscInt  p[3] = { 2, 2, 2};
-  PetscInt  C[3] = {-1,-1,-1};
-  PetscInt  n0=3, n1=3, n2=3, n3=3;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","InputOutput Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dof","dofs/node",__FILE__,dof,&dof,PETSC_NULL);CHKERRQ(ierr);
-  n0 = n1 = n2 = n3 = dim;
-  ierr = PetscOptionsBoolArray("-periodic","periodicity",     __FILE__,b,&n0,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray ("-N","number of elements",     __FILE__,N,&n1,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray ("-p","polynomial order",       __FILE__,p,&n2,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray ("-C","global continuity order",__FILE__,C,&n3,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsEnd();CHKERRQ(ierr);
-  if (n0<3) b[2] = b[0]; if (n0<2) b[1] = b[0];
-  if (n1<3) N[2] = N[0]; if (n1<2) N[1] = N[0];
-  if (n2<3) p[2] = p[0]; if (n2<2) p[1] = p[0];
-  if (n3<3) C[2] = C[0]; if (n3<2) C[1] = C[0];
-  for (i=0; i<dim; i++)  if (C[i] ==-1) C[i] = p[i] - 1;
-
-  MPI_Comm comm = PETSC_COMM_WORLD;
-  
-  IGA iga;
+  comm = PETSC_COMM_WORLD;
   ierr = IGACreate(comm,&iga);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
-  ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
-  for (i=0; i<dim; i++) {
-    IGAAxis axis;
-    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
-    ierr = IGAAxisSetPeriodic(axis,b[i]);CHKERRQ(ierr);
-    ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
-    ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
+  {
+    PetscInt  i;
+    PetscInt  dim = 3;
+    PetscInt  dof = 1;
+    PetscBool b[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};
+    PetscInt  N[3] = {16,16,16}; 
+    PetscInt  p[3] = { 2, 2, 2};
+    PetscInt  C[3] = {-1,-1,-1};
+    PetscInt  n0=3, n1=3, n2=3, n3=3;
+    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","InputOutput Options","IGA");CHKERRQ(ierr);
+    ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsInt("-dof","dofs/node",__FILE__,dof,&dof,PETSC_NULL);CHKERRQ(ierr);
+    n0 = n1 = n2 = n3 = dim;
+    ierr = PetscOptionsBoolArray("-periodic","periodicity",     __FILE__,b,&n0,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray ("-N","number of elements",     __FILE__,N,&n1,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray ("-p","polynomial order",       __FILE__,p,&n2,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray ("-C","global continuity order",__FILE__,C,&n3,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsEnd();CHKERRQ(ierr);
+    if (n0<3) b[2] = b[0]; if (n0<2) b[1] = b[0];
+    if (n1<3) N[2] = N[0]; if (n1<2) N[1] = N[0];
+    if (n2<3) p[2] = p[0]; if (n2<2) p[1] = p[0];
+    if (n3<3) C[2] = C[0]; if (n3<2) C[1] = C[0];
+    for (i=0; i<dim; i++)  if (C[i] ==-1) C[i] = p[i] - 1;
+    ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
+    ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
+    for (i=0; i<dim; i++) {
+      IGAAxis axis;
+      ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+      ierr = IGAAxisSetPeriodic(axis,b[i]);CHKERRQ(ierr);
+      ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
+      ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
+    }
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
-
-  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-  ierr = IGAGetDof(iga,&dof);CHKERRQ(ierr);
   
   ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
   ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr); /* just for testing */
 
-  IGA iga1;
   ierr = IGACreate(comm,&iga1);CHKERRQ(ierr);
   ierr = IGARead(iga1,"iga.dat");CHKERRQ(ierr);
   ierr = IGASetUp(iga1);CHKERRQ(ierr);
   ierr = IGASetUp(iga1);CHKERRQ(ierr);
   ierr = IGADestroy(&iga1);CHKERRQ(ierr);
 
-  PetscViewer viewer;
   ierr = PetscViewerBinaryOpen(comm,"iga.dat",FILE_MODE_WRITE,&viewer);
   ierr = IGASave(iga,viewer);CHKERRQ(ierr);
   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
 
-  IGA iga2;
-  ierr = IGACreate(comm,&iga2);CHKERRQ(ierr);
+  ierr = IGACreate(comm,&iga1);CHKERRQ(ierr);
   ierr = PetscViewerBinaryOpen(comm,"iga.dat",FILE_MODE_READ,&viewer);
-  ierr = IGALoad(iga2,viewer);CHKERRQ(ierr);
+  ierr = IGALoad(iga1,viewer);CHKERRQ(ierr);
   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-  ierr = IGASetUp(iga2);CHKERRQ(ierr);
-  ierr = IGADestroy(&iga2);CHKERRQ(ierr);
+  ierr = IGASetUp(iga1);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga1);CHKERRQ(ierr);
 
-  Vec vec;
-  PetscInt size;
-  PetscScalar value;
-  ierr = IGACreateVec(iga,&vec);CHKERRQ(ierr);
-  ierr = VecSet(vec,1.0);CHKERRQ(ierr);
-  ierr = IGAWriteVec(iga,vec,"igavec.dat");CHKERRQ(ierr);
-  ierr = VecSet(vec,0.0);CHKERRQ(ierr);
-  ierr = IGAReadVec (iga,vec,"igavec.dat");CHKERRQ(ierr);
-  ierr = VecGetSize(vec,&size);CHKERRQ(ierr);
-  ierr = VecSum(vec,&value);CHKERRQ(ierr);
-  if ((PetscReal)size != PetscRealPart(value))
-    SETERRQ(comm,PETSC_ERR_PLIB,"Bad data in file");
-  ierr = VecDestroy(&vec);CHKERRQ(ierr);
+  {
+    Vec         vec;
+    PetscInt    size,bs;
+    PetscScalar value;
+    ierr = IGACreateVec(iga,&vec);CHKERRQ(ierr);
+    ierr = VecSet(vec,1.0);CHKERRQ(ierr);
+    ierr = IGAWriteVec(iga,vec,"igavec.dat");CHKERRQ(ierr);
+    ierr = VecSet(vec,0.0);CHKERRQ(ierr);
+    ierr = IGAReadVec (iga,vec,"igavec.dat");CHKERRQ(ierr);
+    ierr = VecGetSize(vec,&size);CHKERRQ(ierr);
+    ierr = VecSum(vec,&value);CHKERRQ(ierr);
+    if ((PetscReal)size != PetscRealPart(value))
+      SETERRQ(comm,PETSC_ERR_PLIB,"Bad data in file");
+    ierr = VecDestroy(&vec);CHKERRQ(ierr);
 
-  ierr = VecCreate(comm,&vec);CHKERRQ(ierr);
-  ierr = VecSetBlockSize(vec,dof);CHKERRQ(ierr);
-  ierr = PetscViewerBinaryOpen(comm,"igavec.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
-  ierr = VecLoad(vec,viewer);CHKERRQ(ierr);
-  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-  ierr = VecGetSize(vec,&size);CHKERRQ(ierr);
-  ierr = VecSum(vec,&value);CHKERRQ(ierr);
-  if ((PetscReal)size != PetscRealPart(value))
-    SETERRQ(comm,PETSC_ERR_PLIB,"Bad data in file");
-  ierr = VecDestroy(&vec);CHKERRQ(ierr);
+    ierr = IGAGetDof(iga,&bs);CHKERRQ(ierr);
+    ierr = VecCreate(comm,&vec);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(vec,bs);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryOpen(comm,"igavec.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
+    ierr = VecLoad(vec,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+    ierr = VecGetSize(vec,&size);CHKERRQ(ierr);
+    ierr = VecSum(vec,&value);CHKERRQ(ierr);
+    if ((PetscReal)size != PetscRealPart(value))
+      SETERRQ(comm,PETSC_ERR_PLIB,"Bad Vec data in file");
+    ierr = VecDestroy(&vec);CHKERRQ(ierr);
+  }
+
+  {
+    Mat         mat;
+    Vec         diag,diag2;
+    PetscInt    r,rstart,rend;
+    PetscScalar *d;
+    PetscBool   match = PETSC_FALSE;
+
+    ierr = IGACreateVec(iga,&diag);CHKERRQ(ierr);
+    ierr = VecGetOwnershipRange(diag,&rstart,&rend);CHKERRQ(ierr);
+    ierr = VecGetArray(diag,&d);CHKERRQ(ierr);
+    for (r=rstart; r<rend; r++) d[r-rstart] = (PetscReal)r;
+    ierr = VecRestoreArray(diag,&d);CHKERRQ(ierr);
+
+    ierr = IGACreateMat(iga,&mat);CHKERRQ(ierr);
+    ierr = MatDiagonalSet(mat,diag,INSERT_VALUES);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryOpen(comm,"igamat.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
+    ierr = MatView(mat,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+    ierr = MatDestroy(&mat);CHKERRQ(ierr);
+
+    ierr = IGACreateMat(iga,&mat);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryOpen(comm,"igamat.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
+    ierr = MatLoad(mat,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+    ierr = IGACreateVec(iga,&diag2);CHKERRQ(ierr);
+    ierr = MatGetDiagonal(mat,diag2);CHKERRQ(ierr);
+    ierr = MatDestroy(&mat);CHKERRQ(ierr);
+
+    ierr = VecEqual(diag,diag2,&match);;CHKERRQ(ierr);
+    if (!match) SETERRQ(comm,PETSC_ERR_PLIB,"Bad Mat data in file");
+    ierr = VecDestroy(&diag);CHKERRQ(ierr);
+    ierr = VecDestroy(&diag2);CHKERRQ(ierr);
+  }
 
   ierr = IGADestroy(&iga);CHKERRQ(ierr);
   ierr = PetscFinalize();CHKERRQ(ierr);
 runex0a_8:
 	-@${MPIEXEC} -n 8 ./InputOutput ${OPTS} -dim 3 -N 13,11,7 -p 3,2,1
 runex0a.rm:
-	-@${RM} -f iga.dat iga.dat.info igavec.dat igavec.dat.info
+	-@${RM} -f iga*.dat iga*.dat.info
 
 runex1a_1:
 	-@${MPIEXEC} -n 1 ./L2Proj1D ${OPTS}
 			 iga->dim,iga->dof,iga->node_sizes,
 			 iga->node_lstart,iga->node_lwidth,
 			 iga->node_gstart,iga->node_gwidth);CHKERRQ(ierr);
-    /* build the block application ordering */
+    /* build the scalar and block application orderings */
+    ierr = IGA_Grid_GetAO(grid,&iga->ao);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)iga->ao);CHKERRQ(ierr);
     ierr = IGA_Grid_GetAOBlock(grid,&iga->aob);CHKERRQ(ierr);
     ierr = PetscObjectReference((PetscObject)iga->aob);CHKERRQ(ierr);
     /* build the scalar and block local to global mappings */
 #include "private/matimpl.h"
 #endif
 
-EXTERN_C_BEGIN
-extern PetscErrorCode MatView_MPI_DA(Mat,PetscViewer);
-extern PetscErrorCode MatLoad_MPI_DA(Mat,PetscViewer);
-EXTERN_C_END
+#if PETSC_VERSION_LE(3,3,0)
+#undef MatType
+typedef const char* MatType;
+#endif         
 
-PETSC_STATIC_INLINE
-PetscInt Product(const PetscInt a[3]) { return a[0]*a[1]*a[2]; }
+PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat,Mat);
+static       PetscErrorCode MatView_MPI_IGA(Mat,PetscViewer);
+static       PetscErrorCode MatLoad_MPI_IGA(Mat,PetscViewer);
+
+#undef  __FUNCT__
+#define __FUNCT__ "MatView_MPI_IGA"
+static PetscErrorCode MatView_MPI_IGA(Mat A,PetscViewer viewer)
+{
+  PetscViewerFormat format;
+  MPI_Comm          comm;
+  IGA               iga;
+  Mat               Anatural;
+  PetscInt          rstart,rend;
+  IS                is;
+  const char        *prefix;
+  PetscErrorCode    ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
+  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);
+  ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
+  if (format == PETSC_VIEWER_ASCII_INFO       ) PetscFunctionReturn(0);
+  if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0);
+
+  ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
+  ierr = PetscObjectQuery((PetscObject)A,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
+  if (!iga) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Matrix not generated from a IGA");
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,0);
+
+  /* Map natural ordering to PETSc ordering and create IS */
+  ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
+  ierr = ISCreateStride(comm,rend-rstart,rstart,1,&is);CHKERRQ(ierr);
+  ierr = AOApplicationToPetscIS(iga->ao,is);CHKERRQ(ierr);
+
+  /* Do permutation and view matrix */
+  ierr = MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);CHKERRQ(ierr);
+  ierr = PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);CHKERRQ(ierr);
+  ierr = PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);CHKERRQ(ierr);
+  ierr = MatView(Anatural,viewer);CHKERRQ(ierr);
+
+  ierr = ISDestroy(&is);CHKERRQ(ierr);
+  ierr = MatDestroy(&Anatural);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "MatLoad_MPI_IGA"
+static PetscErrorCode MatLoad_MPI_IGA(Mat A,PetscViewer viewer)
+{
+  MPI_Comm       comm;
+  IGA            iga;
+  MatType        mtype;
+  PetscInt       m,n,M,N;
+  Mat            Anatural,Apetsc;
+  PetscInt       rstart,rend;
+  IS             is;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
+  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);
+
+  ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
+  ierr = PetscObjectQuery((PetscObject)A,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
+  if (!iga) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Matrix not generated from a IGA");
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,0);
+
+  /* Create and load the matrix in natural ordering */
+  ierr = MatGetType(A,&mtype);CHKERRQ(ierr);
+  ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
+  ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
+  ierr = MatCreate(comm,&Anatural);CHKERRQ(ierr);
+  ierr = MatSetType(Anatural,mtype);CHKERRQ(ierr);
+  ierr = MatSetSizes(Anatural,m,n,M,N);CHKERRQ(ierr);
+  ierr = MatLoad(Anatural,viewer);CHKERRQ(ierr);
+
+  /* Map PETSc ordering to natural ordering and create IS */
+  ierr = MatGetOwnershipRange(Anatural,&rstart,&rend);CHKERRQ(ierr);
+  ierr = ISCreateStride(comm,rend-rstart,rstart,1,&is);CHKERRQ(ierr);
+  ierr = AOPetscToApplicationIS(iga->ao,is);CHKERRQ(ierr);
+
+  /* Do permutation and replace header */
+  ierr = MatGetSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Apetsc);CHKERRQ(ierr);
+  ierr = MatHeaderReplace(A,Apetsc);CHKERRQ(ierr);
+
+  ierr = ISDestroy(&is);CHKERRQ(ierr);
+  ierr = MatDestroy(&Anatural);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
 
 PETSC_STATIC_INLINE
 PetscInt Index3D(const PetscInt start[3],const PetscInt shape[3],
 
 PETSC_STATIC_INLINE
 PetscInt ColumnIndices(IGA iga,const PetscInt start[3],const PetscInt shape[3],
-                       PetscInt iA,PetscInt jA,PetscInt kA,PetscInt *stencil)
+                       PetscInt iA,PetscInt jA,PetscInt kA,PetscInt stencil[])
 {
   PetscInt first[3] = {0,0,0};
   PetscInt last [3] = {0,0,0};
 #define __FUNCT__ "FilterLowerTriangular"
 PetscErrorCode FilterLowerTriangular(PetscInt row,PetscInt *cnt,PetscInt col[])
 {
-  PetscInt       i,n;
+  PetscInt i,n;
   PetscFunctionBegin;
   for (i=0,n=0; i<*cnt; i++)
     if (col[i] >= row)
   ierr = IGAGetDof(iga,&bs);CHKERRQ(ierr);
   mtype = iga->mattype;
 
-  sizes  = iga->node_sizes;
   lstart = iga->node_lstart;
   lwidth = iga->node_lwidth;
+  sizes  = iga->node_sizes;
   for (i=0; i<dim; i++) {
     PetscInt gfirst, first = lstart[i];
     PetscInt glast,  last  = lstart[i] + lwidth[i] - 1;
     ierr = IGA_Grid_Destroy(&grid);CHKERRQ(ierr);
   }
 
-  n = Product(lwidth);
-  N = Product(sizes);
-  maxnnz = 1;
-  for(i=0; i<dim; i++)
+  maxnnz = n = N = 1;
+  for(i=0; i<dim; i++) {
     maxnnz *= (2*iga->axis[i]->p + 1); /* XXX do better ? */
+    n *= lwidth[i];
+    N *= sizes[i];
+  }
 
   ierr = MatCreate(comm,&A);CHKERRQ(ierr);
 #if PETSC_VERSION_LE(3,2,0)
   ierr = MatSetLocalToGlobalMappingBlock(A,iga->lgmapb,iga->lgmapb);CHKERRQ(ierr);
 
   ierr = PetscObjectCompose((PetscObject)A,"IGA",(PetscObject)iga);CHKERRQ(ierr);
+  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
+  if (size > 1) { /* change viewer to display matrix in natural ordering */
+    ierr = MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MPI_IGA);CHKERRQ(ierr);
+    ierr = MatShellSetOperation(A,MATOP_LOAD,(void (*)(void))MatLoad_MPI_IGA);CHKERRQ(ierr);
+  }
+
   *mat = A;
 
   { /* XXX */
 #else
     ierr = MatSetDM(*mat,iga->node_dm);CHKERRQ(ierr);
 #endif
-    ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
-    if (size > 1) { /* change viewer to display matrix in natural ordering */
-      ierr = MatShellSetOperation(*mat,MATOP_VIEW,(void (*)(void))MatView_MPI_DA);CHKERRQ(ierr);
-      ierr = MatShellSetOperation(*mat,MATOP_LOAD,(void (*)(void))MatLoad_MPI_DA);CHKERRQ(ierr);
-    }
   } /* XXX */
+
   PetscFunctionReturn(0);
 }