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