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:   MatShellSetContext(A, NULL);
112:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_anytype_C", NULL);
113:   return 0;
114: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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