Actual source code: shellcnv.c

  1: #include <petsc/private/matimpl.h>

  3: PetscErrorCode MatConvert_Shell(Mat oldmat, MatType newtype,MatReuse reuse,Mat *newmat)
  4: {
  5:   Mat            mat;
  6:   Vec            in,out;
  7:   MPI_Comm       comm;
  8:   PetscScalar    *array;
  9:   PetscInt       *dnnz,*onnz,*dnnzu,*onnzu;
 10:   PetscInt       cst,Nbs,mbs,nbs,rbs,cbs;
 11:   PetscInt       im,i,m,n,M,N,*rows,start;

 15:   PetscObjectGetComm((PetscObject)oldmat,&comm);

 17:   MatGetOwnershipRange(oldmat,&start,NULL);
 18:   MatGetOwnershipRangeColumn(oldmat,&cst,NULL);
 19:   MatCreateVecs(oldmat,&in,&out);
 20:   MatGetLocalSize(oldmat,&m,&n);
 21:   MatGetSize(oldmat,&M,&N);
 22:   PetscMalloc1(m,&rows);

 24:   MatCreate(comm,&mat);
 25:   MatSetSizes(mat,m,n,M,N);
 26:   MatSetType(mat,newtype);
 27:   MatSetBlockSizesFromMats(mat,oldmat,oldmat);
 28:   MatGetBlockSizes(mat,&rbs,&cbs);
 29:   mbs  = m/rbs;
 30:   nbs  = n/cbs;
 31:   Nbs  = N/cbs;
 32:   cst  = cst/cbs;
 33:   PetscMalloc4(mbs,&dnnz,mbs,&onnz,mbs,&dnnzu,mbs,&onnzu);
 34:   for (i=0; i<mbs; i++) {
 35:     dnnz[i]  = nbs;
 36:     onnz[i]  = Nbs - nbs;
 37:     dnnzu[i] = PetscMax(nbs - i,0);
 38:     onnzu[i] = PetscMax(Nbs - (cst + nbs),0);
 39:   }
 40:   MatXAIJSetPreallocation(mat,PETSC_DECIDE,dnnz,onnz,dnnzu,onnzu);
 41:   PetscFree4(dnnz,onnz,dnnzu,onnzu);
 42:   VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
 43:   MatSetUp(mat);
 44:   for (i=0; i<N; i++) {
 45:     PetscInt j;

 47:     VecZeroEntries(in);
 48:     VecSetValue(in,i,1.,INSERT_VALUES);
 49:     VecAssemblyBegin(in);
 50:     VecAssemblyEnd(in);
 51:     MatMult(oldmat,in,out);
 52:     VecGetArray(out,&array);
 53:     for (j=0, im = 0; j<m; j++) {
 54:       if (PetscAbsScalar(array[j]) == 0.0) continue;
 55:       rows[im]  = j+start;
 56:       array[im] = array[j];
 57:       im++;
 58:     }
 59:     MatSetValues(mat,im,rows,1,&i,array,INSERT_VALUES);
 60:     VecRestoreArray(out,&array);
 61:   }
 62:   PetscFree(rows);
 63:   VecDestroy(&in);
 64:   VecDestroy(&out);
 65:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
 67:   if (reuse == MAT_INPLACE_MATRIX) {
 68:     MatHeaderReplace(oldmat,&mat);
 69:   } else {
 70:     *newmat = mat;
 71:   }
 72:   return(0);
 73: }

 75: static PetscErrorCode MatGetDiagonal_CF(Mat A,Vec X)
 76: {
 77:   Mat            B;

 81:   MatShellGetContext(A,&B);
 82:   if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
 83:   MatGetDiagonal(B,X);
 84:   return(0);
 85: }

 87: static PetscErrorCode MatMult_CF(Mat A,Vec X,Vec Y)
 88: {
 89:   Mat            B;

 93:   MatShellGetContext(A,&B);
 94:   if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
 95:   MatMult(B,X,Y);
 96:   return(0);
 97: }

 99: static PetscErrorCode MatMultTranspose_CF(Mat A,Vec X,Vec Y)
100: {
101:   Mat            B;

105:   MatShellGetContext(A,&B);
106:   if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
107:   MatMultTranspose(B,X,Y);
108:   return(0);
109: }

111: static PetscErrorCode MatDestroy_CF(Mat A)
112: {
113:   Mat            B;

117:   MatShellGetContext(A,&B);
118:   if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
119:   MatDestroy(&B);
120:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_anytype_C",NULL);
121:   return(0);
122: }

124: typedef struct {
125:   void           *userdata;
126:   PetscErrorCode (*userdestroy)(void*);
127:   PetscErrorCode (*numeric)(Mat);
128:   MatProductType ptype;
129:   Mat            Dwork;
130: } MatMatCF;

