Actual source code: transm.c

      2: #include <petsc/private/matimpl.h>

      4: typedef struct {
      5:   Mat A;
      6: } Mat_Transpose;

      8: PetscErrorCode MatMult_Transpose(Mat N,Vec x,Vec y)
      9: {
     10:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

     14:   MatMultTranspose(Na->A,x,y);
     15:   return(0);
     16: }

     18: PetscErrorCode MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
     19: {
     20:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

     24:   MatMultTransposeAdd(Na->A,v1,v2,v3);
     25:   return(0);
     26: }

     28: PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y)
     29: {
     30:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

     34:   MatMult(Na->A,x,y);
     35:   return(0);
     36: }

     38: PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
     39: {
     40:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

     44:   MatMultAdd(Na->A,v1,v2,v3);
     45:   return(0);
     46: }

     48: PetscErrorCode MatDestroy_Transpose(Mat N)
     49: {
     50:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

     54:   MatDestroy(&Na->A);
     55:   PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
     56:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
     57:   PetscFree(N->data);
     58:   return(0);
     59: }

     61: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
     62: {
     63:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

     67:   if (op == MAT_COPY_VALUES) {
     68:     MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);
     69:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
     70:     MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
     71:     MatTranspose(*m,MAT_INPLACE_MATRIX,m);
     72:   } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
     73:   return(0);
     74: }

     76: PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
     77: {
     78:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

     82:   MatCreateVecs(Aa->A,l,r);
     83:   return(0);
     84: }

     86: PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)
     87: {
     88:   Mat_Transpose  *Ya = (Mat_Transpose*)Y->data;
     89:   Mat_Transpose  *Xa = (Mat_Transpose*)X->data;
     90:   Mat              M = Ya->A;
     91:   Mat              N = Xa->A;

     95:   MatAXPY(M,a,N,str);
     96:   return(0);
     97: }

     99: PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
    100: {
    101:   Mat_Transpose  *X = (Mat_Transpose*)mat->data;

    105:   *has = PETSC_FALSE;
    106:   if (op == MATOP_MULT) {
    107:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);
    108:   } else if (op == MATOP_MULT_TRANSPOSE) {
    109:     MatHasOperation(X->A,MATOP_MULT,has);
    110:   } else if (op == MATOP_MULT_ADD) {
    111:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);
    112:   } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
    113:     MatHasOperation(X->A,MATOP_MULT_ADD,has);
    114:   } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
    115:   return(0);
    116: }

    118: /* used by hermitian transpose */
    119: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
    120: {
    121:   Mat            A,B,C,Ain,Bin,Cin;
    122:   PetscBool      Aistrans,Bistrans,Cistrans;
    123:   PetscInt       Atrans,Btrans,Ctrans;
    124:   MatProductType ptype;

    128:   MatCheckProduct(D,1);
    129:   A = D->product->A;
    130:   B = D->product->B;
    131:   C = D->product->C;
    132:   PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);
    133:   PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);
    134:   PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);
    135:   if (!Aistrans && !Bistrans && !Cistrans) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"This should not happen");
    136:   Atrans = 0;
    137:   Ain    = A;
    138:   while (Aistrans) {
    139:     Atrans++;
    140:     MatTransposeGetMat(Ain,&Ain);
    141:     PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);
    142:   }
    143:   Btrans = 0;
    144:   Bin    = B;
    145:   while (Bistrans) {
    146:     Btrans++;
    147:     MatTransposeGetMat(Bin,&Bin);
    148:     PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);
    149:   }
    150:   Ctrans = 0;
    151:   Cin    = C;
    152:   while (Cistrans) {
    153:     Ctrans++;
    154:     MatTransposeGetMat(Cin,&Cin);
    155:     PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&Cistrans);
    156:   }
    157:   Atrans = Atrans%2;
    158:   Btrans = Btrans%2;
    159:   Ctrans = Ctrans%2;
    160:   ptype = D->product->type; /* same product type by default */
    161:   if (Ain->symmetric) Atrans = 0;
    162:   if (Bin->symmetric) Btrans = 0;
    163:   if (Cin && Cin->symmetric) Ctrans = 0;

    165:   if (Atrans || Btrans || Ctrans) {
    166:     ptype = MATPRODUCT_UNSPECIFIED;
    167:     switch (D->product->type) {
    168:     case MATPRODUCT_AB:
    169:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
    170:         /* TODO custom implementation ? */
    171:       } else if (Atrans) { /* At * B */
    172:         ptype = MATPRODUCT_AtB;
    173:       } else { /* A * Bt */
    174:         ptype = MATPRODUCT_ABt;
    175:       }
    176:       break;
    177:     case MATPRODUCT_AtB:
    178:       if (Atrans && Btrans) { /* A * Bt */
    179:         ptype = MATPRODUCT_ABt;
    180:       } else if (Atrans) { /* A * B */
    181:         ptype = MATPRODUCT_AB;
    182:       } else { /* At * Bt we do not have support for this */
    183:         /* TODO custom implementation ? */
    184:       }
    185:       break;
    186:     case MATPRODUCT_ABt:
    187:       if (Atrans && Btrans) { /* At * B */
    188:         ptype = MATPRODUCT_AtB;
    189:       } else if (Atrans) { /* At * Bt we do not have support for this */
    190:         /* TODO custom implementation ? */
    191:       } else {  /* A * B */
    192:         ptype = MATPRODUCT_AB;
    193:       }
    194:       break;
    195:     case MATPRODUCT_PtAP:
    196:       if (Atrans) { /* PtAtP */
    197:         /* TODO custom implementation ? */
    198:       } else { /* RARt */
    199:         ptype = MATPRODUCT_RARt;
    200:       }
    201:       break;
    202:     case MATPRODUCT_RARt:
    203:       if (Atrans) { /* RAtRt */
    204:         /* TODO custom implementation ? */
    205:       } else { /* PtAP */
    206:         ptype = MATPRODUCT_PtAP;
    207:       }
    208:       break;
    209:     case MATPRODUCT_ABC:
    210:       /* TODO custom implementation ? */
    211:       break;
    212:     default: SETERRQ1(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[D->product->type]);
    213:     }
    214:   }
    215:   MatProductReplaceMats(Ain,Bin,Cin,D);
    216:   MatProductSetType(D,ptype);
    217:   MatProductSetFromOptions(D);
    218:   return(0);
    219: }

    221: PetscErrorCode MatGetDiagonal_Transpose(Mat A,Vec v)
    222: {
    223:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

    227:   MatGetDiagonal(Aa->A,v);
    228:   return(0);
    229: }

    231: PetscErrorCode MatConvert_Transpose(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
    232: {
    233:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;
    234:   Mat            B;

    238:   MatTranspose(Aa->A,MAT_INITIAL_MATRIX,&B);
    239:   if (reuse != MAT_INPLACE_MATRIX) {
    240:     MatConvert(B,newtype,reuse,newmat);
    241:     MatDestroy(&B);
    242:   } else {
    243:     MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
    244:     MatHeaderReplace(A,&B);
    245:   }
    246:   return(0);
    247: }

    249: PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
    250: {
    251:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

    254:   *M = Aa->A;
    255:   return(0);
    256: }

    258: /*@
    259:       MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT

    261:    Logically collective on Mat

    263:    Input Parameter:
    264: .   A  - the MATTRANSPOSE matrix

    266:    Output Parameter:
    267: .   M - the matrix object stored inside A

    269:    Level: intermediate

    271: .seealso: MatCreateTranspose()

    273: @*/
    274: PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
    275: {

    282:   PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));
    283:   return(0);
    284: }

    286: /*@
    287:       MatCreateTranspose - Creates a new matrix object that behaves like A'

    289:    Collective on Mat

    291:    Input Parameter:
    292: .   A  - the (possibly rectangular) matrix

    294:    Output Parameter:
    295: .   N - the matrix that represents A'

    297:    Level: intermediate

    299:    Notes:
    300:     The transpose A' is NOT actually formed! Rather the new matrix
    301:           object performs the matrix-vector product by using the MatMultTranspose() on
    302:           the original matrix

    304: .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()

    306: @*/
    307: PetscErrorCode  MatCreateTranspose(Mat A,Mat *N)
    308: {
    310:   PetscInt       m,n;
    311:   Mat_Transpose  *Na;
    312:   VecType        vtype;

    315:   MatGetLocalSize(A,&m,&n);
    316:   MatCreate(PetscObjectComm((PetscObject)A),N);
    317:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
    318:   PetscLayoutSetUp((*N)->rmap);
    319:   PetscLayoutSetUp((*N)->cmap);
    320:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

    322:   PetscNewLog(*N,&Na);
    323:   (*N)->data = (void*) Na;
    324:   PetscObjectReference((PetscObject)A);
    325:   Na->A      = A;

    327:   (*N)->ops->destroy               = MatDestroy_Transpose;
    328:   (*N)->ops->mult                  = MatMult_Transpose;
    329:   (*N)->ops->multadd               = MatMultAdd_Transpose;
    330:   (*N)->ops->multtranspose         = MatMultTranspose_Transpose;
    331:   (*N)->ops->multtransposeadd      = MatMultTransposeAdd_Transpose;
    332:   (*N)->ops->duplicate             = MatDuplicate_Transpose;
    333:   (*N)->ops->getvecs               = MatCreateVecs_Transpose;
    334:   (*N)->ops->axpy                  = MatAXPY_Transpose;
    335:   (*N)->ops->hasoperation          = MatHasOperation_Transpose;
    336:   (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
    337:   (*N)->ops->getdiagonal           = MatGetDiagonal_Transpose;
    338:   (*N)->ops->convert               = MatConvert_Transpose;
    339:   (*N)->assembled                  = PETSC_TRUE;

    341:   PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);
    342:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
    343:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
    344:   MatGetVecType(A,&vtype);
    345:   MatSetVecType(*N,vtype);
    346: #if defined(PETSC_HAVE_DEVICE)
    347:   MatBindToCPU(*N,A->boundtocpu);
    348: #endif
    349:   MatSetUp(*N);
    350:   return(0);
    351: }