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:   PetscScalar    *array;
  8:   PetscInt       *dnnz,*onnz,*dnnzu,*onnzu;
  9:   PetscInt       cst,Nbs,mbs,nbs,rbs,cbs;
 10:   PetscInt       im,i,m,n,M,N,*rows,start;

 12:   MatGetOwnershipRange(oldmat,&start,NULL);
 13:   MatGetOwnershipRangeColumn(oldmat,&cst,NULL);
 14:   MatCreateVecs(oldmat,&in,&out);
 15:   MatGetLocalSize(oldmat,&m,&n);
 16:   MatGetSize(oldmat,&M,&N);
 17:   PetscMalloc1(m,&rows);
 18:   if (reuse != MAT_REUSE_MATRIX) {
 19:     MatCreate(PetscObjectComm((PetscObject)oldmat),&mat);
 20:     MatSetSizes(mat,m,n,M,N);
 21:     MatSetType(mat,newtype);
 22:     MatSetBlockSizesFromMats(mat,oldmat,oldmat);
 23:     MatGetBlockSizes(mat,&rbs,&cbs);
 24:     mbs  = m/rbs;
 25:     nbs  = n/cbs;
 26:     Nbs  = N/cbs;
 27:     cst  = cst/cbs;
 28:     PetscMalloc4(mbs,&dnnz,mbs,&onnz,mbs,&dnnzu,mbs,&onnzu);
 29:     for (i=0; i<mbs; i++) {
 30:       dnnz[i]  = nbs;
 31:       onnz[i]  = Nbs - nbs;
 32:       dnnzu[i] = PetscMax(nbs - i,0);
 33:       onnzu[i] = PetscMax(Nbs - (cst + nbs),0);
 34:     }
 35:     MatXAIJSetPreallocation(mat,PETSC_DECIDE,dnnz,onnz,dnnzu,onnzu);
 36:     PetscFree4(dnnz,onnz,dnnzu,onnzu);
 37:     VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
 38:     MatSetUp(mat);
 39:   } else {
 40:     mat = *newmat;
 41:     MatZeroEntries(mat);
 42:   }
 43:   for (i=0; i<N; i++) {
 44:     PetscInt j;

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

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

 78:   MatShellGetContext(A,&B);
 80:   MatGetDiagonal(B,X);
 81:   return 0;
 82: }

 84: static PetscErrorCode MatMult_CF(Mat A,Vec X,Vec Y)
 85: {
 86:   Mat            B;

 88:   MatShellGetContext(A,&B);
 90:   MatMult(B,X,Y);
 91:   return 0;
 92: }

 94: static PetscErrorCode MatMultTranspose_CF(Mat A,Vec X,Vec Y)
 95: {
 96:   Mat            B;

 98:   MatShellGetContext(A,&B);
100:   MatMultTranspose(B,X,Y);
101:   return 0;
102: }

104: static PetscErrorCode MatDestroy_CF(Mat A)
105: {
106:   Mat            B;

108:   MatShellGetContext(A,&B);
110:   MatDestroy(&B);
111:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_anytype_C",NULL);
112:   return 0;
113: }

115: typedef struct {
116:   void           *userdata;
117:   PetscErrorCode (*userdestroy)(void*);
118:   PetscErrorCode (*numeric)(Mat);
119:   MatProductType ptype;
120:   Mat            Dwork;
121: } MatMatCF;

123: static PetscErrorCode MatProductDestroy_CF(void *data)
124: {
125:   MatMatCF       *mmcfdata = (MatMatCF*)data;

127:   if (mmcfdata->userdestroy) {
128:     (*mmcfdata->userdestroy)(mmcfdata->userdata);
129:   }
130:   MatDestroy(&mmcfdata->Dwork);
131:   PetscFree(mmcfdata);
132:   return 0;
133: }

