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: }