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