Actual source code: ex211.c


  2: static char help[] = "Tests MatCreateSubmatrix() in parallel.";

  4: #include <petscmat.h>
  5: #include <../src/mat/impls/aij/mpi/mpiaij.h>

  7: PetscErrorCode ISGetSeqIS_SameColDist_Private(Mat mat,IS isrow,IS iscol,IS *isrow_d,IS *iscol_d,IS *iscol_o,const PetscInt *garray[])
  8: {
  9:   Vec            x,cmap;
 10:   const PetscInt *is_idx;
 11:   PetscScalar    *xarray,*cmaparray;
 12:   PetscInt       ncols,isstart,*idx,m,rstart,count;
 13:   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)mat->data;
 14:   Mat            B=a->B;
 15:   Vec            lvec=a->lvec,lcmap;
 16:   PetscInt       i,cstart,cend,Bn=B->cmap->N;
 17:   MPI_Comm       comm;
 18:   PetscMPIInt    rank;
 19:   VecScatter     Mvctx;

 21:   PetscObjectGetComm((PetscObject)mat,&comm);
 22:   MPI_Comm_rank(comm,&rank);
 23:   ISGetLocalSize(iscol,&ncols);

 25:   /* (1) iscol is a sub-column vector of mat, pad it with '-1.' to form a full vector x */
 26:   MatCreateVecs(mat,&x,NULL);
 27:   VecDuplicate(x,&cmap);
 28:   VecSet(x,-1.0);
 29:   VecSet(cmap,-1.0);

 31:   VecDuplicate(lvec,&lcmap);

 33:   /* Get start indices */
 34:   MPI_Scan(&ncols,&isstart,1,MPIU_INT,MPI_SUM,comm);
 35:   isstart -= ncols;
 36:   MatGetOwnershipRangeColumn(mat,&cstart,&cend);

 38:   ISGetIndices(iscol,&is_idx);
 39:   VecGetArray(x,&xarray);
 40:   VecGetArray(cmap,&cmaparray);
 41:   PetscMalloc1(ncols,&idx);
 42:   for (i=0; i<ncols; i++) {
 43:     xarray[is_idx[i]-cstart]    = (PetscScalar)is_idx[i];
 44:     cmaparray[is_idx[i]-cstart] = (PetscScalar)(i + isstart);      /* global index of iscol[i] */
 45:     idx[i]                      = is_idx[i]-cstart; /* local index of iscol[i]  */
 46:   }
 47:   VecRestoreArray(x,&xarray);
 48:   VecRestoreArray(cmap,&cmaparray);
 49:   ISRestoreIndices(iscol,&is_idx);

 51:   /* Get iscol_d */
 52:   ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,iscol_d);
 53:   ISGetBlockSize(iscol,&i);
 54:   ISSetBlockSize(*iscol_d,i);

 56:   /* Get isrow_d */
 57:   ISGetLocalSize(isrow,&m);
 58:   rstart = mat->rmap->rstart;
 59:   PetscMalloc1(m,&idx);
 60:   ISGetIndices(isrow,&is_idx);
 61:   for (i=0; i<m; i++) idx[i] = is_idx[i]-rstart;
 62:   ISRestoreIndices(isrow,&is_idx);

 64:   ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,isrow_d);
 65:   ISGetBlockSize(isrow,&i);
 66:   ISSetBlockSize(*isrow_d,i);

 68:   /* (2) Scatter x and cmap using aij->Mvctx to get their off-process portions (see MatMult_MPIAIJ) */
 69: #if 0
 70:   if (!a->Mvctx_mpi1) {
 71:     /* a->Mvctx causes random 'count' in o-build? See src/mat/tests/runex59_2 */
 72:     a->Mvctx_mpi1_flg = PETSC_TRUE;
 73:     MatSetUpMultiply_MPIAIJ(mat);
 74:   }
 75:   Mvctx = a->Mvctx_mpi1;
 76: #endif
 77:   Mvctx = a->Mvctx;
 78:   VecScatterBegin(Mvctx,x,lvec,INSERT_VALUES,SCATTER_FORWARD);
 79:   VecScatterEnd(Mvctx,x,lvec,INSERT_VALUES,SCATTER_FORWARD);

 81:   VecScatterBegin(Mvctx,cmap,lcmap,INSERT_VALUES,SCATTER_FORWARD);
 82:   VecScatterEnd(Mvctx,cmap,lcmap,INSERT_VALUES,SCATTER_FORWARD);

 84:   /* (3) create sequential iscol_o (a subset of iscol) and isgarray */
 85:   /* off-process column indices */
 86:   count = 0;
 87:   PetscInt *cmap1;
 88:   PetscMalloc1(Bn,&idx);
 89:   PetscMalloc1(Bn,&cmap1);

 91:   VecGetArray(lvec,&xarray);
 92:   VecGetArray(lcmap,&cmaparray);
 93:   for (i=0; i<Bn; i++) {
 94:     if (PetscRealPart(xarray[i]) > -1.0) {
 95:       idx[count]   = i;                   /* local column index in off-diagonal part B */
 96:       cmap1[count] = (PetscInt)(PetscRealPart(cmaparray[i]));  /* column index in submat */
 97:       count++;
 98:     }
 99:   }
