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