Actual source code: htransm.c

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

      3: typedef struct {
      4:   PetscErrorCode (*numeric)(Mat);
      5:   PetscCtxDestroyFn *destroy;
      6:   Mat                B, D;
      7:   PetscScalar        scale;
      8:   PetscBool          conjugate;
      9:   void              *data;
     10: } MatProductCtx_HT;

     12: static PetscErrorCode MatProductCtxDestroy_HT(PetscCtxRt ptr)
     13: {
     14:   MatProductCtx_HT *data = *(MatProductCtx_HT **)ptr;
     15:   PetscContainer    container;

     17:   PetscFunctionBegin;
     18:   if (data->data) PetscCall((*data->destroy)(&data->data));
     19:   if (data->conjugate) PetscCall(MatDestroy(&data->B));
     20:   PetscCall(PetscObjectQuery((PetscObject)data->D, "MatProductCtx_HT", (PetscObject *)&container));
     21:   PetscCall(PetscContainerDestroy(&container));
     22:   PetscCall(PetscObjectCompose((PetscObject)data->D, "MatProductCtx_HT", NULL));
     23:   PetscCall(PetscFree(data));
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

     27: static PetscErrorCode MatProductNumeric_HT(Mat D)
     28: {
     29:   Mat_Product      *product;
     30:   Mat               B;
     31:   MatProductCtx_HT *data;
     32:   PetscContainer    container;

     34:   PetscFunctionBegin;
     35:   MatCheckProduct(D, 1);
     36:   PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
     37:   product = D->product;
     38:   PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_HT", (PetscObject *)&container));
     39:   PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_HT missing");
     40:   PetscCall(PetscContainerGetPointer(container, &data));
     41:   B    = product->B;
     42:   data = (MatProductCtx_HT *)product->data;
     43:   if (data->conjugate) {
     44:     PetscCall(MatCopy(product->B, data->B, SAME_NONZERO_PATTERN));
     45:     PetscCall(MatConjugate(data->B));
     46:     product->B = data->B;
     47:   }
     48:   product->data = data->data;
     49:   PetscCall((*data->numeric)(D));
     50:   if (data->conjugate) {
     51:     PetscCall(MatConjugate(D));
     52:     product->B = B;
     53:   }
     54:   PetscCall(MatScale(D, data->scale));
     55:   product->data = data;
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: static PetscErrorCode MatProductSymbolic_HT(Mat D)
     60: {
     61:   Mat_Product      *product;
     62:   Mat               B;
     63:   MatProductCtx_HT *data;
     64:   PetscContainer    container;

     66:   PetscFunctionBegin;
     67:   MatCheckProduct(D, 1);
     68:   product = D->product;
     69:   B       = product->B;
     70:   if (D->ops->productsymbolic == MatProductSymbolic_HT) {
     71:     PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
     72:     PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_HT", (PetscObject *)&container));
     73:     PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_HT missing");
     74:     PetscCall(PetscContainerGetPointer(container, &data));
     75:     PetscCall(MatProductSetFromOptions(D));
     76:     if (data->conjugate) {
     77:       PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &data->B));
     78:       product->B = data->B;
     79:     }
     80:     PetscCall(MatProductSymbolic(D));
     81:     data->numeric          = D->ops->productnumeric;
     82:     data->destroy          = product->destroy;
     83:     data->data             = product->data;
     84:     D->ops->productnumeric = MatProductNumeric_HT;
     85:     product->destroy       = MatProductCtxDestroy_HT;
     86:     if (data->conjugate) product->B = B;
     87:     product->data = data;
     88:   }
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: static PetscErrorCode MatProductSetFromOptions_HT(Mat D)
     93: {
     94:   Mat               A, B, C, Ain, Bin, Cin;
     95:   PetscScalar       scale = 1.0, vscale;
     96:   PetscBool         Aistrans, Bistrans, Cistrans, conjugate = PETSC_FALSE;
     97:   PetscInt          Atrans, Btrans, Ctrans;
     98:   PetscContainer    container = NULL;
     99:   MatProductCtx_HT *data;
    100:   MatProductType    ptype;

    102:   PetscFunctionBegin;
    103:   MatCheckProduct(D, 1);
    104:   A = D->product->A;
    105:   B = D->product->B;
    106:   C = D->product->C;
    107:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
    108:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
    109:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
    110:   PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
    111:   Atrans = 0;
    112:   Ain    = A;
    113:   while (Aistrans) {
    114:     Atrans++;
    115:     PetscCall(MatShellGetScalingShifts(Ain, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &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));
    116:     conjugate = (PetscBool)!conjugate;
    117:     scale *= vscale;
    118:     PetscCall(MatHermitianTransposeGetMat(Ain, &Ain));
    119:     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
    120:   }
    121:   Btrans = 0;
    122:   Bin    = B;
    123:   while (Bistrans) {
    124:     Btrans++;
    125:     PetscCall(MatShellGetScalingShifts(Bin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &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));
    126:     scale *= vscale;
    127:     PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
    128:     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
    129:   }
    130:   Ctrans = 0;
    131:   Cin    = C;
    132:   while (Cistrans) {
    133:     Ctrans++;
    134:     PetscCall(MatShellGetScalingShifts(Cin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &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));
    135:     scale *= vscale;
    136:     PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
    137:     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
    138:   }
    139:   Atrans = Atrans % 2;
    140:   Btrans = Btrans % 2;
    141:   Ctrans = Ctrans % 2;
    142:   ptype  = D->product->type; /* same product type by default */
    143:   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
    144:   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
    145:   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;

    147:   if (Atrans || Btrans || Ctrans) {
    148:     PetscCheck(!PetscDefined(USE_COMPLEX) || (!Btrans && !Ctrans), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for complex Hermitian transpose matrices");
    149:     if ((PetscDefined(USE_COMPLEX) && Atrans) || scale != 1.0) {
    150:       PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_HT", (PetscObject *)&container));
    151:       if (!container) {
    152:         PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)D), &container));
    153:         PetscCall(PetscNew(&data));
    154:         PetscCall(PetscContainerSetPointer(container, data));
    155:         PetscCall(PetscObjectCompose((PetscObject)D, "MatProductCtx_HT", (PetscObject)container));
    156:       } else PetscCall(PetscContainerGetPointer(container, &data));
    157:       data->scale     = scale;
    158:       data->conjugate = (PetscBool)Atrans;
    159:       data->D         = D;
    160:     }
    161:     ptype = MATPRODUCT_UNSPECIFIED;
    162:     switch (D->product->type) {
    163:     case MATPRODUCT_AB:
    164:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
    165:         /* TODO custom implementation ? */
    166:       } else if (Atrans) { /* At * B */
    167:         ptype = MATPRODUCT_AtB;
    168:       } else { /* A * Bt */
    169:         ptype = MATPRODUCT_ABt;
    170:       }
    171:       break;
    172:     case MATPRODUCT_AtB:
    173:       if (Atrans && Btrans) { /* A * Bt */
    174:         ptype = MATPRODUCT_ABt;
    175:       } else if (Atrans) { /* A * B */
    176:         ptype = MATPRODUCT_AB;
    177:       } else { /* At * Bt we do not have support for this */
    178:         /* TODO custom implementation ? */
    179:       }
    180:       break;
    181:     case MATPRODUCT_ABt:
    182:       if (Atrans && Btrans) { /* At * B */
    183:         ptype = MATPRODUCT_AtB;
    184:       } else if (Atrans) { /* At * Bt we do not have support for this */
    185:         /* TODO custom implementation ? */
    186:       } else { /* A * B */
    187:         ptype = MATPRODUCT_AB;
    188:       }
    189:       break;
    190:     case MATPRODUCT_PtAP:
    191:       if (Atrans) { /* PtAtP */
    192:         /* TODO custom implementation ? */
    193:       } else { /* RARt */
    194:         ptype = MATPRODUCT_RARt;
    195:       }
    196:       break;
    197:     case MATPRODUCT_RARt:
    198:       if (Atrans) { /* RAtRt */
    199:         /* TODO custom implementation ? */
    200:       } else { /* PtAP */
    201:         ptype = MATPRODUCT_PtAP;
    202:       }
    203:       break;
    204:     case MATPRODUCT_ABC:
    205:       /* TODO custom implementation ? */
    206:       break;
    207:     default:
    208:       SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
    209:     }
    210:   }
    211:   PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
    212:   PetscCall(MatProductSetType(D, ptype));
    213:   if (container) D->ops->productsymbolic = MatProductSymbolic_HT;
    214:   else PetscCall(MatProductSetFromOptions(D));
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: static PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
    219: {
    220:   Mat A;

    222:   PetscFunctionBegin;
    223:   PetscCall(MatShellGetContext(N, &A));
    224:   PetscCall(MatMultHermitianTranspose(A, x, y));
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }

    228: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
    229: {
    230:   Mat A;

    232:   PetscFunctionBegin;
    233:   PetscCall(MatShellGetContext(N, &A));
    234:   PetscCall(MatMult(A, x, y));
    235:   PetscFunctionReturn(PETSC_SUCCESS);
    236: }

    238: static PetscErrorCode MatSolve_HT_LU(Mat N, Vec b, Vec x)
    239: {
    240:   Mat A;
    241:   Vec w;

    243:   PetscFunctionBegin;
    244:   PetscCall(MatShellGetContext(N, &A));
    245:   PetscCall(VecDuplicate(b, &w));
    246:   PetscCall(VecCopy(b, w));
    247:   PetscCall(VecConjugate(w));
    248:   PetscCall(MatSolveTranspose(A, w, x));
    249:   PetscCall(VecConjugate(x));
    250:   PetscCall(VecDestroy(&w));
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: static PetscErrorCode MatSolveAdd_HT_LU(Mat N, Vec b, Vec y, Vec x)
    255: {
    256:   Mat A;
    257:   Vec v, w;

    259:   PetscFunctionBegin;
    260:   PetscCall(MatShellGetContext(N, &A));
    261:   PetscCall(VecDuplicate(b, &v));
    262:   PetscCall(VecDuplicate(b, &w));
    263:   PetscCall(VecCopy(y, v));
    264:   PetscCall(VecCopy(b, w));
    265:   PetscCall(VecConjugate(v));
    266:   PetscCall(VecConjugate(w));
    267:   PetscCall(MatSolveTransposeAdd(A, w, v, x));
    268:   PetscCall(VecConjugate(x));
    269:   PetscCall(VecDestroy(&v));
    270:   PetscCall(VecDestroy(&w));
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: static PetscErrorCode MatMatSolve_HT_LU(Mat N, Mat B, Mat X)
    275: {
    276:   Mat A, W;

    278:   PetscFunctionBegin;
    279:   PetscCall(MatShellGetContext(N, &A));
    280:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &W));
    281:   PetscCall(MatConjugate(W));
    282:   PetscCall(MatMatSolveTranspose(A, W, X));
    283:   PetscCall(MatConjugate(X));
    284:   PetscCall(MatDestroy(&W));
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: static PetscErrorCode MatLUFactor_HT(Mat N, IS row, IS col, const MatFactorInfo *minfo)
    289: {
    290:   Mat A;

    292:   PetscFunctionBegin;
    293:   PetscCall(MatShellGetContext(N, &A));
    294:   PetscCall(MatLUFactor(A, col, row, minfo));
    295:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_HT_LU));
    296:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_HT_LU));
    297:   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_HT_LU));
    298:   PetscFunctionReturn(PETSC_SUCCESS);
    299: }

    301: static PetscErrorCode MatSolve_HT_Cholesky(Mat N, Vec b, Vec x)
    302: {
    303:   Mat A;

    305:   PetscFunctionBegin;
    306:   PetscCall(MatShellGetContext(N, &A));
    307:   PetscCall(MatSolve(A, b, x));
    308:   PetscFunctionReturn(PETSC_SUCCESS);
    309: }

    311: static PetscErrorCode MatSolveAdd_HT_Cholesky(Mat N, Vec b, Vec y, Vec x)
    312: {
    313:   Mat A;
    314:   Vec v, w;

    316:   PetscFunctionBegin;
    317:   PetscCall(MatShellGetContext(N, &A));
    318:   PetscCall(VecDuplicate(b, &v));
    319:   PetscCall(VecDuplicate(b, &w));
    320:   PetscCall(VecCopy(y, v));
    321:   PetscCall(VecCopy(b, w));
    322:   PetscCall(VecConjugate(v));
    323:   PetscCall(VecConjugate(w));
    324:   PetscCall(MatSolveTransposeAdd(A, w, v, x));
    325:   PetscCall(VecConjugate(x));
    326:   PetscCall(VecDestroy(&v));
    327:   PetscCall(VecDestroy(&w));
    328:   PetscFunctionReturn(PETSC_SUCCESS);
    329: }

    331: static PetscErrorCode MatMatSolve_HT_Cholesky(Mat N, Mat B, Mat X)
    332: {
    333:   Mat A, W;

    335:   PetscFunctionBegin;
    336:   PetscCall(MatShellGetContext(N, &A));
    337:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &W));
    338:   PetscCall(MatConjugate(W));
    339:   PetscCall(MatMatSolveTranspose(A, W, X));
    340:   PetscCall(MatConjugate(X));
    341:   PetscCall(MatDestroy(&W));
    342:   PetscFunctionReturn(PETSC_SUCCESS);
    343: }

    345: static PetscErrorCode MatCholeskyFactor_HT(Mat N, IS perm, const MatFactorInfo *minfo)
    346: {
    347:   Mat A;

    349:   PetscFunctionBegin;
    350:   PetscCall(MatShellGetContext(N, &A));
    351:   PetscCheck(!PetscDefined(USE_COMPLEX) || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cholesky supported only if original matrix is Hermitian");
    352:   PetscCall(MatCholeskyFactor(A, perm, minfo));
    353:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_HT_Cholesky));
    354:   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_HT_Cholesky));
    355:   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_HT_Cholesky));
    356:   PetscFunctionReturn(PETSC_SUCCESS);
    357: }

    359: static PetscErrorCode MatLUFactorNumeric_HT(Mat F, Mat N, const MatFactorInfo *info)
    360: {
    361:   Mat A, FA;

    363:   PetscFunctionBegin;
    364:   PetscCall(MatShellGetContext(N, &A));
    365:   PetscCall(MatShellGetContext(F, &FA));
    366:   PetscCall(MatLUFactorNumeric(FA, A, info));
    367:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_HT_LU));
    368:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_HT_LU));
    369:   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_HT_LU));
    370:   PetscFunctionReturn(PETSC_SUCCESS);
    371: }

    373: static PetscErrorCode MatLUFactorSymbolic_HT(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info)
    374: {
    375:   Mat A, FA;

    377:   PetscFunctionBegin;
    378:   PetscCall(MatShellGetContext(N, &A));
    379:   PetscCall(MatShellGetContext(F, &FA));
    380:   PetscCall(MatLUFactorSymbolic(FA, A, row, col, info));
    381:   PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (PetscErrorCodeFn *)MatLUFactorNumeric_HT));
    382:   PetscFunctionReturn(PETSC_SUCCESS);
    383: }

    385: static PetscErrorCode MatCholeskyFactorNumeric_HT(Mat F, Mat N, const MatFactorInfo *info)
    386: {
    387:   Mat A, FA;

    389:   PetscFunctionBegin;
    390:   PetscCall(MatShellGetContext(N, &A));
    391:   PetscCall(MatShellGetContext(F, &FA));
    392:   PetscCall(MatCholeskyFactorNumeric(FA, A, info));
    393:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_HT_Cholesky));
    394:   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_HT_Cholesky));
    395:   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_HT_Cholesky));
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: static PetscErrorCode MatCholeskyFactorSymbolic_HT(Mat F, Mat N, IS perm, const MatFactorInfo *info)
    400: {
    401:   Mat A, FA;

    403:   PetscFunctionBegin;
    404:   PetscCall(MatShellGetContext(N, &A));
    405:   PetscCall(MatShellGetContext(F, &FA));
    406:   PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info));
    407:   PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (PetscErrorCodeFn *)MatCholeskyFactorNumeric_HT));
    408:   PetscFunctionReturn(PETSC_SUCCESS);
    409: }

    411: static PetscErrorCode MatGetFactor_HT(Mat N, MatSolverType type, MatFactorType ftype, Mat *F)
    412: {
    413:   Mat A, FA;

    415:   PetscFunctionBegin;
    416:   PetscCall(MatShellGetContext(N, &A));
    417:   PetscCall(MatGetFactor(A, type, ftype, &FA));
    418:   PetscCall(MatCreateTranspose(FA, F));
    419:   if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatLUFactorSymbolic_HT));
    420:   else if (ftype == MAT_FACTOR_CHOLESKY) {
    421:     PetscCheck(!PetscDefined(USE_COMPLEX) || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cholesky supported only if original matrix is Hermitian");
    422:     PetscCall(MatPropagateSymmetryOptions(A, FA));
    423:     PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatCholeskyFactorSymbolic_HT));
    424:   } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]);
    425:   (*F)->factortype = ftype;
    426:   PetscCall(MatDestroy(&FA));
    427:   PetscFunctionReturn(PETSC_SUCCESS);
    428: }

    430: static PetscErrorCode MatDestroy_HT(Mat N)
    431: {
    432:   Mat A;

    434:   PetscFunctionBegin;
    435:   PetscCall(MatShellGetContext(N, &A));
    436:   PetscCall(MatDestroy(&A));
    437:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
    438: #if !defined(PETSC_USE_COMPLEX)
    439:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
    440: #endif
    441:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
    442:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
    443:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL));
    444:   PetscFunctionReturn(PETSC_SUCCESS);
    445: }

    447: static PetscErrorCode MatGetInfo_HT(Mat N, MatInfoType flag, MatInfo *info)
    448: {
    449:   Mat A;

    451:   PetscFunctionBegin;
    452:   PetscCall(MatShellGetContext(N, &A));
    453:   PetscCall(MatGetInfo(A, flag, info));
    454:   PetscFunctionReturn(PETSC_SUCCESS);
    455: }

    457: static PetscErrorCode MatFactorGetSolverType_HT(Mat N, MatSolverType *type)
    458: {
    459:   Mat A;

    461:   PetscFunctionBegin;
    462:   PetscCall(MatShellGetContext(N, &A));
    463:   PetscCall(MatFactorGetSolverType(A, type));
    464:   PetscFunctionReturn(PETSC_SUCCESS);
    465: }

    467: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
    468: {
    469:   Mat A, C;

    471:   PetscFunctionBegin;
    472:   PetscCall(MatShellGetContext(N, &A));
    473:   PetscCall(MatDuplicate(A, op, &C));
    474:   PetscCall(MatCreateHermitianTranspose(C, m));
    475:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
    476:   PetscCall(MatDestroy(&C));
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }

    480: static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has)
    481: {
    482:   Mat A;

    484:   PetscFunctionBegin;
    485:   PetscCall(MatShellGetContext(mat, &A));
    486:   *has = PETSC_FALSE;
    487:   if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
    488:     PetscCall(MatHasOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, has));
    489:     if (!*has) PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
    490:   } else if (op == MATOP_MULT_HERMITIAN_TRANSPOSE || op == MATOP_MULT_HERMITIAN_TRANS_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
    491:     PetscCall(MatHasOperation(A, MATOP_MULT, has));
    492:   } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
    493:   PetscFunctionReturn(PETSC_SUCCESS);
    494: }

    496: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
    497: {
    498:   PetscFunctionBegin;
    499:   PetscCall(MatShellGetContext(N, M));
    500:   PetscFunctionReturn(PETSC_SUCCESS);
    501: }

    503: /*@
    504:   MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`

    506:   Logically Collective

    508:   Input Parameter:
    509: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix

    511:   Output Parameter:
    512: . M - the matrix object stored inside A

    514:   Level: intermediate

    516: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
    517: @*/
    518: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
    519: {
    520:   PetscFunctionBegin;
    523:   PetscAssertPointer(M, 2);
    524:   PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
    525:   PetscFunctionReturn(PETSC_SUCCESS);
    526: }

    528: static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v)
    529: {
    530:   Mat A;

    532:   PetscFunctionBegin;
    533:   PetscCall(MatShellGetContext(N, &A));
    534:   PetscCall(MatGetDiagonal(A, v));
    535:   PetscCall(VecConjugate(v));
    536:   PetscFunctionReturn(PETSC_SUCCESS);
    537: }

    539: static PetscErrorCode MatCopy_HT(Mat A, Mat B, MatStructure str)
    540: {
    541:   Mat a, b;

    543:   PetscFunctionBegin;
    544:   PetscCall(MatShellGetContext(A, &a));
    545:   PetscCall(MatShellGetContext(B, &b));
    546:   PetscCall(MatCopy(a, b, str));
    547:   PetscFunctionReturn(PETSC_SUCCESS);
    548: }

    550: static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
    551: {
    552:   Mat         A;
    553:   PetscScalar vscale = 1.0, vshift = 0.0;
    554:   PetscBool   flg;

    556:   PetscFunctionBegin;
    557:   PetscCall(MatShellGetContext(N, &A));
    558:   PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg));
    559:   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 */
    560:     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));
    561:   }
    562:   if (flg) {
    563:     Mat B;

    565:     PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B));
    566:     if (reuse != MAT_INPLACE_MATRIX) {
    567:       PetscCall(MatConvert(B, newtype, reuse, newmat));
    568:       PetscCall(MatDestroy(&B));
    569:     } else {
    570:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
    571:       PetscCall(MatHeaderReplace(N, &B));
    572:     }
    573:   } else { /* use basic converter as fallback */
    574:     flg = (PetscBool)(N->ops->getrow != NULL);
    575:     PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
    576:   }
    577:   if (flg) {
    578:     PetscCall(MatScale(*newmat, vscale));
    579:     PetscCall(MatShift(*newmat, vshift));
    580:   }
    581:   PetscFunctionReturn(PETSC_SUCCESS);
    582: }

    584: /*MC
    585:    MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix

    587:   Level: advanced

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

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

    594: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`
    595: M*/

    597: /*@
    598:   MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'*

    600:   Collective

    602:   Input Parameter:
    603: . A - the (possibly rectangular) matrix

    605:   Output Parameter:
    606: . N - the matrix that represents A'*

    608:   Level: intermediate

    610:   Note:
    611:   The Hermitian transpose A' is NOT actually formed! Rather the new matrix
    612:   object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
    613:   the original matrix

    615: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
    616:           `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
    617: @*/
    618: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
    619: {
    620:   VecType vtype;

    622:   PetscFunctionBegin;
    623:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
    624:   PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
    625:   PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
    626:   PetscCall(MatSetType(*N, MATSHELL));
    627:   PetscCall(MatShellSetContext(*N, A));
    628:   PetscCall(PetscObjectReference((PetscObject)A));

    630:   PetscCall(MatSetBlockSizes(*N, A->cmap->bs, A->rmap->bs));
    631:   PetscCall(MatGetVecType(A, &vtype));
    632:   PetscCall(MatSetVecType(*N, vtype));
    633: #if defined(PETSC_HAVE_DEVICE)
    634:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
    635: #endif
    636:   PetscCall(MatSetUp(*N));

    638:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_HT));
    639:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_HT));
    640:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (PetscErrorCodeFn *)MatMultHermitianTranspose_HT));
    641: #if !defined(PETSC_USE_COMPLEX)
    642:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultHermitianTranspose_HT));
    643: #endif
    644:   PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (PetscErrorCodeFn *)MatLUFactor_HT));
    645:   PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (PetscErrorCodeFn *)MatCholeskyFactor_HT));
    646:   PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (PetscErrorCodeFn *)MatGetFactor_HT));
    647:   PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (PetscErrorCodeFn *)MatGetInfo_HT));
    648:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_HT));
    649:   PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (PetscErrorCodeFn *)MatHasOperation_HT));
    650:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_HT));
    651:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_HT));
    652:   PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (PetscErrorCodeFn *)MatConvert_HT));

    654:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
    655: #if !defined(PETSC_USE_COMPLEX)
    656:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
    657: #endif
    658:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT));
    659:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_HT));
    660:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
    661:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
    662:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
    663:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
    664:   PetscFunctionReturn(PETSC_SUCCESS);
    665: }