Actual source code: normmh.c

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

      3: typedef struct {
      4:   Mat A;
      5:   Mat D; /* local submatrix for diagonal part */
      6:   Vec w;
      7: } Mat_NormalHermitian;

      9: static PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
     10: {
     11:   Mat_NormalHermitian *a;
     12:   Mat                  B, *suba;
     13:   IS                  *row;
     14:   PetscScalar          shift, scale;
     15:   PetscInt             M;

     17:   PetscFunctionBegin;
     18:   PetscCheck(irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
     19:   PetscCall(MatShellGetScalingShifts(mat, &shift, &scale, (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));
     20:   PetscCall(MatShellGetContext(mat, &a));
     21:   B = a->A;
     22:   if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
     23:   PetscCall(MatGetSize(B, &M, NULL));
     24:   PetscCall(PetscMalloc1(n, &row));
     25:   PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
     26:   PetscCall(ISSetIdentity(row[0]));
     27:   for (M = 1; M < n; ++M) row[M] = row[0];
     28:   PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
     29:   for (M = 0; M < n; ++M) {
     30:     PetscCall(MatCreateNormalHermitian(suba[M], *submat + M));
     31:     PetscCall(MatShift((*submat)[M], shift));
     32:     PetscCall(MatScale((*submat)[M], scale));
     33:   }
     34:   PetscCall(ISDestroy(&row[0]));
     35:   PetscCall(PetscFree(row));
     36:   PetscCall(MatDestroySubMatrices(n, &suba));
     37:   PetscFunctionReturn(PETSC_SUCCESS);
     38: }

     40: static PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B)
     41: {
     42:   Mat_NormalHermitian *a;
     43:   Mat                  C, Aa;
     44:   IS                   row;
     45:   PetscScalar          shift, scale;

     47:   PetscFunctionBegin;
     48:   PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
     49:   PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (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));
     50:   PetscCall(MatShellGetContext(A, &a));
     51:   Aa = a->A;
     52:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
     53:   PetscCall(ISSetIdentity(row));
     54:   PetscCall(MatPermute(Aa, row, colp, &C));
     55:   PetscCall(ISDestroy(&row));
     56:   PetscCall(MatCreateNormalHermitian(C, B));
     57:   PetscCall(MatDestroy(&C));
     58:   PetscCall(MatShift(*B, shift));
     59:   PetscCall(MatScale(*B, scale));
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: static PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
     64: {
     65:   Mat_NormalHermitian *a;
     66:   Mat                  C;

     68:   PetscFunctionBegin;
     69:   PetscCall(MatShellGetContext(A, &a));
     70:   PetscCall(MatDuplicate(a->A, op, &C));
     71:   PetscCall(MatCreateNormalHermitian(C, B));
     72:   PetscCall(MatDestroy(&C));
     73:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: static PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
     78: {
     79:   Mat_NormalHermitian *a, *b;

     81:   PetscFunctionBegin;
     82:   PetscCall(MatShellGetContext(A, &a));
     83:   PetscCall(MatShellGetContext(B, &b));
     84:   PetscCall(MatCopy(a->A, b->A, str));
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: static PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
     89: {
     90:   Mat_NormalHermitian *Na;

     92:   PetscFunctionBegin;
     93:   PetscCall(MatShellGetContext(N, &Na));
     94:   PetscCall(MatMult(Na->A, x, Na->w));
     95:   PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode MatDestroy_NormalHermitian(Mat N)
    100: {
    101:   Mat_NormalHermitian *Na;

    103:   PetscFunctionBegin;
    104:   PetscCall(MatShellGetContext(N, &Na));
    105:   PetscCall(MatDestroy(&Na->A));
    106:   PetscCall(MatDestroy(&Na->D));
    107:   PetscCall(VecDestroy(&Na->w));
    108:   PetscCall(PetscFree(Na));
    109:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalHermitianGetMat_C", NULL));
    110: #if !defined(PETSC_USE_COMPLEX)
    111:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL));
    112: #endif
    113:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL));
    114:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL));
    115: #if defined(PETSC_HAVE_HYPRE)
    116:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_hypre_C", NULL));
    117: #endif
    118:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: /*
    123:       Slow, nonscalable version
    124: */
    125: static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
    126: {
    127:   Mat_NormalHermitian *Na;
    128:   Mat                  A;
    129:   PetscInt             i, j, rstart, rend, nnz;
    130:   const PetscInt      *cols;
    131:   PetscScalar         *work, *values;
    132:   const PetscScalar   *mvalues;

    134:   PetscFunctionBegin;
    135:   PetscCall(MatShellGetContext(N, &Na));
    136:   A = Na->A;
    137:   PetscCall(PetscMalloc1(A->cmap->N, &work));
    138:   PetscCall(PetscArrayzero(work, A->cmap->N));
    139:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
    140:   for (i = rstart; i < rend; i++) {
    141:     PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
    142:     for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
    143:     PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
    144:   }
    145:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, work, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
    146:   rstart = N->cmap->rstart;
    147:   rend   = N->cmap->rend;
    148:   PetscCall(VecGetArray(v, &values));
    149:   PetscCall(PetscArraycpy(values, work + rstart, rend - rstart));
    150:   PetscCall(VecRestoreArray(v, &values));
    151:   PetscCall(PetscFree(work));
    152:   PetscFunctionReturn(PETSC_SUCCESS);
    153: }

    155: static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
    156: {
    157:   Mat_NormalHermitian *Na;
    158:   Mat                  M, A;

    160:   PetscFunctionBegin;
    161:   PetscCall(MatShellGetContext(N, &Na));
    162:   A = Na->A;
    163:   PetscCall(MatGetDiagonalBlock(A, &M));
    164:   PetscCall(MatCreateNormalHermitian(M, &Na->D));
    165:   *D = Na->D;
    166:   PetscFunctionReturn(PETSC_SUCCESS);
    167: }

    169: static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M)
    170: {
    171:   Mat_NormalHermitian *Aa;

    173:   PetscFunctionBegin;
    174:   PetscCall(MatShellGetContext(A, &Aa));
    175:   *M = Aa->A;
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: /*@
    180:   MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`

    182:   Logically Collective

    184:   Input Parameter:
    185: . A - the `MATNORMALHERMITIAN` matrix

    187:   Output Parameter:
    188: . M - the matrix object stored inside `A`

    190:   Level: intermediate

    192: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
    193: @*/
    194: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
    195: {
    196:   PetscFunctionBegin;
    199:   PetscAssertPointer(M, 2);
    200:   PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M));
    201:   PetscFunctionReturn(PETSC_SUCCESS);
    202: }

    204: static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
    205: {
    206:   Mat_NormalHermitian *Aa;
    207:   Mat                  B, conjugate;
    208:   Vec                  left, right, dshift;
    209:   PetscScalar          scale, shift;
    210:   PetscInt             m, n, M, N;

    212:   PetscFunctionBegin;
    213:   PetscCall(MatShellGetContext(A, &Aa));
    214:   PetscCall(MatShellGetScalingShifts(A, &shift, &scale, &dshift, &left, &right, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
    215:   PetscCall(MatGetSize(A, &M, &N));
    216:   PetscCall(MatGetLocalSize(A, &m, &n));
    217:   if (reuse == MAT_REUSE_MATRIX) {
    218:     B = *newmat;
    219:     PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
    220:   } else {
    221:     PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
    222:     PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
    223:     PetscCall(MatProductSetFromOptions(B));
    224:     PetscCall(MatProductSymbolic(B));
    225:     PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
    226:   }
    227:   if (PetscDefined(USE_COMPLEX)) {
    228:     PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
    229:     PetscCall(MatConjugate(conjugate));
    230:     PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
    231:   }
    232:   PetscCall(MatProductNumeric(B));
    233:   if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
    234:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
    235:   else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
    236:   PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
    237:   PetscCall(MatDiagonalScale(*newmat, left, right));
    238:   PetscCall(MatScale(*newmat, scale));
    239:   PetscCall(MatShift(*newmat, shift));
    240:   if (dshift) PetscCall(MatDiagonalSet(*newmat, dshift, ADD_VALUES));
    241:   PetscFunctionReturn(PETSC_SUCCESS);
    242: }

    244: #if defined(PETSC_HAVE_HYPRE)
    245: static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
    246: {
    247:   PetscFunctionBegin;
    248:   if (reuse == MAT_INITIAL_MATRIX) {
    249:     PetscCall(MatConvert(A, MATAIJ, reuse, B));
    250:     PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
    251:   } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
    252:   PetscFunctionReturn(PETSC_SUCCESS);
    253: }
    254: #endif

    256: /*MC
    257:   MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A

    259:   Level: intermediate

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

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

    266: .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
    267: M*/

    269: /*@
    270:   MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like $A^* A$.

    272:   Collective

    274:   Input Parameter:
    275: . A - the (possibly rectangular complex) matrix

    277:   Output Parameter:
    278: . N - the matrix that represents $ A^* A$

    280:   Level: intermediate

    282:   Note:
    283:   The product $ A^* A$ is NOT actually formed! Rather the new matrix
    284:   object performs the matrix-vector product, `MatMult()`, by first multiplying by
    285:   $A$ and then $A^*$

    287:   If `MatGetFactor()` is called on this matrix with `MAT_FACTOR_QR` then the inner matrix `A` is used for the factorization

    289: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
    290: @*/
    291: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
    292: {
    293:   Mat_NormalHermitian *Na;
    294:   VecType              vtype;

    296:   PetscFunctionBegin;
    297:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
    298:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
    299:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
    300:   PetscCall(MatSetType(*N, MATSHELL));
    301:   PetscCall(PetscNew(&Na));
    302:   PetscCall(MatShellSetContext(*N, Na));
    303:   PetscCall(PetscObjectReference((PetscObject)A));
    304:   Na->A = A;
    305:   PetscCall(MatCreateVecs(A, NULL, &Na->w));

    307:   PetscCall(MatSetBlockSize(*N, A->cmap->bs));
    308:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_NormalHermitian));
    309:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_NormalHermitian));
    310:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (PetscErrorCodeFn *)MatMult_NormalHermitian));
    311: #if !defined(PETSC_USE_COMPLEX)
    312:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMult_NormalHermitian));
    313: #endif
    314:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_NormalHermitian));
    315:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_NormalHermitian));
    316:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (PetscErrorCodeFn *)MatGetDiagonalBlock_NormalHermitian));
    317:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_NormalHermitian));
    318:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
    319:   (*N)->ops->permute           = MatPermute_NormalHermitian;

    321:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalHermitianGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
    322: #if !defined(PETSC_USE_COMPLEX)
    323:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
    324: #endif
    325:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ));
    326:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ));
    327: #if defined(PETSC_HAVE_HYPRE)
    328:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE));
    329: #endif
    330:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
    331:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
    332:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
    333:   PetscCall(MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE));
    334:   PetscCall(MatGetVecType(A, &vtype));
    335:   PetscCall(MatSetVecType(*N, vtype));
    336: #if defined(PETSC_HAVE_DEVICE)
    337:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
    338: #endif
    339:   PetscCall(MatSetUp(*N));
    340:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN));
    341:   PetscFunctionReturn(PETSC_SUCCESS);
    342: }