132: static PetscErrorCode MatProductDestroy_CF(void *data)
133: {
135:   MatMatCF       *mmcfdata = (MatMatCF*)data;

138:   if (mmcfdata->userdestroy) {
139:     (*mmcfdata->userdestroy)(mmcfdata->userdata);
140:   }
141:   MatDestroy(&mmcfdata->Dwork);
142:   PetscFree(mmcfdata);
143:   return(0);
144: }

146: static PetscErrorCode MatProductNumericPhase_CF(Mat A, Mat B, Mat C, void *data)
147: {
149:   MatMatCF       *mmcfdata = (MatMatCF*)data;

152:   if (!mmcfdata) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing data");
153:   if (!mmcfdata->numeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric operation");
154:   /* the MATSHELL interface allows us to play with the product data */
155:   PetscNew(&C->product);
156:   C->product->type  = mmcfdata->ptype;
157:   C->product->data  = mmcfdata->userdata;
158:   C->product->Dwork = mmcfdata->Dwork;
159:   MatShellGetContext(A,&C->product->A);
160:   C->product->B = B;
161:   (*mmcfdata->numeric)(C);
162:   PetscFree(C->product);
163:   return(0);
164: }

166: static PetscErrorCode MatProductSymbolicPhase_CF(Mat A, Mat B, Mat C, void **data)
167: {
169:   MatMatCF       *mmcfdata;

172:   MatShellGetContext(A,&C->product->A);
173:   MatProductSetFromOptions(C);
174:   MatProductSymbolic(C);
175:   /* the MATSHELL interface does not allow non-empty product data */
176:   PetscNew(&mmcfdata);

178:   mmcfdata->numeric     = C->ops->productnumeric;
179:   mmcfdata->ptype       = C->product->type;
180:   mmcfdata->userdata    = C->product->data;
181:   mmcfdata->userdestroy = C->product->destroy;
182:   mmcfdata->Dwork       = C->product->Dwork;

184:   C->product->Dwork   = NULL;
185:   C->product->data    = NULL;
186:   C->product->destroy = NULL;
187:   C->product->A       = A;

189:   *data = mmcfdata;
190:   return(0);
191: }

193: /* only for A of type shell, mainly used for MatMat operations of shells with AXPYs */
194: static PetscErrorCode MatProductSetFromOptions_CF(Mat D)
195: {
196:   Mat            A,B,Ain;
197:   void           (*Af)(void) = NULL;
198:   PetscBool      flg;

202:   MatCheckProduct(D,1);
203:   if (D->product->type == MATPRODUCT_ABC) return(0);
204:   A = D->product->A;
205:   B = D->product->B;
206:   MatIsShell(A,&flg);
207:   if (!flg) return(0);
208:   PetscObjectQueryFunction((PetscObject)A,"MatProductSetFromOptions_anytype_C",&Af);
209:   if (Af == (void(*)(void))MatProductSetFromOptions_CF) {
210:     MatShellGetContext(A,&Ain);
211:   } else return(0);
212:   D->product->A = Ain;
213:   MatProductSetFromOptions(D);
214:   D->product->A = A;
215:   if (D->ops->productsymbolic) { /* we have a symbolic match, now populate the MATSHELL operations */
216:     MatShellSetMatProductOperation(A,D->product->type,MatProductSymbolicPhase_CF,MatProductNumericPhase_CF,MatProductDestroy_CF,((PetscObject)B)->type_name,NULL);
217:     MatProductSetFromOptions(D);
218:   }
219:   return(0);
220: }

222: PetscErrorCode MatConvertFrom_Shell(Mat A, MatType newtype,MatReuse reuse,Mat *B)
223: {
224:   Mat            M;
225:   PetscBool      flg;

229:   PetscStrcmp(newtype,MATSHELL,&flg);
230:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only conversion to MATSHELL");
231:   if (reuse == MAT_INITIAL_MATRIX) {
232:     PetscObjectReference((PetscObject)A);
233:     MatCreateShell(PetscObjectComm((PetscObject)A),A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,A,&M);
234:     MatSetBlockSizesFromMats(M,A,A);
235:     MatShellSetOperation(M,MATOP_MULT,          (void (*)(void))MatMult_CF);
236:     MatShellSetOperation(M,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_CF);
237:     MatShellSetOperation(M,MATOP_GET_DIAGONAL,  (void (*)(void))MatGetDiagonal_CF);
238:     MatShellSetOperation(M,MATOP_DESTROY,       (void (*)(void))MatDestroy_CF);
239:     PetscObjectComposeFunction((PetscObject)M,"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_CF);
240:     PetscFree(M->defaultvectype);
241:     PetscStrallocpy(A->defaultvectype,&M->defaultvectype);
242: #if defined(PETSC_HAVE_DEVICE)
243:     MatBindToCPU(M,A->boundtocpu);
244: #endif
245:     *B = M;
246:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
247:   return(0);
248: }