135: static PetscErrorCode MatProductNumericPhase_CF(Mat A, Mat B, Mat C, void *data)
136: {
137:   MatMatCF       *mmcfdata = (MatMatCF*)data;

141:   /* the MATSHELL interface allows us to play with the product data */
142:   PetscNew(&C->product);
143:   C->product->type  = mmcfdata->ptype;
144:   C->product->data  = mmcfdata->userdata;
145:   C->product->Dwork = mmcfdata->Dwork;
146:   MatShellGetContext(A,&C->product->A);
147:   C->product->B = B;
148:   (*mmcfdata->numeric)(C);
149:   PetscFree(C->product);
150:   return 0;
151: }

153: static PetscErrorCode MatProductSymbolicPhase_CF(Mat A, Mat B, Mat C, void **data)
154: {
155:   MatMatCF       *mmcfdata;

157:   MatShellGetContext(A,&C->product->A);
158:   MatProductSetFromOptions(C);
159:   MatProductSymbolic(C);
160:   /* the MATSHELL interface does not allow non-empty product data */
161:   PetscNew(&mmcfdata);

163:   mmcfdata->numeric     = C->ops->productnumeric;
164:   mmcfdata->ptype       = C->product->type;
165:   mmcfdata->userdata    = C->product->data;
166:   mmcfdata->userdestroy = C->product->destroy;
167:   mmcfdata->Dwork       = C->product->Dwork;

169:   C->product->Dwork   = NULL;
170:   C->product->data    = NULL;
171:   C->product->destroy = NULL;
172:   C->product->A       = A;

174:   *data = mmcfdata;
175:   return 0;
176: }

178: /* only for A of type shell, mainly used for MatMat operations of shells with AXPYs */
179: static PetscErrorCode MatProductSetFromOptions_CF(Mat D)
180: {
181:   Mat            A,B,Ain;
182:   void           (*Af)(void) = NULL;
183:   PetscBool      flg;

185:   MatCheckProduct(D,1);
186:   if (D->product->type == MATPRODUCT_ABC) return 0;
187:   A = D->product->A;
188:   B = D->product->B;
189:   MatIsShell(A,&flg);
190:   if (!flg) return 0;
191:   PetscObjectQueryFunction((PetscObject)A,"MatProductSetFromOptions_anytype_C",&Af);
192:   if (Af == (void(*)(void))MatProductSetFromOptions_CF) {
193:     MatShellGetContext(A,&Ain);
194:   } else return 0;
195:   D->product->A = Ain;
196:   MatProductSetFromOptions(D);
197:   D->product->A = A;
198:   if (D->ops->productsymbolic) { /* we have a symbolic match, now populate the MATSHELL operations */
199:     MatShellSetMatProductOperation(A,D->product->type,MatProductSymbolicPhase_CF,MatProductNumericPhase_CF,MatProductDestroy_CF,((PetscObject)B)->type_name,NULL);
200:     MatProductSetFromOptions(D);
201:   }
202:   return 0;
203: }

205: PetscErrorCode MatConvertFrom_Shell(Mat A,MatType newtype,MatReuse reuse,Mat *B)
206: {
207:   Mat            M;
208:   PetscBool      flg;

210:   PetscStrcmp(newtype,MATSHELL,&flg);
212:   if (reuse == MAT_INITIAL_MATRIX) {
213:     PetscObjectReference((PetscObject)A);
214:     MatCreateShell(PetscObjectComm((PetscObject)A),A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,A,&M);
215:     MatSetBlockSizesFromMats(M,A,A);
216:     MatShellSetOperation(M,MATOP_MULT,          (void (*)(void))MatMult_CF);
217:     MatShellSetOperation(M,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_CF);
218:     MatShellSetOperation(M,MATOP_GET_DIAGONAL,  (void (*)(void))MatGetDiagonal_CF);
219:     MatShellSetOperation(M,MATOP_DESTROY,       (void (*)(void))MatDestroy_CF);
220:     PetscObjectComposeFunction((PetscObject)M,"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_CF);
221:     PetscFree(M->defaultvectype);
222:     PetscStrallocpy(A->defaultvectype,&M->defaultvectype);
223: #if defined(PETSC_HAVE_DEVICE)
224:     MatBindToCPU(M,A->boundtocpu);
225: #endif
226:     *B = M;
227:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
228:   return 0;
229: }