Actual source code: axpy.c

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

      3: static PetscErrorCode MatTransposeAXPY_Private(Mat Y, PetscScalar a, Mat X, MatStructure str, Mat T)
      4: {
      5:   Mat         A, F;
      6:   PetscScalar vshift, vscale;
      7:   PetscErrorCode (*f)(Mat, Mat *);

      9:   PetscFunctionBegin;
     10:   if (T == X) PetscCall(MatShellGetScalingShifts(T, &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));
     11:   else {
     12:     vshift = 0.0;
     13:     vscale = 1.0;
     14:   }
     15:   PetscCall(PetscObjectQueryFunction((PetscObject)T, "MatTransposeGetMat_C", &f));
     16:   if (f) {
     17:     PetscCall(MatTransposeGetMat(T, &A));
     18:     if (T == X) {
     19:       PetscCall(PetscInfo(NULL, "Explicitly transposing X of type MATTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
     20:       PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &F));
     21:       A = Y;
     22:     } else {
     23:       PetscCall(PetscInfo(NULL, "Transposing X because Y of type MATTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
     24:       PetscCall(MatTranspose(X, MAT_INITIAL_MATRIX, &F));
     25:     }
     26:   } else {
     27:     PetscCall(MatHermitianTransposeGetMat(T, &A));
     28:     if (T == X) {
     29:       PetscCall(PetscInfo(NULL, "Explicitly Hermitian transposing X of type MATHERMITIANTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
     30:       PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &F));
     31:       A = Y;
     32:     } else {
     33:       PetscCall(PetscInfo(NULL, "Hermitian transposing X because Y of type MATHERMITIANTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
     34:       PetscCall(MatHermitianTranspose(X, MAT_INITIAL_MATRIX, &F));
     35:     }
     36:   }
     37:   PetscCall(MatAXPY(A, a * vscale, F, str));
     38:   PetscCall(MatShift(A, a * vshift));
     39:   PetscCall(MatDestroy(&F));
     40:   PetscFunctionReturn(PETSC_SUCCESS);
     41: }

     43: static PetscErrorCode MatAXPY_BasicWithTypeCompare(Mat Y, PetscScalar a, Mat X, MatStructure str)
     44: {
     45:   PetscBool flg;

     47:   PetscFunctionBegin;
     48:   PetscCall(MatIsShell(Y, &flg));
     49:   if (flg) { /* MatShell has special support for AXPY */
     50:     PetscErrorCode (*f)(Mat, PetscScalar, Mat, MatStructure);

     52:     PetscCall(MatGetOperation(Y, MATOP_AXPY, (void (**)(void))&f));
     53:     if (f) {
     54:       PetscCall((*f)(Y, a, X, str));
     55:       PetscFunctionReturn(PETSC_SUCCESS);
     56:     }
     57:   } else {
     58:     /* no need to preallocate if Y is dense */
     59:     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &flg, MATSEQDENSE, MATMPIDENSE, ""));
     60:     if (flg) {
     61:       PetscCall(PetscObjectTypeCompare((PetscObject)X, MATNEST, &flg));
     62:       if (flg) {
     63:         PetscCall(MatAXPY_Dense_Nest(Y, a, X));
     64:         PetscFunctionReturn(PETSC_SUCCESS);
     65:       }
     66:     }
     67:     PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSCALAPACK, MATELEMENTAL, ""));
     68:     if (flg) { /* Avoid MatAXPY_Basic() due to missing MatGetRow() */
     69:       Mat C;

     71:       PetscCall(MatConvert(X, ((PetscObject)Y)->type_name, MAT_INITIAL_MATRIX, &C));
     72:       PetscCall(MatAXPY(Y, a, C, str));
     73:       PetscCall(MatDestroy(&C));
     74:       PetscFunctionReturn(PETSC_SUCCESS);
     75:     }
     76:   }
     77:   PetscCall(MatAXPY_Basic(Y, a, X, str));
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: /*@
     82:   MatAXPY - Computes Y = a*X + Y.

     84:   Logically Collective

     86:   Input Parameters:
     87: + a   - the scalar multiplier
     88: . X   - the first matrix
     89: . Y   - the second matrix
     90: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)

     92:   Level: intermediate

     94: .seealso: [](ch_matrices), `Mat`, `MatAYPX()`
     95:  @*/
     96: PetscErrorCode MatAXPY(Mat Y, PetscScalar a, Mat X, MatStructure str)
     97: {
     98:   PetscInt  M1, M2, N1, N2;
     99:   PetscInt  m1, m2, n1, n2;
    100:   PetscBool sametype, transpose;

    102:   PetscFunctionBegin;
    106:   PetscCheckSameComm(Y, 1, X, 3);
    107:   PetscCall(MatGetSize(X, &M1, &N1));
    108:   PetscCall(MatGetSize(Y, &M2, &N2));
    109:   PetscCall(MatGetLocalSize(X, &m1, &n1));
    110:   PetscCall(MatGetLocalSize(Y, &m2, &n2));
    111:   PetscCheck(M1 == M2 && N1 == N2, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_SIZ, "Non conforming matrix add: global sizes X %" PetscInt_FMT " x %" PetscInt_FMT ", Y %" PetscInt_FMT " x %" PetscInt_FMT, M1, N1, M2, N2);
    112:   PetscCheck(m1 == m2 && n1 == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Non conforming matrix add: local sizes X %" PetscInt_FMT " x %" PetscInt_FMT ", Y %" PetscInt_FMT " x %" PetscInt_FMT, m1, n1, m2, n2);
    113:   PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix (Y)");
    114:   PetscCheck(X->assembled, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix (X)");
    115:   if (a == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
    116:   if (Y == X) {
    117:     PetscCall(MatScale(Y, 1.0 + a));
    118:     PetscFunctionReturn(PETSC_SUCCESS);
    119:   }
    120:   PetscCall(PetscObjectObjectTypeCompare((PetscObject)X, (PetscObject)Y, &sametype));
    121:   PetscCall(PetscLogEventBegin(MAT_AXPY, Y, 0, 0, 0));
    122:   if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) {
    123:     PetscUseTypeMethod(Y, axpy, a, X, str);
    124:   } else {
    125:     PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &transpose, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
    126:     if (transpose) {
    127:       PetscCall(MatTransposeAXPY_Private(Y, a, X, str, X));
    128:     } else {
    129:       PetscCall(PetscObjectTypeCompareAny((PetscObject)Y, &transpose, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
    130:       if (transpose) {
    131:         PetscCall(MatTransposeAXPY_Private(Y, a, X, str, Y));
    132:       } else {
    133:         PetscCall(MatAXPY_BasicWithTypeCompare(Y, a, X, str));
    134:       }
    135:     }
    136:   }
    137:   PetscCall(PetscLogEventEnd(MAT_AXPY, Y, 0, 0, 0));
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B)
    142: {
    143:   PetscErrorCode (*preall)(Mat, Mat, Mat *) = NULL;

    145:   PetscFunctionBegin;
    146:   /* look for any available faster alternative to the general preallocator */
    147:   PetscCall(PetscObjectQueryFunction((PetscObject)Y, "MatAXPYGetPreallocation_C", &preall));
    148:   if (preall) {
    149:     PetscCall((*preall)(Y, X, B));
    150:   } else { /* Use MatPrellocator, assumes same row-col distribution */
    151:     Mat      preallocator;
    152:     PetscInt r, rstart, rend;
    153:     PetscInt m, n, M, N;

    155:     PetscCall(MatGetRowUpperTriangular(Y));
    156:     PetscCall(MatGetRowUpperTriangular(X));
    157:     PetscCall(MatGetSize(Y, &M, &N));
    158:     PetscCall(MatGetLocalSize(Y, &m, &n));
    159:     PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &preallocator));
    160:     PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
    161:     PetscCall(MatSetLayouts(preallocator, Y->rmap, Y->cmap));
    162:     PetscCall(MatSetUp(preallocator));
    163:     PetscCall(MatGetOwnershipRange(preallocator, &rstart, &rend));
    164:     for (r = rstart; r < rend; ++r) {
    165:       PetscInt           ncols;
    166:       const PetscInt    *row;
    167:       const PetscScalar *vals;

    169:       PetscCall(MatGetRow(Y, r, &ncols, &row, &vals));
    170:       PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, vals, INSERT_VALUES));
    171:       PetscCall(MatRestoreRow(Y, r, &ncols, &row, &vals));
    172:       PetscCall(MatGetRow(X, r, &ncols, &row, &vals));
    173:       PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, vals, INSERT_VALUES));
    174:       PetscCall(MatRestoreRow(X, r, &ncols, &row, &vals));
    175:     }
    176:     PetscCall(MatSetOption(preallocator, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
    177:     PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
    178:     PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
    179:     PetscCall(MatRestoreRowUpperTriangular(Y));
    180:     PetscCall(MatRestoreRowUpperTriangular(X));

    182:     PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), B));
    183:     PetscCall(MatSetType(*B, ((PetscObject)Y)->type_name));
    184:     PetscCall(MatSetLayouts(*B, Y->rmap, Y->cmap));
    185:     PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_FALSE, *B));
    186:     PetscCall(MatDestroy(&preallocator));
    187:   }
    188:   PetscFunctionReturn(PETSC_SUCCESS);
    189: }

    191: PetscErrorCode MatAXPY_Basic(Mat Y, PetscScalar a, Mat X, MatStructure str)
    192: {
    193:   PetscFunctionBegin;
    194:   if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) {
    195:     PetscBool isdense;

    197:     /* no need to preallocate if Y is dense */
    198:     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
    199:     if (isdense) str = SUBSET_NONZERO_PATTERN;
    200:   }
    201:   if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) {
    202:     PetscInt           i, start, end, j, ncols, m, n;
    203:     const PetscInt    *row;
    204:     PetscScalar       *val;
    205:     const PetscScalar *vals;

    207:     PetscCall(MatGetSize(X, &m, &n));
    208:     PetscCall(MatGetOwnershipRange(X, &start, &end));
    209:     PetscCall(MatGetRowUpperTriangular(X));
    210:     if (a == 1.0) {
    211:       for (i = start; i < end; i++) {
    212:         PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
    213:         PetscCall(MatSetValues(Y, 1, &i, ncols, row, vals, ADD_VALUES));
    214:         PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
    215:       }
    216:     } else {
    217:       PetscInt vs = 100;
    218:       /* realloc if needed, as this function may be used in parallel */
    219:       PetscCall(PetscMalloc1(vs, &val));
    220:       for (i = start; i < end; i++) {
    221:         PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
    222:         if (vs < ncols) {
    223:           vs = PetscMin(2 * ncols, n);
    224:           PetscCall(PetscRealloc(vs * sizeof(*val), &val));
    225:         }
    226:         for (j = 0; j < ncols; j++) val[j] = a * vals[j];
    227:         PetscCall(MatSetValues(Y, 1, &i, ncols, row, val, ADD_VALUES));
    228:         PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
    229:       }
    230:       PetscCall(PetscFree(val));
    231:     }
    232:     PetscCall(MatRestoreRowUpperTriangular(X));
    233:     PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
    234:     PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
    235:   } else {
    236:     Mat B;

    238:     PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
    239:     PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
    240:     PetscCall(MatHeaderMerge(Y, &B));
    241:   }
    242:   PetscFunctionReturn(PETSC_SUCCESS);
    243: }

    245: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B, Mat Y, PetscScalar a, Mat X, MatStructure str)
    246: {
    247:   PetscInt           i, start, end, j, ncols, m, n;
    248:   const PetscInt    *row;
    249:   PetscScalar       *val;
    250:   const PetscScalar *vals;

    252:   PetscFunctionBegin;
    253:   PetscCall(MatGetSize(X, &m, &n));
    254:   PetscCall(MatGetOwnershipRange(X, &start, &end));
    255:   PetscCall(MatGetRowUpperTriangular(Y));
    256:   PetscCall(MatGetRowUpperTriangular(X));
    257:   if (a == 1.0) {
    258:     for (i = start; i < end; i++) {
    259:       PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
    260:       PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
    261:       PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));

    263:       PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
    264:       PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
    265:       PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
    266:     }
    267:   } else {
    268:     PetscInt vs = 100;
    269:     /* realloc if needed, as this function may be used in parallel */
    270:     PetscCall(PetscMalloc1(vs, &val));
    271:     for (i = start; i < end; i++) {
    272:       PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
    273:       PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
    274:       PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));

    276:       PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
    277:       if (vs < ncols) {
    278:         vs = PetscMin(2 * ncols, n);
    279:         PetscCall(PetscRealloc(vs * sizeof(*val), &val));
    280:       }
    281:       for (j = 0; j < ncols; j++) val[j] = a * vals[j];
    282:       PetscCall(MatSetValues(B, 1, &i, ncols, row, val, ADD_VALUES));
    283:       PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
    284:     }
    285:     PetscCall(PetscFree(val));
    286:   }
    287:   PetscCall(MatRestoreRowUpperTriangular(Y));
    288:   PetscCall(MatRestoreRowUpperTriangular(X));
    289:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    290:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    291:   PetscFunctionReturn(PETSC_SUCCESS);
    292: }

    294: /*@
    295:   MatShift - Computes `Y =  Y + a I`, where `a` is a `PetscScalar`

    297:   Neighbor-wise Collective

    299:   Input Parameters:
    300: + Y - the matrix
    301: - a - the `PetscScalar`

    303:   Level: intermediate

    305:   Notes:
    306:   If `Y` is a rectangular matrix, the shift is done on the main diagonal of the matrix (https://en.wikipedia.org/wiki/Main_diagonal)

    308:   If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
    309:   fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
    310:   entry). No operation is performed when a is zero.

    312:   To form Y = Y + diag(V) use `MatDiagonalSet()`

    314: .seealso: [](ch_matrices), `Mat`, `MatDiagonalSet()`, `MatScale()`, `MatDiagonalScale()`
    315:  @*/
    316: PetscErrorCode MatShift(Mat Y, PetscScalar a)
    317: {
    318:   PetscFunctionBegin;
    320:   PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    321:   PetscCheck(!Y->factortype, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
    322:   MatCheckPreallocated(Y, 1);
    323:   if (a == 0.0) PetscFunctionReturn(PETSC_SUCCESS);

    325:   if (Y->ops->shift) PetscUseTypeMethod(Y, shift, a);
    326:   else PetscCall(MatShift_Basic(Y, a));

    328:   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
    329:   PetscFunctionReturn(PETSC_SUCCESS);
    330: }

    332: PetscErrorCode MatDiagonalSet_Default(Mat Y, Vec D, InsertMode is)
    333: {
    334:   PetscInt           i, start, end;
    335:   const PetscScalar *v;

    337:   PetscFunctionBegin;
    338:   PetscCall(MatGetOwnershipRange(Y, &start, &end));
    339:   PetscCall(VecGetArrayRead(D, &v));
    340:   for (i = start; i < end; i++) PetscCall(MatSetValues(Y, 1, &i, 1, &i, v + i - start, is));
    341:   PetscCall(VecRestoreArrayRead(D, &v));
    342:   PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
    343:   PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
    344:   PetscFunctionReturn(PETSC_SUCCESS);
    345: }

    347: /*@
    348:   MatDiagonalSet - Computes `Y` = `Y` + `D`, where `D` is a diagonal matrix
    349:   that is represented as a vector. Or Y[i,i] = D[i] if `InsertMode` is
    350:   `INSERT_VALUES`.

    352:   Neighbor-wise Collective

    354:   Input Parameters:
    355: + Y  - the input matrix
    356: . D  - the diagonal matrix, represented as a vector
    357: - is - `INSERT_VALUES` or `ADD_VALUES`

    359:   Level: intermediate

    361:   Note:
    362:   If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
    363:   fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
    364:   entry).

    366: .seealso: [](ch_matrices), `Mat`, `MatShift()`, `MatScale()`, `MatDiagonalScale()`
    367: @*/
    368: PetscErrorCode MatDiagonalSet(Mat Y, Vec D, InsertMode is)
    369: {
    370:   PetscInt matlocal, veclocal;

    372:   PetscFunctionBegin;
    375:   MatCheckPreallocated(Y, 1);
    376:   PetscCall(MatGetLocalSize(Y, &matlocal, NULL));
    377:   PetscCall(VecGetLocalSize(D, &veclocal));
    378:   PetscCheck(matlocal == veclocal, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number local rows of matrix %" PetscInt_FMT " does not match that of vector for diagonal %" PetscInt_FMT, matlocal, veclocal);
    379:   if (Y->ops->diagonalset) PetscUseTypeMethod(Y, diagonalset, D, is);
    380:   else PetscCall(MatDiagonalSet_Default(Y, D, is));
    381:   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
    382:   PetscFunctionReturn(PETSC_SUCCESS);
    383: }

    385: /*@
    386:   MatAYPX - Computes Y = a*Y + X.

    388:   Logically Collective

    390:   Input Parameters:
    391: + a   - the `PetscScalar` multiplier
    392: . Y   - the first matrix
    393: . X   - the second matrix
    394: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)

    396:   Level: intermediate

    398: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`
    399:  @*/
    400: PetscErrorCode MatAYPX(Mat Y, PetscScalar a, Mat X, MatStructure str)
    401: {
    402:   PetscFunctionBegin;
    403:   PetscCall(MatScale(Y, a));
    404:   PetscCall(MatAXPY(Y, 1.0, X, str));
    405:   PetscFunctionReturn(PETSC_SUCCESS);
    406: }

    408: /*@
    409:   MatComputeOperator - Computes the explicit matrix

    411:   Collective

    413:   Input Parameters:
    414: + inmat   - the matrix
    415: - mattype - the matrix type for the explicit operator

    417:   Output Parameter:
    418: . mat - the explicit  operator

    420:   Level: advanced

    422:   Note:
    423:   This computation is done by applying the operator to columns of the identity matrix.
    424:   This routine is costly in general, and is recommended for use only with relatively small systems.
    425:   Currently, this routine uses a dense matrix format if `mattype` == `NULL`.

    427: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperatorTranspose()`
    428: @*/
    429: PetscErrorCode MatComputeOperator(Mat inmat, MatType mattype, Mat *mat)
    430: {
    431:   PetscFunctionBegin;
    433:   PetscAssertPointer(mat, 3);
    434:   PetscCall(MatConvert_Shell(inmat, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
    435:   PetscFunctionReturn(PETSC_SUCCESS);
    436: }

    438: /*@
    439:   MatComputeOperatorTranspose - Computes the explicit matrix representation of
    440:   a give matrix that can apply `MatMultTranspose()`

    442:   Collective

    444:   Input Parameters:
    445: + inmat   - the matrix
    446: - mattype - the matrix type for the explicit operator

    448:   Output Parameter:
    449: . mat - the explicit  operator transposed

    451:   Level: advanced

    453:   Note:
    454:   This computation is done by applying the transpose of the operator to columns of the identity matrix.
    455:   This routine is costly in general, and is recommended for use only with relatively small systems.
    456:   Currently, this routine uses a dense matrix format if `mattype` == `NULL`.

    458: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperator()`
    459: @*/
    460: PetscErrorCode MatComputeOperatorTranspose(Mat inmat, MatType mattype, Mat *mat)
    461: {
    462:   Mat A;

    464:   PetscFunctionBegin;
    466:   PetscAssertPointer(mat, 3);
    467:   PetscCall(MatCreateTranspose(inmat, &A));
    468:   PetscCall(MatConvert_Shell(A, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
    469:   PetscCall(MatDestroy(&A));
    470:   PetscFunctionReturn(PETSC_SUCCESS);
    471: }

    473: /*@
    474:   MatFilter - Set all values in the matrix with an absolute value less than or equal to the tolerance to zero, and optionally compress the underlying storage

    476:   Input Parameters:
    477: + A        - The matrix
    478: . tol      - The zero tolerance
    479: . compress - Whether the storage from the input matrix `A` should be compressed once values less than or equal to `tol` are set to zero
    480: - keep     - If `compress` is true and for a given row of `A`, the diagonal coefficient is less than or equal to `tol`, indicates whether it should be left in the structure or eliminated as well

    482:   Level: intermediate

    484: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatZeroEntries()`, `MatEliminateZeros()`, `VecFilter()`
    485:  @*/
    486: PetscErrorCode MatFilter(Mat A, PetscReal tol, PetscBool compress, PetscBool keep)
    487: {
    488:   Mat          a;
    489:   PetscScalar *newVals;
    490:   PetscInt    *newCols, rStart, rEnd, maxRows, r, colMax = 0, nnz0 = 0, nnz1 = 0;
    491:   PetscBool    flg;

    493:   PetscFunctionBegin;
    494:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, ""));
    495:   if (flg) {
    496:     PetscCall(MatDenseGetLocalMatrix(A, &a));
    497:     PetscCall(MatDenseGetLDA(a, &r));
    498:     PetscCall(MatGetSize(a, &rStart, &rEnd));
    499:     PetscCall(MatDenseGetArray(a, &newVals));
    500:     for (; colMax < rEnd; ++colMax) {
    501:       for (maxRows = 0; maxRows < rStart; ++maxRows) newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) <= tol ? 0.0 : newVals[maxRows + colMax * r];
    502:     }
    503:     PetscCall(MatDenseRestoreArray(a, &newVals));
    504:   } else {
    505:     const PetscInt *ranges;
    506:     PetscMPIInt     rank, size;

    508:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
    509:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    510:     PetscCall(MatGetOwnershipRanges(A, &ranges));
    511:     rStart = ranges[rank];
    512:     rEnd   = ranges[rank + 1];
    513:     PetscCall(MatGetRowUpperTriangular(A));
    514:     for (r = rStart; r < rEnd; ++r) {
    515:       PetscInt ncols;

    517:       PetscCall(MatGetRow(A, r, &ncols, NULL, NULL));
    518:       colMax = PetscMax(colMax, ncols);
    519:       PetscCall(MatRestoreRow(A, r, &ncols, NULL, NULL));
    520:     }
    521:     maxRows = 0;
    522:     for (r = 0; r < size; ++r) maxRows = PetscMax(maxRows, ranges[r + 1] - ranges[r]);
    523:     PetscCall(PetscCalloc2(colMax, &newCols, colMax, &newVals));
    524:     PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &flg)); /* cache user-defined value */
    525:     PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
    526:     /* short-circuit code in MatAssemblyBegin() and MatAssemblyEnd()             */
    527:     /* that are potentially called many times depending on the distribution of A */
    528:     for (r = rStart; r < rStart + maxRows; ++r) {
    529:       if (r < rEnd) {
    530:         const PetscScalar *vals;
    531:         const PetscInt    *cols;
    532:         PetscInt           ncols, newcols = 0, c;

    534:         PetscCall(MatGetRow(A, r, &ncols, &cols, &vals));
    535:         nnz0 += ncols - 1;
    536:         for (c = 0; c < ncols; ++c) {
    537:           if (PetscUnlikely(PetscAbsScalar(vals[c]) <= tol)) newCols[newcols++] = cols[c];
    538:         }
    539:         nnz1 += ncols - newcols - 1;
    540:         PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals));
    541:         PetscCall(MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES));
    542:       }
    543:       PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    544:       PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    545:     }
    546:     PetscCall(MatRestoreRowUpperTriangular(A));
    547:     PetscCall(PetscFree2(newCols, newVals));
    548:     PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, flg)); /* reset option to its user-defined value */
    549:     if (nnz0 > 0) PetscCall(PetscInfo(NULL, "Filtering left %g%% edges in graph\n", 100 * (double)nnz1 / (double)nnz0));
    550:     else PetscCall(PetscInfo(NULL, "Warning: %" PetscInt_FMT " edges to filter with %" PetscInt_FMT " rows\n", nnz0, maxRows));
    551:   }
    552:   if (compress && A->ops->eliminatezeros) {
    553:     Mat       B;
    554:     PetscBool flg;

    556:     PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
    557:     if (!flg) {
    558:       PetscCall(MatEliminateZeros(A, keep));
    559:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
    560:       PetscCall(MatHeaderReplace(A, &B));
    561:     }
    562:   }
    563:   PetscFunctionReturn(PETSC_SUCCESS);
    564: }