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