1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2 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; 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 ierr = PetscObjectGetComm((PetscObject)oldmat,&comm);CHKERRQ(ierr); 16 17 ierr = MatGetOwnershipRange(oldmat,&start,NULL);CHKERRQ(ierr); 18 ierr = MatGetOwnershipRangeColumn(oldmat,&cst,NULL);CHKERRQ(ierr); 19 ierr = MatCreateVecs(oldmat,&in,&out);CHKERRQ(ierr); 20 ierr = MatGetLocalSize(oldmat,&m,&n);CHKERRQ(ierr); 21 ierr = MatGetSize(oldmat,&M,&N);CHKERRQ(ierr); 22 ierr = PetscMalloc1(m,&rows);CHKERRQ(ierr); 23 24 ierr = MatCreate(comm,&mat);CHKERRQ(ierr); 25 ierr = MatSetSizes(mat,m,n,M,N);CHKERRQ(ierr); 26 ierr = MatSetType(mat,newtype);CHKERRQ(ierr); 27 ierr = MatSetBlockSizesFromMats(mat,oldmat,oldmat);CHKERRQ(ierr); 28 ierr = MatGetBlockSizes(mat,&rbs,&cbs);CHKERRQ(ierr); 29 mbs = m/rbs; 30 nbs = n/cbs; 31 Nbs = N/cbs; 32 cst = cst/cbs; 33 ierr = PetscMalloc4(mbs,&dnnz,mbs,&onnz,mbs,&dnnzu,mbs,&onnzu);CHKERRQ(ierr); 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 ierr = MatXAIJSetPreallocation(mat,PETSC_DECIDE,dnnz,onnz,dnnzu,onnzu);CHKERRQ(ierr); 41 ierr = PetscFree4(dnnz,onnz,dnnzu,onnzu);CHKERRQ(ierr); 42 ierr = VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 43 ierr = MatSetUp(mat);CHKERRQ(ierr); 44 for (i=0; i<N; i++) { 45 PetscInt j; 46 47 ierr = VecZeroEntries(in);CHKERRQ(ierr); 48 ierr = VecSetValue(in,i,1.,INSERT_VALUES);CHKERRQ(ierr); 49 ierr = VecAssemblyBegin(in);CHKERRQ(ierr); 50 ierr = VecAssemblyEnd(in);CHKERRQ(ierr); 51 ierr = MatMult(oldmat,in,out);CHKERRQ(ierr); 52 ierr = VecGetArray(out,&array);CHKERRQ(ierr); 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 ierr = MatSetValues(mat,im,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 60 ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); 61 } 62 ierr = PetscFree(rows);CHKERRQ(ierr); 63 ierr = VecDestroy(&in);CHKERRQ(ierr); 64 ierr = VecDestroy(&out);CHKERRQ(ierr); 65 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 66 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 67 if (reuse == MAT_INPLACE_MATRIX) { 68 ierr = MatHeaderReplace(oldmat,&mat);CHKERRQ(ierr); 69 } else { 70 *newmat = mat; 71 } 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode MatGetDiagonal_CF(Mat A,Vec X) 76 { 77 Mat B; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 82 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 83 ierr = MatGetDiagonal(B,X);CHKERRQ(ierr); 84 PetscFunctionReturn(0); 85 } 86 87 static PetscErrorCode MatMult_CF(Mat A,Vec X,Vec Y) 88 { 89 Mat B; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 94 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 95 ierr = MatMult(B,X,Y);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 static PetscErrorCode MatMultTranspose_CF(Mat A,Vec X,Vec Y) 100 { 101 Mat B; 102 PetscErrorCode ierr; 103 104 PetscFunctionBegin; 105 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 106 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 107 ierr = MatMultTranspose(B,X,Y);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 static PetscErrorCode MatDestroy_CF(Mat A) 112 { 113 Mat B; 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 118 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 119 ierr = MatDestroy(&B);CHKERRQ(ierr); 120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_anytype_C",NULL);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 typedef struct { 125 void *userdata; 126 PetscErrorCode (*userdestroy)(void*); 127 PetscErrorCode (*numeric)(Mat); 128 MatProductType ptype; 129 Mat Dwork; 130 } MatMatCF; 131 132 static PetscErrorCode MatProductDestroy_CF(void *data) 133 { 134 PetscErrorCode ierr; 135 MatMatCF *mmcfdata = (MatMatCF*)data; 136 137 PetscFunctionBegin; 138 if (mmcfdata->userdestroy) { 139 ierr = (*mmcfdata->userdestroy)(mmcfdata->userdata);CHKERRQ(ierr); 140 } 141 ierr = MatDestroy(&mmcfdata->Dwork);CHKERRQ(ierr); 142 ierr = PetscFree(mmcfdata);CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 146 static PetscErrorCode MatProductNumericPhase_CF(Mat A, Mat B, Mat C, void *data) 147 { 148 PetscErrorCode ierr; 149 MatMatCF *mmcfdata = (MatMatCF*)data; 150 151 PetscFunctionBegin; 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 ierr = PetscNew(&C->product);CHKERRQ(ierr); 156 C->product->type = mmcfdata->ptype; 157 C->product->data = mmcfdata->userdata; 158 C->product->Dwork = mmcfdata->Dwork; 159 ierr = MatShellGetContext(A,&C->product->A);CHKERRQ(ierr); 160 C->product->B = B; 161 ierr = (*mmcfdata->numeric)(C);CHKERRQ(ierr); 162 ierr = PetscFree(C->product);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 static PetscErrorCode MatProductSymbolicPhase_CF(Mat A, Mat B, Mat C, void **data) 167 { 168 PetscErrorCode ierr; 169 MatMatCF *mmcfdata; 170 171 PetscFunctionBegin; 172 ierr = MatShellGetContext(A,&C->product->A);CHKERRQ(ierr); 173 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 174 ierr = MatProductSymbolic(C);CHKERRQ(ierr); 175 /* the MATSHELL interface does not allow non-empty product data */ 176 ierr = PetscNew(&mmcfdata);CHKERRQ(ierr); 177 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; 183 184 C->product->Dwork = NULL; 185 C->product->data = NULL; 186 C->product->destroy = NULL; 187 C->product->A = A; 188 189 *data = mmcfdata; 190 PetscFunctionReturn(0); 191 } 192 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; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 MatCheckProduct(D,1); 203 if (D->product->type == MATPRODUCT_ABC) PetscFunctionReturn(0); 204 A = D->product->A; 205 B = D->product->B; 206 ierr = MatIsShell(A,&flg);CHKERRQ(ierr); 207 if (!flg) PetscFunctionReturn(0); 208 ierr = PetscObjectQueryFunction((PetscObject)A,"MatProductSetFromOptions_anytype_C",&Af);CHKERRQ(ierr); 209 if (Af == (void(*)(void))MatProductSetFromOptions_CF) { 210 ierr = MatShellGetContext(A,&Ain);CHKERRQ(ierr); 211 } else PetscFunctionReturn(0); 212 D->product->A = Ain; 213 ierr = MatProductSetFromOptions(D);CHKERRQ(ierr); 214 D->product->A = A; 215 if (D->ops->productsymbolic) { /* we have a symbolic match, now populate the MATSHELL operations */ 216 ierr = MatShellSetMatProductOperation(A,D->product->type,MatProductSymbolicPhase_CF,MatProductNumericPhase_CF,MatProductDestroy_CF,((PetscObject)B)->type_name,NULL);CHKERRQ(ierr); 217 ierr = MatProductSetFromOptions(D);CHKERRQ(ierr); 218 } 219 PetscFunctionReturn(0); 220 } 221 222 PetscErrorCode MatConvertFrom_Shell(Mat A, MatType newtype,MatReuse reuse,Mat *B) 223 { 224 Mat M; 225 PetscBool flg; 226 PetscErrorCode ierr; 227 228 PetscFunctionBegin; 229 ierr = PetscStrcmp(newtype,MATSHELL,&flg);CHKERRQ(ierr); 230 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only conversion to MATSHELL"); 231 if (reuse == MAT_INITIAL_MATRIX) { 232 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 233 ierr = MatCreateShell(PetscObjectComm((PetscObject)A),A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,A,&M);CHKERRQ(ierr); 234 ierr = MatSetBlockSizesFromMats(M,A,A);CHKERRQ(ierr); 235 ierr = MatShellSetOperation(M,MATOP_MULT, (void (*)(void))MatMult_CF);CHKERRQ(ierr); 236 ierr = MatShellSetOperation(M,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_CF);CHKERRQ(ierr); 237 ierr = MatShellSetOperation(M,MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_CF);CHKERRQ(ierr); 238 ierr = MatShellSetOperation(M,MATOP_DESTROY, (void (*)(void))MatDestroy_CF);CHKERRQ(ierr); 239 ierr = PetscObjectComposeFunction((PetscObject)M,"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_CF);CHKERRQ(ierr); 240 ierr = PetscFree(M->defaultvectype);CHKERRQ(ierr); 241 ierr = PetscStrallocpy(A->defaultvectype,&M->defaultvectype);CHKERRQ(ierr); 242 *B = M; 243 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented"); 244 PetscFunctionReturn(0); 245 } 246