Actual source code: transm.c

      1: #include <../src/mat/impls/shell/shell.h>

      3: static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y)
      4: {
      5:   Mat A;

      7:   PetscFunctionBegin;
      8:   PetscCall(MatShellGetContext(N, &A));
      9:   PetscCall(MatMultTranspose(A, x, y));
     10:   PetscFunctionReturn(PETSC_SUCCESS);
     11: }

     13: static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y)
     14: {
     15:   Mat A;

     17:   PetscFunctionBegin;
     18:   PetscCall(MatShellGetContext(N, &A));
     19:   PetscCall(MatMult(A, x, y));
     20:   PetscFunctionReturn(PETSC_SUCCESS);
     21: }

     23: static PetscErrorCode MatSolve_Transpose_LU(Mat N, Vec b, Vec x)
     24: {
     25:   Mat A;

     27:   PetscFunctionBegin;
     28:   PetscCall(MatShellGetContext(N, &A));
     29:   PetscCall(MatSolveTranspose(A, b, x));
     30:   PetscFunctionReturn(PETSC_SUCCESS);
     31: }

     33: static PetscErrorCode MatSolveAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x)
     34: {
     35:   Mat A;

     37:   PetscFunctionBegin;
     38:   PetscCall(MatShellGetContext(N, &A));
     39:   PetscCall(MatSolveTransposeAdd(A, b, y, x));
     40:   PetscFunctionReturn(PETSC_SUCCESS);
     41: }

     43: static PetscErrorCode MatSolveTranspose_Transpose_LU(Mat N, Vec b, Vec x)
     44: {
     45:   Mat A;

     47:   PetscFunctionBegin;
     48:   PetscCall(MatShellGetContext(N, &A));
     49:   PetscCall(MatSolve(A, b, x));
     50:   PetscFunctionReturn(PETSC_SUCCESS);
     51: }

     53: static PetscErrorCode MatSolveTransposeAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x)
     54: {
     55:   Mat A;

     57:   PetscFunctionBegin;
     58:   PetscCall(MatShellGetContext(N, &A));
     59:   PetscCall(MatSolveAdd(A, b, y, x));
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: static PetscErrorCode MatMatSolve_Transpose_LU(Mat N, Mat B, Mat X)
     64: {
     65:   Mat A;

     67:   PetscFunctionBegin;
     68:   PetscCall(MatShellGetContext(N, &A));
     69:   PetscCall(MatMatSolveTranspose(A, B, X));
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: static PetscErrorCode MatMatSolveTranspose_Transpose_LU(Mat N, Mat B, Mat X)
     74: {
     75:   Mat A;

     77:   PetscFunctionBegin;
     78:   PetscCall(MatShellGetContext(N, &A));
     79:   PetscCall(MatMatSolve(A, B, X));
     80:   PetscFunctionReturn(PETSC_SUCCESS);
     81: }

     83: static PetscErrorCode MatLUFactor_Transpose(Mat N, IS row, IS col, const MatFactorInfo *minfo)
     84: {
     85:   Mat A;

     87:   PetscFunctionBegin;
     88:   PetscCall(MatShellGetContext(N, &A));
     89:   PetscCall(MatLUFactor(A, col, row, minfo));
     90:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_LU));
     91:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_LU));
     92:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_LU));
     93:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_LU));
     94:   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_LU));
     95:   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_LU));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode MatSolve_Transpose_Cholesky(Mat N, Vec b, Vec x)
    100: {
    101:   Mat A;

    103:   PetscFunctionBegin;
    104:   PetscCall(MatShellGetContext(N, &A));
    105:   PetscCall(MatSolveTranspose(A, b, x));
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: static PetscErrorCode MatSolveAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x)
    110: {
    111:   Mat A;

    113:   PetscFunctionBegin;
    114:   PetscCall(MatShellGetContext(N, &A));
    115:   PetscCall(MatSolveTransposeAdd(A, b, y, x));
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: static PetscErrorCode MatSolveTranspose_Transpose_Cholesky(Mat N, Vec b, Vec x)
    120: {
    121:   Mat A;

    123:   PetscFunctionBegin;
    124:   PetscCall(MatShellGetContext(N, &A));
    125:   PetscCall(MatSolve(A, b, x));
    126:   PetscFunctionReturn(PETSC_SUCCESS);
    127: }

    129: static PetscErrorCode MatSolveTransposeAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x)
    130: {
    131:   Mat A;

    133:   PetscFunctionBegin;
    134:   PetscCall(MatShellGetContext(N, &A));
    135:   PetscCall(MatSolveAdd(A, b, y, x));
    136:   PetscFunctionReturn(PETSC_SUCCESS);
    137: }

    139: static PetscErrorCode MatMatSolve_Transpose_Cholesky(Mat N, Mat B, Mat X)
    140: {
    141:   Mat A;

    143:   PetscFunctionBegin;
    144:   PetscCall(MatShellGetContext(N, &A));
    145:   PetscCall(MatMatSolveTranspose(A, B, X));
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: static PetscErrorCode MatMatSolveTranspose_Transpose_Cholesky(Mat N, Mat B, Mat X)
    150: {
    151:   Mat A;

    153:   PetscFunctionBegin;
    154:   PetscCall(MatShellGetContext(N, &A));
    155:   PetscCall(MatMatSolve(A, B, X));
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: static PetscErrorCode MatCholeskyFactor_Transpose(Mat N, IS perm, const MatFactorInfo *minfo)
    160: {
    161:   Mat A;

    163:   PetscFunctionBegin;
    164:   PetscCall(MatShellGetContext(N, &A));
    165:   PetscCall(MatCholeskyFactor(A, perm, minfo));
    166:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_Cholesky));
    167:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_Cholesky));
    168:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_Cholesky));
    169:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_Cholesky));
    170:   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_Cholesky));
    171:   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_Cholesky));
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }

    175: static PetscErrorCode MatLUFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info)
    176: {
    177:   Mat A, FA;

    179:   PetscFunctionBegin;
    180:   PetscCall(MatShellGetContext(N, &A));
    181:   PetscCall(MatShellGetContext(F, &FA));
    182:   PetscCall(MatLUFactorNumeric(FA, A, info));
    183:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_LU));
    184:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_LU));
    185:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_LU));
    186:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_LU));
    187:   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_LU));
    188:   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_LU));
    189:   PetscFunctionReturn(PETSC_SUCCESS);
    190: }

    192: static PetscErrorCode MatLUFactorSymbolic_Transpose(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info)
    193: {
    194:   Mat A, FA;

    196:   PetscFunctionBegin;
    197:   PetscCall(MatShellGetContext(N, &A));
    198:   PetscCall(MatShellGetContext(F, &FA));
    199:   PetscCall(MatLUFactorSymbolic(FA, A, row, col, info));
    200:   PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (PetscErrorCodeFn *)MatLUFactorNumeric_Transpose));
    201:   PetscFunctionReturn(PETSC_SUCCESS);
    202: }

    204: static PetscErrorCode MatCholeskyFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info)
    205: {
    206:   Mat A, FA;

    208:   PetscFunctionBegin;
    209:   PetscCall(MatShellGetContext(N, &A));
    210:   PetscCall(MatShellGetContext(F, &FA));
    211:   PetscCall(MatCholeskyFactorNumeric(FA, A, info));
    212:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_Cholesky));
    213:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_Cholesky));
    214:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_Cholesky));
    215:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_Cholesky));
    216:   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_Cholesky));
    217:   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_Cholesky));
    218:   PetscFunctionReturn(PETSC_SUCCESS);
    219: }

    221: static PetscErrorCode MatCholeskyFactorSymbolic_Transpose(Mat F, Mat N, IS perm, const MatFactorInfo *info)
    222: {
    223:   Mat A, FA;

    225:   PetscFunctionBegin;
    226:   PetscCall(MatShellGetContext(N, &A));
    227:   PetscCall(MatShellGetContext(F, &FA));
    228:   PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info));
    229:   PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (PetscErrorCodeFn *)MatCholeskyFactorNumeric_Transpose));
    230:   PetscFunctionReturn(PETSC_SUCCESS);
    231: }

    233: static PetscErrorCode MatGetFactor_Transpose(Mat N, MatSolverType type, MatFactorType ftype, Mat *F)
    234: {
    235:   Mat A, FA;

    237:   PetscFunctionBegin;
    238:   PetscCall(MatShellGetContext(N, &A));
    239:   PetscCall(MatGetFactor(A, type, ftype, &FA));
    240:   PetscCall(MatCreateTranspose(FA, F));
    241:   if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatLUFactorSymbolic_Transpose));
    242:   else if (ftype == MAT_FACTOR_CHOLESKY) {
    243:     PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatCholeskyFactorSymbolic_Transpose));
    244:     PetscCall(MatPropagateSymmetryOptions(A, FA));
    245:   } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]);
    246:   (*F)->factortype = ftype;
    247:   PetscCall(MatDestroy(&FA));
    248:   PetscFunctionReturn(PETSC_SUCCESS);
    249: }

    251: static PetscErrorCode MatDestroy_Transpose(Mat N)
    252: {
    253:   Mat A;

    255:   PetscFunctionBegin;
    256:   PetscCall(MatShellGetContext(N, &A));
    257:   PetscCall(MatDestroy(&A));
    258:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
    259:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
    260:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
    261:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL));
    262:   PetscFunctionReturn(PETSC_SUCCESS);
    263: }

    265: static PetscErrorCode MatGetInfo_Transpose(Mat N, MatInfoType flag, MatInfo *info)
    266: {
    267:   Mat A;

    269:   PetscFunctionBegin;
    270:   PetscCall(MatShellGetContext(N, &A));
    271:   PetscCall(MatGetInfo(A, flag, info));
    272:   PetscFunctionReturn(PETSC_SUCCESS);
    273: }

    275: static PetscErrorCode MatFactorGetSolverType_Transpose(Mat N, MatSolverType *type)
    276: {
    277:   Mat A;

    279:   PetscFunctionBegin;
    280:   PetscCall(MatShellGetContext(N, &A));
    281:   PetscCall(MatFactorGetSolverType(A, type));
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
    286: {
    287:   Mat A, C;

    289:   PetscFunctionBegin;
    290:   PetscCall(MatShellGetContext(N, &A));
    291:   PetscCall(MatDuplicate(A, op, &C));
    292:   PetscCall(MatCreateTranspose(C, m));
    293:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
    294:   PetscCall(MatDestroy(&C));
    295:   PetscFunctionReturn(PETSC_SUCCESS);
    296: }

    298: static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
    299: {
    300:   Mat A;

    302:   PetscFunctionBegin;
    303:   PetscCall(MatShellGetContext(mat, &A));
    304:   *has = PETSC_FALSE;
    305:   if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
    306:     PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
    307:   } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
    308:     PetscCall(MatHasOperation(A, MATOP_MULT, has));
    309:   } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: static PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
    314: {
    315:   Mat            A, B, C, Ain, Bin, Cin;
    316:   PetscBool      Aistrans, Bistrans, Cistrans;
    317:   PetscInt       Atrans, Btrans, Ctrans;
    318:   MatProductType ptype;

    320:   PetscFunctionBegin;
    321:   MatCheckProduct(D, 1);
    322:   A = D->product->A;
    323:   B = D->product->B;
    324:   C = D->product->C;
    325:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
    326:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
    327:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
    328:   PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
    329:   Atrans = 0;
    330:   Ain    = A;
    331:   while (Aistrans) {
    332:     Atrans++;
    333:     PetscCall(MatTransposeGetMat(Ain, &Ain));
    334:     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
    335:   }
    336:   Btrans = 0;
    337:   Bin    = B;
    338:   while (Bistrans) {
    339:     Btrans++;
    340:     PetscCall(MatTransposeGetMat(Bin, &Bin));
    341:     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
    342:   }
    343:   Ctrans = 0;
    344:   Cin    = C;
    345:   while (Cistrans) {
    346:     Ctrans++;
    347:     PetscCall(MatTransposeGetMat(Cin, &Cin));
    348:     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
    349:   }
    350:   Atrans = Atrans % 2;
    351:   Btrans = Btrans % 2;
    352:   Ctrans = Ctrans % 2;
    353:   ptype  = D->product->type; /* same product type by default */
    354:   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
    355:   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
    356:   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;

    358:   if (Atrans || Btrans || Ctrans) {
    359:     ptype = MATPRODUCT_UNSPECIFIED;
    360:     switch (D->product->type) {
    361:     case MATPRODUCT_AB:
    362:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
    363:         /* TODO custom implementation ? */
    364:       } else if (Atrans) { /* At * B */
    365:         ptype = MATPRODUCT_AtB;
    366:       } else { /* A * Bt */
    367:         ptype = MATPRODUCT_ABt;
    368:       }
    369:       break;
    370:     case MATPRODUCT_AtB:
    371:       if (Atrans && Btrans) { /* A * Bt */
    372:         ptype = MATPRODUCT_ABt;
    373:       } else if (Atrans) { /* A * B */
    374:         ptype = MATPRODUCT_AB;
    375:       } else { /* At * Bt we do not have support for this */
    376:         /* TODO custom implementation ? */
    377:       }
    378:       break;
    379:     case MATPRODUCT_ABt:
    380:       if (Atrans && Btrans) { /* At * B */
    381:         ptype = MATPRODUCT_AtB;
    382:       } else if (Atrans) { /* At * Bt we do not have support for this */
    383:         /* TODO custom implementation ? */
    384:       } else { /* A * B */
    385:         ptype = MATPRODUCT_AB;
    386:       }
    387:       break;
    388:     case MATPRODUCT_PtAP:
    389:       if (Atrans) { /* PtAtP */
    390:         /* TODO custom implementation ? */
    391:       } else { /* RARt */
    392:         ptype = MATPRODUCT_RARt;
    393:       }
    394:       break;
    395:     case MATPRODUCT_RARt:
    396:       if (Atrans) { /* RAtRt */
    397:         /* TODO custom implementation ? */
    398:       } else { /* PtAP */
    399:         ptype = MATPRODUCT_PtAP;
    400:       }
    401:       break;
    402:     case MATPRODUCT_ABC:
    403:       /* TODO custom implementation ? */
    404:       break;
    405:     default:
    406:       SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
    407:     }
    408:   }
    409:   PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
    410:   PetscCall(MatProductSetType(D, ptype));
    411:   PetscCall(MatProductSetFromOptions(D));
    412:   PetscFunctionReturn(PETSC_SUCCESS);
    413: }

    415: static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
    416: {
    417:   Mat A;

    419:   PetscFunctionBegin;
    420:   PetscCall(MatShellGetContext(N, &A));
    421:   PetscCall(MatGetDiagonal(A, v));
    422:   PetscFunctionReturn(PETSC_SUCCESS);
    423: }

    425: static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
    426: {
    427:   Mat a, b;

    429:   PetscFunctionBegin;
    430:   PetscCall(MatShellGetContext(A, &a));
    431:   PetscCall(MatShellGetContext(B, &b));
    432:   PetscCall(MatCopy(a, b, str));
    433:   PetscFunctionReturn(PETSC_SUCCESS);
    434: }

    436: static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
    437: {
    438:   Mat         A;
    439:   PetscScalar vscale = 1.0, vshift = 0.0;
    440:   PetscBool   flg;

    442:   PetscFunctionBegin;
    443:   PetscCall(MatShellGetContext(N, &A));
    444:   PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
    445:   if (flg || N->ops->getrow) { /* if this condition is false, MatConvert_Shell() will be called in MatConvert_Basic(), so the following checks are not needed */
    446:     PetscCall(MatShellGetScalingShifts(N, &vshift, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
    447:   }
    448:   if (flg) {
    449:     Mat B;

    451:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
    452:     if (reuse != MAT_INPLACE_MATRIX) {
    453:       PetscCall(MatConvert(B, newtype, reuse, newmat));
    454:       PetscCall(MatDestroy(&B));
    455:     } else {
    456:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
    457:       PetscCall(MatHeaderReplace(N, &B));
    458:     }
    459:   } else { /* use basic converter as fallback */
    460:     flg = (PetscBool)(N->ops->getrow != NULL);
    461:     PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
    462:   }
    463:   if (flg) {
    464:     PetscCall(MatScale(*newmat, vscale));
    465:     PetscCall(MatShift(*newmat, vshift));
    466:   }
    467:   PetscFunctionReturn(PETSC_SUCCESS);
    468: }

    470: static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
    471: {
    472:   PetscFunctionBegin;
    473:   PetscCall(MatShellGetContext(N, M));
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: /*@
    478:   MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`

    480:   Logically Collective

    482:   Input Parameter:
    483: . A - the `MATTRANSPOSEVIRTUAL` matrix

    485:   Output Parameter:
    486: . M - the matrix object stored inside `A`

    488:   Level: intermediate

    490: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
    491: @*/
    492: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
    493: {
    494:   PetscFunctionBegin;
    497:   PetscAssertPointer(M, 2);
    498:   PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: /*MC
    503:    MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix

    505:   Level: advanced

    507:   Developer Notes:
    508:   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

    510:   Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage

    512: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
    513:           `MATNORMALHERMITIAN`, `MATNORMAL`
    514: M*/

    516: /*@
    517:   MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'

    519:   Collective

    521:   Input Parameter:
    522: . A - the (possibly rectangular) matrix

    524:   Output Parameter:
    525: . N - the matrix that represents A'

    527:   Level: intermediate

    529:   Note:
    530:   The transpose A' is NOT actually formed! Rather the new matrix
    531:   object performs the matrix-vector product by using the `MatMultTranspose()` on
    532:   the original matrix

    534: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
    535:           `MATNORMALHERMITIAN`
    536: @*/
    537: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
    538: {
    539:   VecType vtype;

    541:   PetscFunctionBegin;
    542:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
    543:   PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
    544:   PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
    545:   PetscCall(MatSetType(*N, MATSHELL));
    546:   PetscCall(MatShellSetContext(*N, A));
    547:   PetscCall(PetscObjectReference((PetscObject)A));

    549:   PetscCall(MatSetBlockSizes(*N, A->cmap->bs, A->rmap->bs));
    550:   PetscCall(MatGetVecType(A, &vtype));
    551:   PetscCall(MatSetVecType(*N, vtype));
    552: #if defined(PETSC_HAVE_DEVICE)
    553:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
    554: #endif
    555:   PetscCall(MatSetUp(*N));

    557:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Transpose));
    558:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Transpose));
    559:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Transpose));
    560:   PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (PetscErrorCodeFn *)MatLUFactor_Transpose));
    561:   PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (PetscErrorCodeFn *)MatCholeskyFactor_Transpose));
    562:   PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (PetscErrorCodeFn *)MatGetFactor_Transpose));
    563:   PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (PetscErrorCodeFn *)MatGetInfo_Transpose));
    564:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_Transpose));
    565:   PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (PetscErrorCodeFn *)MatHasOperation_Transpose));
    566:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Transpose));
    567:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_Transpose));
    568:   PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (PetscErrorCodeFn *)MatConvert_Transpose));

    570:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
    571:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
    572:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_Transpose));
    573:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
    574:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
    575:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
    576:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
    577:   PetscFunctionReturn(PETSC_SUCCESS);
    578: }