100:   printf("[%d] Bn %d, count %d\n",rank,Bn,count);
101:   VecRestoreArray(lvec,&xarray);
102:   VecRestoreArray(lcmap,&cmaparray);
103:   if (count != 6) {
104:     printf("[%d] count %d != 6 lvec:\n",rank,count);
105:     VecView(lvec,0);

107:     printf("[%d] count %d != 6 lcmap:\n",rank,count);
108:     VecView(lcmap,0);
109:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"count %d != 6",count);
110:   }

112:   ISCreateGeneral(PETSC_COMM_SELF,count,idx,PETSC_COPY_VALUES,iscol_o);
113:   /* cannot ensure iscol_o has same blocksize as iscol! */

115:   PetscFree(idx);

117:   *garray = cmap1;

119:   VecDestroy(&x);
120:   VecDestroy(&cmap);
121:   VecDestroy(&lcmap);
122:   return 0;
123: }

125: int main(int argc,char **args)
126: {
127:   Mat            C,A;
128:   PetscInt       i,j,m = 3,n = 2,rstart,rend;
129:   PetscMPIInt    size,rank;
130:   PetscScalar    v;
131:   IS             isrow,iscol;

133:   PetscInitialize(&argc,&args,(char*)0,help);
134:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
135:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
136:   n    = 2*size;

138:   MatCreate(PETSC_COMM_WORLD,&C);
139:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
140:   MatSetFromOptions(C);
141:   MatSetUp(C);

143:   /*
144:         This is JUST to generate a nice test matrix, all processors fill up
145:     the entire matrix. This is not something one would ever do in practice.
146:   */
147:   MatGetOwnershipRange(C,&rstart,&rend);
148:   for (i=rstart; i<rend; i++) {
149:     for (j=0; j<m*n; j++) {
150:       v    = i + j + 1;
151:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
152:     }
153:   }

155:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
156:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

158:   /*
159:      Generate a new matrix consisting of every second row and column of
160:    the original matrix
161:   */
162:   MatGetOwnershipRange(C,&rstart,&rend);
163:   /* Create parallel IS with the rows we want on THIS processor */
164:   ISCreateStride(PETSC_COMM_WORLD,(rend-rstart)/2,rstart,2,&isrow);
165:   /* Create parallel IS with the rows we want on THIS processor (same as rows for now) */
166:   ISCreateStride(PETSC_COMM_WORLD,(rend-rstart)/2,rstart,2,&iscol);

168:   IS             iscol_d,isrow_d,iscol_o;
169:   const PetscInt *garray;
170:   ISGetSeqIS_SameColDist_Private(C,isrow,iscol,&isrow_d,&iscol_d,&iscol_o,&garray);

172:   ISDestroy(&isrow_d);
173:   ISDestroy(&iscol_d);
174:   ISDestroy(&iscol_o);
175:   PetscFree(garray);

177:   MatCreateSubMatrix(C,isrow,iscol,MAT_INITIAL_MATRIX,&A);
178:   MatCreateSubMatrix(C,isrow,iscol,MAT_REUSE_MATRIX,&A);

180:   ISDestroy(&isrow);
181:   ISDestroy(&iscol);
182:   MatDestroy(&A);
183:   MatDestroy(&C);
184:   PetscFinalize();
185:   return 0;
186: }