Actual source code: schurm.c

      1: #include <../src/ksp/ksp/utils/schurm/schurm.h>

      3: const char *const MatSchurComplementAinvTypes[] = {"DIAG", "LUMP", "BLOCKDIAG", "FULL", "MatSchurComplementAinvType", "MAT_SCHUR_COMPLEMENT_AINV_", NULL};

      5: PetscErrorCode MatCreateVecs_SchurComplement(Mat N, Vec *right, Vec *left)
      6: {
      7:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

      9:   PetscFunctionBegin;
     10:   if (Na->D) {
     11:     PetscCall(MatCreateVecs(Na->D, right, left));
     12:     PetscFunctionReturn(PETSC_SUCCESS);
     13:   }
     14:   if (right) PetscCall(MatCreateVecs(Na->B, right, NULL));
     15:   if (left) PetscCall(MatCreateVecs(Na->C, NULL, left));
     16:   PetscFunctionReturn(PETSC_SUCCESS);
     17: }

     19: PetscErrorCode MatView_SchurComplement(Mat N, PetscViewer viewer)
     20: {
     21:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

     23:   PetscFunctionBegin;
     24:   PetscCall(PetscViewerASCIIPrintf(viewer, "Schur complement A11 - A10 inv(A00) A01\n"));
     25:   if (Na->D) {
     26:     PetscCall(PetscViewerASCIIPrintf(viewer, "A11\n"));
     27:     PetscCall(PetscViewerASCIIPushTab(viewer));
     28:     PetscCall(MatView(Na->D, viewer));
     29:     PetscCall(PetscViewerASCIIPopTab(viewer));
     30:   } else {
     31:     PetscCall(PetscViewerASCIIPrintf(viewer, "A11 = 0\n"));
     32:   }
     33:   PetscCall(PetscViewerASCIIPrintf(viewer, "A10\n"));
     34:   PetscCall(PetscViewerASCIIPushTab(viewer));
     35:   PetscCall(MatView(Na->C, viewer));
     36:   PetscCall(PetscViewerASCIIPopTab(viewer));
     37:   PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block viewable with the additional option -%sksp_view\n", ((PetscObject)Na->ksp)->prefix ? ((PetscObject)Na->ksp)->prefix : NULL));
     38:   PetscCall(PetscViewerASCIIPrintf(viewer, "A01\n"));
     39:   PetscCall(PetscViewerASCIIPushTab(viewer));
     40:   PetscCall(MatView(Na->B, viewer));
     41:   PetscCall(PetscViewerASCIIPopTab(viewer));
     42:   PetscFunctionReturn(PETSC_SUCCESS);
     43: }

     45: /*
     46:            A11^T - A01^T ksptrans(A00,Ap00) A10^T
     47: */
     48: PetscErrorCode MatMultTranspose_SchurComplement(Mat N, Vec x, Vec y)
     49: {
     50:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

     52:   PetscFunctionBegin;
     53:   if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
     54:   if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
     55:   PetscCall(MatMultTranspose(Na->C, x, Na->work1));
     56:   PetscCall(KSPSolveTranspose(Na->ksp, Na->work1, Na->work2));
     57:   PetscCall(MatMultTranspose(Na->B, Na->work2, y));
     58:   PetscCall(VecScale(y, -1.0));
     59:   if (Na->D) PetscCall(MatMultTransposeAdd(Na->D, x, y, y));
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: /*
     64:            A11 - A10 ksp(A00,Ap00) A01
     65: */
     66: PetscErrorCode MatMult_SchurComplement(Mat N, Vec x, Vec y)
     67: {
     68:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

     70:   PetscFunctionBegin;
     71:   if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
     72:   if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
     73:   PetscCall(MatMult(Na->B, x, Na->work1));
     74:   PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
     75:   PetscCall(MatMult(Na->C, Na->work2, y));
     76:   PetscCall(VecScale(y, -1.0));
     77:   if (Na->D) PetscCall(MatMultAdd(Na->D, x, y, y));
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: /*
     82:            A11 - A10 ksp(A00,Ap00) A01
     83: */
     84: PetscErrorCode MatMultAdd_SchurComplement(Mat N, Vec x, Vec y, Vec z)
     85: {
     86:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

     88:   PetscFunctionBegin;
     89:   if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
     90:   if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
     91:   PetscCall(MatMult(Na->B, x, Na->work1));
     92:   PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
     93:   if (y == z) {
     94:     PetscCall(VecScale(Na->work2, -1.0));
     95:     PetscCall(MatMultAdd(Na->C, Na->work2, z, z));
     96:   } else {
     97:     PetscCall(MatMult(Na->C, Na->work2, z));
     98:     PetscCall(VecAYPX(z, -1.0, y));
     99:   }
    100:   if (Na->D) PetscCall(MatMultAdd(Na->D, x, z, z));
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: PetscErrorCode MatSetFromOptions_SchurComplement(Mat N, PetscOptionItems *PetscOptionsObject)
    105: {
    106:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

    108:   PetscFunctionBegin;
    109:   PetscOptionsHeadBegin(PetscOptionsObject, "MatSchurComplementOptions");
    110:   Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
    111:   PetscCall(PetscOptionsEnum("-mat_schur_complement_ainv_type", "Type of approximation for DIAGFORM(A00) used when assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01", "MatSchurComplementSetAinvType", MatSchurComplementAinvTypes, (PetscEnum)Na->ainvtype,
    112:                              (PetscEnum *)&Na->ainvtype, NULL));
    113:   PetscOptionsHeadEnd();
    114:   PetscCall(KSPSetFromOptions(Na->ksp));
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: PetscErrorCode MatDestroy_SchurComplement(Mat N)
    119: {
    120:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

    122:   PetscFunctionBegin;
    123:   PetscCall(MatDestroy(&Na->A));
    124:   PetscCall(MatDestroy(&Na->Ap));
    125:   PetscCall(MatDestroy(&Na->B));
    126:   PetscCall(MatDestroy(&Na->C));
    127:   PetscCall(MatDestroy(&Na->D));
    128:   PetscCall(VecDestroy(&Na->work1));
    129:   PetscCall(VecDestroy(&Na->work2));
    130:   PetscCall(KSPDestroy(&Na->ksp));
    131:   PetscCall(PetscFree(N->data));
    132:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", NULL));
    133:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", NULL));
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: /*@
    138:   MatCreateSchurComplement - Creates a new `Mat` that behaves like the Schur complement of a matrix

    140:   Collective

    142:   Input Parameters:
    143: + A00  - the upper-left block of the original matrix $A = [A00 A01; A10 A11]$
    144: . Ap00 - preconditioning matrix for use in $ksp(A00,Ap00)$ to approximate the action of $A00^{-1}$
    145: . A01  - the upper-right block of the original matrix $A = [A00 A01; A10 A11]$
    146: . A10  - the lower-left block of the original matrix $A = [A00 A01; A10 A11]$
    147: - A11  - (optional) the lower-right block of the original matrix $A = [A00 A01; A10 A11]$

    149:   Output Parameter:
    150: . S - the matrix that behaves as the Schur complement $S = A11 - A10 ksp(A00,Ap00) A01$

    152:   Level: intermediate

    154:   Notes:
    155:   The Schur complement is NOT explicitly formed! Rather, this function returns a virtual Schur complement
    156:   that can compute the matrix-vector product by using formula $S = A11 - A10 A^{-1} A01$
    157:   for Schur complement `S` and a `KSP` solver to approximate the action of $A^{-1}$.

    159:   All four matrices must have the same MPI communicator.

    161:   `A00` and  `A11` must be square matrices.

    163:   `MatGetSchurComplement()` takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
    164:   a sparse approximation to the Schur complement (useful for building a preconditioner for the Schur complement) which can be obtained from this
    165:   matrix with `MatSchurComplementGetPmat()`

    167:   Developer Notes:
    168:   The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
    169:   remove redundancy and be clearer and simpler.

    171: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatGetSchurComplement()`,
    172:           `MatSchurComplementGetPmat()`, `MatSchurComplementSetSubMatrices()`
    173: @*/
    174: PetscErrorCode MatCreateSchurComplement(Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11, Mat *S)
    175: {
    176:   PetscFunctionBegin;
    177:   PetscCall(KSPInitializePackage());
    178:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A00), S));
    179:   PetscCall(MatSetType(*S, MATSCHURCOMPLEMENT));
    180:   PetscCall(MatSchurComplementSetSubMatrices(*S, A00, Ap00, A01, A10, A11));
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: /*@
    185:   MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement

    187:   Collective

    189:   Input Parameters:
    190: + S    - matrix obtained with `MatSetType`(S,`MATSCHURCOMPLEMENT`)
    191: . A00  - the upper-left block of the original matrix $A = [A00 A01; A10 A11]$
    192: . Ap00 - preconditioning matrix for use in $ksp(A00,Ap00)$ to approximate the action of $A00^{-1}$
    193: . A01  - the upper-right block of the original matrix $A = [A00 A01; A10 A11]$
    194: . A10  - the lower-left block of the original matrix $A = [A00 A01; A10 A11]$
    195: - A11  - (optional) the lower-right block of the original matrix $A = [A00 A01; A10 A11]$

    197:   Level: intermediate

    199:   Notes:
    200:   The Schur complement is NOT explicitly formed! Rather, this
    201:   object performs the matrix-vector product of the Schur complement by using formula $S = A11 - A10 ksp(A00,Ap00) A01$

    203:   All four matrices must have the same MPI communicator.

    205:   `A00` and `A11` must be square matrices.

    207:   This is to be used in the context of code such as
    208: .vb
    209:      MatSetType(S,MATSCHURCOMPLEMENT);
    210:      MatSchurComplementSetSubMatrices(S,...);
    211: .ve
    212:   while `MatSchurComplementUpdateSubMatrices()` should only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`

    214: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`
    215: @*/
    216: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
    217: {
    218:   Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
    219:   PetscBool            isschur;

    221:   PetscFunctionBegin;
    222:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    223:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    224:   PetscCheck(!S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementUpdateSubMatrices() for already used matrix");
    229:   PetscCheckSameComm(A00, 2, Ap00, 3);
    230:   PetscCheckSameComm(A00, 2, A01, 4);
    231:   PetscCheckSameComm(A00, 2, A10, 5);
    232:   PetscCheck(A00->rmap->n == A00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, A00->rmap->n, A00->cmap->n);
    233:   PetscCheck(A00->rmap->n == Ap00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local rows of Ap00 %" PetscInt_FMT, A00->rmap->n, Ap00->rmap->n);
    234:   PetscCheck(Ap00->rmap->n == Ap00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of Ap00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, Ap00->rmap->n, Ap00->cmap->n);
    235:   PetscCheck(A00->cmap->n == A01->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A00 %" PetscInt_FMT " do not equal local rows of A01 %" PetscInt_FMT, A00->cmap->n, A01->rmap->n);
    236:   PetscCheck(A10->cmap->n == A00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A10 %" PetscInt_FMT " do not equal local rows of A00 %" PetscInt_FMT, A10->cmap->n, A00->rmap->n);
    237:   if (A11) {
    239:     PetscCheckSameComm(A00, 2, A11, 6);
    240:     PetscCheck(A10->rmap->n == A11->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A10 %" PetscInt_FMT " do not equal local rows A11 %" PetscInt_FMT, A10->rmap->n, A11->rmap->n);
    241:   }

    243:   PetscCall(MatSetSizes(S, A10->rmap->n, A01->cmap->n, A10->rmap->N, A01->cmap->N));
    244:   PetscCall(PetscObjectReference((PetscObject)A00));
    245:   PetscCall(PetscObjectReference((PetscObject)Ap00));
    246:   PetscCall(PetscObjectReference((PetscObject)A01));
    247:   PetscCall(PetscObjectReference((PetscObject)A10));
    248:   Na->A  = A00;
    249:   Na->Ap = Ap00;
    250:   Na->B  = A01;
    251:   Na->C  = A10;
    252:   Na->D  = A11;
    253:   if (A11) PetscCall(PetscObjectReference((PetscObject)A11));
    254:   PetscCall(MatSetUp(S));
    255:   PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
    256:   S->assembled = PETSC_TRUE;
    257:   PetscFunctionReturn(PETSC_SUCCESS);
    258: }

    260: /*@
    261:   MatSchurComplementGetKSP - Gets the `KSP` object that is used to solve with `A00` in the Schur complement matrix $S = A11 - A10 ksp(A00,Ap00) A01$

    263:   Not Collective

    265:   Input Parameter:
    266: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of $ A11 - A10 ksp(A00,Ap00) A01 $

    268:   Output Parameter:
    269: . ksp - the linear solver object

    271:   Options Database Key:
    272: . -fieldsplit_<splitname_0>_XXX - sets `KSP` and `PC` options for the 0-split solver inside the Schur complement used in `PCFIELDSPLIT`; default <splitname_0> is 0.

    274:   Level: intermediate

    276: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementSetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`
    277: @*/
    278: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
    279: {
    280:   Mat_SchurComplement *Na;
    281:   PetscBool            isschur;

    283:   PetscFunctionBegin;
    285:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    286:   PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
    287:   PetscAssertPointer(ksp, 2);
    288:   Na   = (Mat_SchurComplement *)S->data;
    289:   *ksp = Na->ksp;
    290:   PetscFunctionReturn(PETSC_SUCCESS);
    291: }

    293: /*@
    294:   MatSchurComplementSetKSP - Sets the `KSP` object that is used to solve with `A00` in the Schur complement matrix $ S = A11 - A10 ksp(A00,Ap00) A01$

    296:   Not Collective

    298:   Input Parameters:
    299: + S   - matrix created with `MatCreateSchurComplement()`
    300: - ksp - the linear solver object

    302:   Level: developer

    304:   Developer Notes:
    305:   This is used in `PCFIELDSPLIT` to reuse the 0-split `KSP` to implement $ksp(A00,Ap00)$ in `S`.
    306:   The `KSP` operators are overwritten with `A00` and `Ap00` currently set in `S`.

    308: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MATSCHURCOMPLEMENT`
    309: @*/
    310: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
    311: {
    312:   Mat_SchurComplement *Na;
    313:   PetscBool            isschur;

    315:   PetscFunctionBegin;
    317:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    318:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    320:   Na = (Mat_SchurComplement *)S->data;
    321:   PetscCall(PetscObjectReference((PetscObject)ksp));
    322:   PetscCall(KSPDestroy(&Na->ksp));
    323:   Na->ksp = ksp;
    324:   PetscCall(KSPSetOperators(Na->ksp, Na->A, Na->Ap));
    325:   PetscFunctionReturn(PETSC_SUCCESS);
    326: }

    328: /*@
    329:   MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

    331:   Collective

    333:   Input Parameters:
    334: + S    - matrix obtained with `MatCreateSchurComplement()` (or `MatSchurSetSubMatrices()`) and implementing the action of $A11 - A10 ksp(A00,Ap00) A01$
    335: . A00  - the upper-left block of the original matrix $A = [A00 A01; A10 A11]$
    336: . Ap00 - preconditioning matrix for use in $ksp(A00,Ap00)$ to approximate the action of $A00^{-1}$
    337: . A01  - the upper-right block of the original matrix $A = [A00 A01; A10 A11]$
    338: . A10  - the lower-left block of the original matrix $A = [A00 A01; A10 A11]$
    339: - A11  - (optional) the lower-right block of the original matrix $A = [A00 A01; A10 A11]$

    341:   Level: intermediate

    343:   Notes:
    344:   All four matrices must have the same MPI communicator

    346:   `A00` and  `A11` must be square matrices

    348:   All of the matrices provided must have the same sizes as was used with `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`
    349:   though they need not be the same matrices.

    351:   This can only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`, it cannot be called immediately after `MatSetType`(S,`MATSCHURCOMPLEMENT`);

    353:   Developer Notes:
    354:   This code is almost identical to `MatSchurComplementSetSubMatrices()`. The API should be refactored.

    356: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`
    357: @*/
    358: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
    359: {
    360:   Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
    361:   PetscBool            isschur;

    363:   PetscFunctionBegin;
    365:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    366:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    367:   PetscCheck(S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementSetSubMatrices() for a new matrix");
    372:   PetscCheckSameComm(A00, 2, Ap00, 3);
    373:   PetscCheckSameComm(A00, 2, A01, 4);
    374:   PetscCheckSameComm(A00, 2, A10, 5);
    375:   PetscCheck(A00->rmap->n == A00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, A00->rmap->n, A00->cmap->n);
    376:   PetscCheck(A00->rmap->n == Ap00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local rows of Ap00 %" PetscInt_FMT, A00->rmap->n, Ap00->rmap->n);
    377:   PetscCheck(Ap00->rmap->n == Ap00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of Ap00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, Ap00->rmap->n, Ap00->cmap->n);
    378:   PetscCheck(A00->cmap->n == A01->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A00 %" PetscInt_FMT " do not equal local rows of A01 %" PetscInt_FMT, A00->cmap->n, A01->rmap->n);
    379:   PetscCheck(A10->cmap->n == A00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A10 %" PetscInt_FMT " do not equal local rows of A00 %" PetscInt_FMT, A10->cmap->n, A00->rmap->n);
    380:   if (A11) {
    382:     PetscCheckSameComm(A00, 2, A11, 6);
    383:     PetscCheck(A10->rmap->n == A11->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A10 %" PetscInt_FMT " do not equal local rows A11 %" PetscInt_FMT, A10->rmap->n, A11->rmap->n);
    384:   }

    386:   PetscCall(PetscObjectReference((PetscObject)A00));
    387:   PetscCall(PetscObjectReference((PetscObject)Ap00));
    388:   PetscCall(PetscObjectReference((PetscObject)A01));
    389:   PetscCall(PetscObjectReference((PetscObject)A10));
    390:   if (A11) PetscCall(PetscObjectReference((PetscObject)A11));

    392:   PetscCall(MatDestroy(&Na->A));
    393:   PetscCall(MatDestroy(&Na->Ap));
    394:   PetscCall(MatDestroy(&Na->B));
    395:   PetscCall(MatDestroy(&Na->C));
    396:   PetscCall(MatDestroy(&Na->D));

    398:   Na->A  = A00;
    399:   Na->Ap = Ap00;
    400:   Na->B  = A01;
    401:   Na->C  = A10;
    402:   Na->D  = A11;

    404:   PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
    405:   PetscFunctionReturn(PETSC_SUCCESS);
    406: }

    408: /*@
    409:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

    411:   Collective

    413:   Input Parameter:
    414: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of $A11 - A10 ksp(A00,Ap00) A01$

    416:   Output Parameters:
    417: + A00  - the upper-left block of the original matrix $A = [A00 A01; A10 A11]$
    418: . Ap00 - preconditioning matrix for use in $ksp(A00,Ap00)$ to approximate the action of $A^{-1}$
    419: . A01  - the upper-right block of the original matrix $A = [A00 A01; A10 A11]$
    420: . A10  - the lower-left block of the original matrix $A = [A00 A01; A10 A11]$
    421: - A11  - (optional) the lower-right block of the original matrix $A = [A00 A01; A10 A11]$

    423:   Level: intermediate

    425:   Note:
    426:   Use `NULL` for any unneeded output argument.

    428:   The reference counts of the submatrices are not increased before they are returned and the matrices should not be modified or destroyed.

    430: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`
    431: @*/
    432: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S, Mat *A00, Mat *Ap00, Mat *A01, Mat *A10, Mat *A11)
    433: {
    434:   Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
    435:   PetscBool            flg;

    437:   PetscFunctionBegin;
    439:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &flg));
    440:   PetscCheck(flg, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
    441:   if (A00) *A00 = Na->A;
    442:   if (Ap00) *Ap00 = Na->Ap;
    443:   if (A01) *A01 = Na->B;
    444:   if (A10) *A10 = Na->C;
    445:   if (A11) *A11 = Na->D;
    446:   PetscFunctionReturn(PETSC_SUCCESS);
    447: }

    449: #include <petsc/private/kspimpl.h>

    451: /*@
    452:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

    454:   Collective

    456:   Input Parameter:
    457: . A - the matrix obtained with `MatCreateSchurComplement()`

    459:   Output Parameter:
    460: . S - the Schur complement matrix

    462:   Level: advanced

    464:   Notes:
    465:   This can be expensive when `S` is large, so it is mainly for testing

    467:   Use `MatSchurComplementGetPmat()` to get a sparse approximation for the Schur complement suitable for use in building a preconditioner

    469:   `S` will automatically have the same prefix as `A` appended by `explicit_operator_`,
    470:   there are three options available: `-fieldsplit_1_explicit_operator_mat_type`,
    471:   `-fieldsplit_1_explicit_operator_mat_symmetric`, and `-fieldsplit_1_explicit_operator_mat_hermitian`

    473:   Developer Note:
    474:   The three aforementioned should not be parsed and used in this routine, but rather in `MatSetFromOptions()`

    476: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`, `MatSchurComplementGetPmat()`
    477: @*/
    478: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
    479: {
    480:   Mat       P, B, C, D, E = NULL, Bd, AinvBd, sub = NULL;
    481:   MatType   mtype;
    482:   VecType   vtype;
    483:   KSP       ksp;
    484:   PetscInt  n, N, m, M;
    485:   PetscBool flg = PETSC_FALSE, set, symm;
    486:   char      prefix[256], type[256];

    488:   PetscFunctionBegin;
    489:   PetscCall(PetscObjectQuery((PetscObject)A, "AinvB", (PetscObject *)&AinvBd));
    490:   set = (PetscBool)(AinvBd != NULL);
    491:   if (set && AinvBd->cmap->N == -1) PetscFunctionReturn(PETSC_SUCCESS); // early bail out if composed Mat is uninitialized
    492:   PetscCall(MatSchurComplementGetSubMatrices(A, &P, NULL, &B, &C, &D));
    493:   PetscCall(MatGetVecType(B, &vtype));
    494:   PetscCall(MatGetLocalSize(B, &m, &n));
    495:   PetscCall(MatSchurComplementGetKSP(A, &ksp));
    496:   PetscCall(KSPSetUp(ksp));
    497:   if (set) {
    498:     PetscCheck(AinvBd->cmap->N >= A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Composed Mat should have at least as many columns as the Schur complement (%" PetscInt_FMT " >= %" PetscInt_FMT ")", AinvBd->cmap->N, A->cmap->N);
    499:     PetscCall(MatGetType(AinvBd, &mtype));
    500:     if (AinvBd->cmap->N > A->cmap->N) {
    501:       Mat s[2];

    503:       PetscCall(MatDuplicate(AinvBd, MAT_DO_NOT_COPY_VALUES, &Bd));
    504:       PetscCall(MatDenseGetSubMatrix(Bd, PETSC_DECIDE, PETSC_DECIDE, A->cmap->N, AinvBd->cmap->N, s));
    505:       PetscCall(MatDenseGetSubMatrix(AinvBd, PETSC_DECIDE, PETSC_DECIDE, A->cmap->N, AinvBd->cmap->N, s + 1));
    506:       PetscCall(MatCopy(s[1], s[0], SAME_NONZERO_PATTERN)); // copy the last columns of the composed Mat, which are likely the input columns of PCApply_FieldSplit_Schur()
    507:       PetscCall(MatDenseRestoreSubMatrix(AinvBd, s + 1));
    508:       PetscCall(MatDenseRestoreSubMatrix(Bd, s));
    509:       PetscCall(MatDenseGetSubMatrix(Bd, PETSC_DECIDE, PETSC_DECIDE, 0, A->cmap->N, &sub));
    510:       PetscCall(MatConvert(B, mtype, MAT_REUSE_MATRIX, &sub)); // copy A01 into the first columns of the block of RHS of KSPMatSolve()
    511:       PetscCall(MatDenseRestoreSubMatrix(Bd, &sub));
    512:     } else PetscCall(MatConvert(B, mtype, MAT_INITIAL_MATRIX, &Bd));
    513:   } else {
    514:     PetscCall(MatGetSize(B, &M, &N));
    515:     PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)A), vtype, m, n, M, N, -1, NULL, &AinvBd));
    516:     PetscCall(MatGetType(AinvBd, &mtype));
    517:     PetscCall(MatConvert(B, mtype, MAT_INITIAL_MATRIX, &Bd));
    518:   }
    519:   PetscCall(KSPMatSolve(ksp, Bd, AinvBd));
    520:   if (set && AinvBd->cmap->N > A->cmap->N) {
    521:     Mat          AinvB;
    522:     PetscScalar *v;
    523:     PetscBool    match;

    525:     PetscCall(PetscObjectTypeCompareAny((PetscObject)AinvBd, &match, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
    526:     if (match) {
    527: #if PetscDefined(HAVE_CUDA)
    528:       PetscCall(MatDenseCUDAGetArrayWrite(AinvBd, &v));
    529:       PetscCall(MatCreateDenseCUDA(PetscObjectComm((PetscObject)A), AinvBd->rmap->n, A->cmap->n, AinvBd->rmap->N, A->cmap->N, v, &AinvB));
    530:       PetscCall(MatDenseCUDAReplaceArray(AinvB, v));
    531:       PetscCall(MatDenseCUDARestoreArrayWrite(AinvBd, &v));
    532: #endif
    533:     } else {
    534:       PetscCall(PetscObjectTypeCompareAny((PetscObject)AinvBd, &match, MATSEQDENSEHIP, MATMPIDENSEHIP, ""));
    535:       if (match) {
    536: #if PetscDefined(HAVE_HIP)
    537:         PetscCall(MatDenseHIPGetArrayWrite(AinvBd, &v));
    538:         PetscCall(MatCreateDenseHIP(PetscObjectComm((PetscObject)A), AinvBd->rmap->n, A->cmap->n, AinvBd->rmap->N, A->cmap->N, v, &AinvB));
    539:         PetscCall(MatDenseHIPReplaceArray(AinvB, v));
    540:         PetscCall(MatDenseHIPRestoreArrayWrite(AinvBd, &v));
    541: #endif
    542:       } else {
    543:         PetscCall(MatDenseGetArrayWrite(AinvBd, &v)); // no easy way to resize a Mat, so create a new one with the same data pointer
    544:         PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), AinvBd->rmap->n, A->cmap->n, AinvBd->rmap->N, A->cmap->N, v, &AinvB));
    545:         PetscCall(MatDenseReplaceArray(AinvB, v)); // let MatDestroy() free the data pointer
    546:         PetscCall(MatDenseRestoreArrayWrite(AinvBd, &v));
    547:       }
    548:     }
    549:     PetscCall(MatHeaderReplace(AinvBd, &AinvB)); // replace the input composed Mat with just A00^-1 A01 (trailing columns are removed)
    550:   }
    551:   PetscCall(MatDestroy(&Bd));
    552:   if (!set) PetscCall(MatFilter(AinvBd, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
    553:   if (D) {
    554:     PetscCall(MatGetLocalSize(D, &m, &n));
    555:     PetscCall(MatGetSize(D, &M, &N));
    556:     PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)A), vtype, m, n, M, N, -1, NULL, S));
    557:   }
    558:   PetscCall(MatMatMult(C, AinvBd, D ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DETERMINE, S));
    559:   if (!set) PetscCall(MatDestroy(&AinvBd));
    560:   else {
    561:     PetscCall(MatScale(AinvBd, -1.0));
    562:     PetscCall(MatFilter(AinvBd, PETSC_MACHINE_EPSILON, PETSC_FALSE, PETSC_FALSE));
    563:     PetscCall(MatFilter(*S, PETSC_MACHINE_EPSILON, PETSC_FALSE, PETSC_FALSE));
    564:   }
    565:   if (D) {
    566:     PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
    567:     if (flg) {
    568:       PetscCall(MatIsSymmetricKnown(A, &set, &symm));
    569:       if (!set || !symm) PetscCall(MatConvert(D, MATBAIJ, MAT_INITIAL_MATRIX, &E)); /* convert the (1,1) block to nonsymmetric storage for MatAXPY() */
    570:     }
    571:     PetscCall(MatAXPY(*S, -1.0, E ? E : D, DIFFERENT_NONZERO_PATTERN));        /* calls Mat[Get|Restore]RowUpperTriangular(), so only the upper triangular part is valid with symmetric storage */
    572:     if (!E && flg) PetscCall(MatConvert(*S, MATSBAIJ, MAT_INPLACE_MATRIX, S)); /* if A is symmetric and the (1,1) block is a MatSBAIJ, return S as a MatSBAIJ since the lower triangular part is invalid */
    573:   }
    574:   PetscCall(MatDestroy(&E));
    575:   PetscCall(MatScale(*S, -1.0));
    576:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sexplicit_operator_", ((PetscObject)A)->prefix ? ((PetscObject)A)->prefix : ""));
    577:   PetscCall(MatSetOptionsPrefix(*S, prefix));
    578:   PetscObjectOptionsBegin((PetscObject)*S);
    579:   PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, !E && flg ? MATSBAIJ : mtype, type, 256, &set));
    580:   if (set) PetscCall(MatConvert(*S, type, MAT_INPLACE_MATRIX, S));
    581:   flg = PETSC_FALSE;
    582:   PetscCall(PetscOptionsBool("-mat_symmetric", "Sets the MAT_SYMMETRIC option", "MatSetOption", flg, &flg, &set));
    583:   if (set) PetscCall(MatSetOption(*S, MAT_SYMMETRIC, flg));
    584:   if (PetscDefined(USE_COMPLEX)) {
    585:     flg = PETSC_FALSE;
    586:     PetscCall(PetscOptionsBool("-mat_hermitian", "Sets the MAT_HERMITIAN option", "MatSetOption", flg, &flg, &set));
    587:     if (set) PetscCall(MatSetOption(*S, MAT_HERMITIAN, flg));
    588:   }
    589:   PetscOptionsEnd();
    590:   PetscFunctionReturn(PETSC_SUCCESS);
    591: }

    593: /* Developer Notes:
    594:     This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
    595: PetscErrorCode MatGetSchurComplement_Basic(Mat mat, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
    596: {
    597:   Mat      A = NULL, Ap = NULL, B = NULL, C = NULL, D = NULL;
    598:   MatReuse reuse;

    600:   PetscFunctionBegin;
    610:   if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);

    614:   PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

    616:   reuse = MAT_INITIAL_MATRIX;
    617:   if (mreuse == MAT_REUSE_MATRIX) {
    618:     PetscCall(MatSchurComplementGetSubMatrices(*S, &A, &Ap, &B, &C, &D));
    619:     PetscCheck(A && Ap && B && C, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Attempting to reuse matrix but Schur complement matrices unset");
    620:     PetscCheck(A == Ap, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Preconditioning matrix does not match operator");
    621:     PetscCall(MatDestroy(&Ap)); /* get rid of extra reference */
    622:     reuse = MAT_REUSE_MATRIX;
    623:   }
    624:   PetscCall(MatCreateSubMatrix(mat, isrow0, iscol0, reuse, &A));
    625:   PetscCall(MatCreateSubMatrix(mat, isrow0, iscol1, reuse, &B));
    626:   PetscCall(MatCreateSubMatrix(mat, isrow1, iscol0, reuse, &C));
    627:   PetscCall(MatCreateSubMatrix(mat, isrow1, iscol1, reuse, &D));
    628:   switch (mreuse) {
    629:   case MAT_INITIAL_MATRIX:
    630:     PetscCall(MatCreateSchurComplement(A, A, B, C, D, S));
    631:     break;
    632:   case MAT_REUSE_MATRIX:
    633:     PetscCall(MatSchurComplementUpdateSubMatrices(*S, A, A, B, C, D));
    634:     break;
    635:   default:
    636:     PetscCheck(mreuse == MAT_IGNORE_MATRIX, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unrecognized value of mreuse %d", (int)mreuse);
    637:   }
    638:   if (preuse != MAT_IGNORE_MATRIX) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, ainvtype, preuse, Sp));
    639:   PetscCall(MatDestroy(&A));
    640:   PetscCall(MatDestroy(&B));
    641:   PetscCall(MatDestroy(&C));
    642:   PetscCall(MatDestroy(&D));
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }

    646: /*@
    647:   MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.

    649:   Collective

    651:   Input Parameters:
    652: + A        - matrix in which the complement is to be taken
    653: . isrow0   - rows to eliminate
    654: . iscol0   - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
    655: . isrow1   - rows in which the Schur complement is formed
    656: . iscol1   - columns in which the Schur complement is formed
    657: . mreuse   - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in `S`
    658: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming `Sp`:
    659:              `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
    660: - preuse   - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in `Sp`

    662:   Output Parameters:
    663: + S  - exact Schur complement, often of type `MATSCHURCOMPLEMENT` which is difficult to use for preconditioning
    664: - Sp - approximate Schur complement from which a preconditioner can be built $A11 - A10 inv(DIAGFORM(A00)) A01$

    666:   Level: advanced

    668:   Notes:
    669:   Since the real Schur complement is usually dense, providing a good approximation to `Sp` usually requires
    670:   application-specific information.

    672:   Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
    673:   and column index sets.  In that case, the user should call `PetscObjectComposeFunction()` on the *S matrix and pass mreuse of `MAT_REUSE_MATRIX` to set
    674:   "MatGetSchurComplement_C" to their function.  If their function needs to fall back to the default implementation, it
    675:   should call `MatGetSchurComplement_Basic()`.

    677:   `MatCreateSchurComplement()` takes as arguments the four submatrices and returns the virtual Schur complement (what this function returns in S).

    679:   `MatSchurComplementGetPmat()` takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).

    681:   In other words calling `MatCreateSchurComplement()` followed by `MatSchurComplementGetPmat()` produces the same output as this function but with slightly different
    682:   inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.

    684:   Developer Notes:
    685:   The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
    686:   remove redundancy and be clearer and simpler.

    688: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatCreateSchurComplement()`, `MatSchurComplementAinvType`
    689: @*/
    690: PetscErrorCode MatGetSchurComplement(Mat A, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
    691: {
    692:   PetscErrorCode (*f)(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatReuse, Mat *) = NULL;

    694:   PetscFunctionBegin;
    706:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
    707:   if (mreuse == MAT_REUSE_MATRIX) { /* This is the only situation, in which we can demand that the user pass a non-NULL pointer to non-garbage in S. */
    708:     PetscCall(PetscObjectQueryFunction((PetscObject)*S, "MatGetSchurComplement_C", &f));
    709:   }
    710:   if (f) PetscCall((*f)(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, preuse, Sp));
    711:   else PetscCall(MatGetSchurComplement_Basic(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, ainvtype, preuse, Sp));
    712:   PetscFunctionReturn(PETSC_SUCCESS);
    713: }

    715: /*@
    716:   MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming `Sp` in `MatSchurComplementGetPmat()`

    718:   Not Collective

    720:   Input Parameters:
    721: + S        - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of $A11 - A10 ksp(A00,Ap00) A01$
    722: - ainvtype - type of approximation to be used to form approximate Schur complement $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$:
    723:              `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`

    725:   Options Database Key:
    726: . -mat_schur_complement_ainv_type diag | lump | blockdiag | full - set schur complement type

    728:   Level: advanced

    730: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementGetAinvType()`
    731: @*/
    732: PetscErrorCode MatSchurComplementSetAinvType(Mat S, MatSchurComplementAinvType ainvtype)
    733: {
    734:   PetscBool            isschur;
    735:   Mat_SchurComplement *schur;

    737:   PetscFunctionBegin;
    739:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    740:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    742:   schur = (Mat_SchurComplement *)S->data;
    743:   PetscCheck(ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_FULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", (int)ainvtype);
    744:   schur->ainvtype = ainvtype;
    745:   PetscFunctionReturn(PETSC_SUCCESS);
    746: }

    748: /*@
    749:   MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming `Sp` in `MatSchurComplementGetPmat()`

    751:   Not Collective

    753:   Input Parameter:
    754: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of $A11 - A10 ksp(A00,Ap00) A01$

    756:   Output Parameter:
    757: . ainvtype - type of approximation used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
    758:              `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`

    760:   Level: advanced

    762: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementSetAinvType()`
    763: @*/
    764: PetscErrorCode MatSchurComplementGetAinvType(Mat S, MatSchurComplementAinvType *ainvtype)
    765: {
    766:   PetscBool            isschur;
    767:   Mat_SchurComplement *schur;

    769:   PetscFunctionBegin;
    771:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    772:   PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
    773:   schur = (Mat_SchurComplement *)S->data;
    774:   if (ainvtype) *ainvtype = schur->ainvtype;
    775:   PetscFunctionReturn(PETSC_SUCCESS);
    776: }

    778: /*@
    779:   MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by explicitly assembling the sparse matrix
    780:   $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$

    782:   Collective

    784:   Input Parameters:
    785: + A00      - the upper-left part of the original matrix $A = [A00 A01; A10 A11]$
    786: . A01      - (optional) the upper-right part of the original matrix $A = [A00 A01; A10 A11]$
    787: . A10      - (optional) the lower-left part of the original matrix $A = [A00 A01; A10 A11]$
    788: . A11      - (optional) the lower-right part of the original matrix $A = [A00 A01; A10 A11]$
    789: . ainvtype - type of approximation for DIAGFORM(A00) used when forming $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$. See `MatSchurComplementAinvType`.
    790: - preuse   - `MAT_INITIAL_MATRIX` for a new `Sp`, or `MAT_REUSE_MATRIX` to reuse an existing `Sp`, or `MAT_IGNORE_MATRIX` to put nothing in `Sp`

    792:   Output Parameter:
    793: . Sp - approximate Schur complement suitable for constructing a preconditioner for the true Schur complement $S = A11 - A10 inv(A00) A01$

    795:   Level: advanced

    797: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementAinvType`
    798: @*/
    799: PetscErrorCode MatCreateSchurComplementPmat(Mat A00, Mat A01, Mat A10, Mat A11, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
    800: {
    801:   PetscInt N00;

    803:   PetscFunctionBegin;
    804:   /* Use an appropriate approximate inverse of A00 to form A11 - A10 inv(DIAGFORM(A00)) A01; a NULL A01, A10 or A11 indicates a zero matrix. */
    805:   /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
    807:   if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);

    809:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
    810:   PetscCall(MatGetSize(A00, &N00, NULL));
    811:   if (!A01 || !A10 || !N00) {
    812:     if (preuse == MAT_INITIAL_MATRIX) {
    813:       PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, Sp));
    814:     } else { /* MAT_REUSE_MATRIX */
    815:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
    816:       PetscCall(MatCopy(A11, *Sp, DIFFERENT_NONZERO_PATTERN));
    817:     }
    818:   } else {
    819:     Mat       AdB, T;
    820:     Vec       diag;
    821:     PetscBool flg;

    823:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
    824:       PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
    825:       if (flg) {
    826:         PetscCall(MatTransposeGetMat(A01, &T));
    827:         PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &AdB));
    828:       } else {
    829:         PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
    830:         if (flg) {
    831:           PetscCall(MatHermitianTransposeGetMat(A01, &T));
    832:           PetscCall(MatHermitianTranspose(T, MAT_INITIAL_MATRIX, &AdB));
    833:         }
    834:       }
    835:       if (!flg) PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &AdB));
    836:       else {
    837:         PetscScalar shift, scale;

    839:         PetscCall(MatShellGetScalingShifts(A01, &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));
    840:         PetscCall(MatShift(AdB, shift));
    841:         PetscCall(MatScale(AdB, scale));
    842:       }
    843:       PetscCall(MatCreateVecs(A00, &diag, NULL));
    844:       if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
    845:         PetscCall(MatGetRowSum(A00, diag));
    846:       } else {
    847:         PetscCall(MatGetDiagonal(A00, diag));
    848:       }
    849:       PetscCall(VecReciprocal(diag));
    850:       PetscCall(MatDiagonalScale(AdB, diag, NULL));
    851:       PetscCall(VecDestroy(&diag));
    852:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
    853:       Mat      A00_inv;
    854:       MatType  type;
    855:       MPI_Comm comm;

    857:       PetscCall(PetscObjectGetComm((PetscObject)A00, &comm));
    858:       PetscCall(MatGetType(A00, &type));
    859:       PetscCall(MatCreate(comm, &A00_inv));
    860:       PetscCall(MatSetType(A00_inv, type));
    861:       PetscCall(MatInvertBlockDiagonalMat(A00, A00_inv));
    862:       PetscCall(MatMatMult(A00_inv, A01, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AdB));
    863:       PetscCall(MatDestroy(&A00_inv));
    864:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", ainvtype);
    865:     /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
    866:          MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult()  */
    867:     if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
    868:     PetscCall(MatMatMult(A10, AdB, MAT_INITIAL_MATRIX, PETSC_DETERMINE, Sp));
    869:     if (!A11) {
    870:       PetscCall(MatScale(*Sp, -1.0));
    871:     } else {
    872:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
    873:       PetscCall(MatAYPX(*Sp, -1, A11, DIFFERENT_NONZERO_PATTERN));
    874:     }
    875:     PetscCall(MatDestroy(&AdB));
    876:   }
    877:   PetscFunctionReturn(PETSC_SUCCESS);
    878: }

    880: static PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S, MatReuse preuse, Mat *Sp)
    881: {
    882:   Mat                  A, B, C, D;
    883:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;

    885:   PetscFunctionBegin;
    886:   if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
    887:   PetscCall(MatSchurComplementGetSubMatrices(S, &A, NULL, &B, &C, &D));
    888:   PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Schur complement component matrices unset");
    889:   if (schur->ainvtype != MAT_SCHUR_COMPLEMENT_AINV_FULL) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, schur->ainvtype, preuse, Sp));
    890:   else {
    891:     if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
    892:     PetscCall(MatSchurComplementComputeExplicitOperator(S, Sp));
    893:   }
    894:   PetscFunctionReturn(PETSC_SUCCESS);
    895: }

    897: /*@
    898:   MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$

    900:   Collective

    902:   Input Parameters:
    903: + S      - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of $A11 - A10 ksp(A00,Ap00) A01$
    904: - preuse - `MAT_INITIAL_MATRIX` for a new `Sp`, or `MAT_REUSE_MATRIX` to reuse an existing `Sp`, or `MAT_IGNORE_MATRIX` to put nothing in `Sp`

    906:   Output Parameter:
    907: . Sp - approximate Schur complement suitable for preconditioning the exact Schur complement $S = A11 - A10 inv(A00) A01$

    909:   Level: advanced

    911:   Notes:
    912:   The approximation of `Sp` depends on the argument passed to `MatSchurComplementSetAinvType()`
    913:   `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
    914:   -mat_schur_complement_ainv_type <diag,lump,blockdiag,full>

    916:   Sometimes users would like to provide problem-specific data in the Schur complement, usually only
    917:   for special row and column index sets.  In that case, the user should call `PetscObjectComposeFunction()` to set
    918:   "MatSchurComplementGetPmat_C" to their function.  If their function needs to fall back to the default implementation,
    919:   it should call `MatSchurComplementGetPmat_Basic()`.

    921:   Developer Notes:
    922:   The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
    923:   remove redundancy and be clearer and simpler.

    925:   This routine should be called `MatSchurComplementCreatePmat()`

    927: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetAinvType()`
    928: @*/
    929: PetscErrorCode MatSchurComplementGetPmat(Mat S, MatReuse preuse, Mat *Sp)
    930: {
    931:   PetscErrorCode (*f)(Mat, MatReuse, Mat *);

    933:   PetscFunctionBegin;
    937:   if (preuse != MAT_IGNORE_MATRIX) {
    938:     PetscAssertPointer(Sp, 3);
    939:     if (preuse == MAT_INITIAL_MATRIX) *Sp = NULL;
    941:   }
    942:   PetscCheck(!S->factortype, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

    944:   PetscCall(PetscObjectQueryFunction((PetscObject)S, "MatSchurComplementGetPmat_C", &f));
    945:   if (f) PetscCall((*f)(S, preuse, Sp));
    946:   else PetscCall(MatSchurComplementGetPmat_Basic(S, preuse, Sp));
    947:   PetscFunctionReturn(PETSC_SUCCESS);
    948: }

    950: static PetscErrorCode MatProductNumeric_SchurComplement_Dense(Mat C)
    951: {
    952:   Mat_Product         *product = C->product;
    953:   Mat_SchurComplement *Na      = (Mat_SchurComplement *)product->A->data;
    954:   Mat                  work1, work2;
    955:   PetscScalar         *v;
    956:   PetscInt             lda;

    958:   PetscFunctionBegin;
    959:   PetscCall(MatMatMult(Na->B, product->B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &work1));
    960:   PetscCall(MatDuplicate(work1, MAT_DO_NOT_COPY_VALUES, &work2));
    961:   PetscCall(KSPMatSolve(Na->ksp, work1, work2));
    962:   PetscCall(MatDestroy(&work1));
    963:   PetscCall(MatDenseGetArrayWrite(C, &v));
    964:   PetscCall(MatDenseGetLDA(C, &lda));
    965:   PetscCall(MatCreateDense(PetscObjectComm((PetscObject)C), C->rmap->n, C->cmap->n, C->rmap->N, C->cmap->N, v, &work1));
    966:   PetscCall(MatDenseSetLDA(work1, lda));
    967:   PetscCall(MatMatMult(Na->C, work2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &work1));
    968:   PetscCall(MatDenseRestoreArrayWrite(C, &v));
    969:   PetscCall(MatDestroy(&work2));
    970:   PetscCall(MatDestroy(&work1));
    971:   if (Na->D) {
    972:     PetscCall(MatMatMult(Na->D, product->B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &work1));
    973:     PetscCall(MatAYPX(C, -1.0, work1, SAME_NONZERO_PATTERN));
    974:     PetscCall(MatDestroy(&work1));
    975:   } else PetscCall(MatScale(C, -1.0));
    976:   PetscFunctionReturn(PETSC_SUCCESS);
    977: }

    979: static PetscErrorCode MatProductSymbolic_SchurComplement_Dense(Mat C)
    980: {
    981:   Mat_Product *product = C->product;
    982:   Mat          A = product->A, B = product->B;
    983:   PetscInt     m = A->rmap->n, n = B->cmap->n, M = A->rmap->N, N = B->cmap->N;
    984:   PetscBool    flg;

    986:   PetscFunctionBegin;
    987:   PetscCall(MatSetSizes(C, m, n, M, N));
    988:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &flg, MATSEQDENSE, MATMPIDENSE, ""));
    989:   if (!flg) {
    990:     PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
    991:     C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
    992:   }
    993:   PetscCall(MatSetUp(C));
    994:   C->ops->productnumeric = MatProductNumeric_SchurComplement_Dense;
    995:   PetscFunctionReturn(PETSC_SUCCESS);
    996: }

    998: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Dense(Mat C)
    999: {
    1000:   Mat_Product *product = C->product;

    1002:   PetscFunctionBegin;
    1003:   PetscCheck(product->type == MATPRODUCT_AB, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Not for product type %s", MatProductTypes[product->type]);
    1004:   C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
    1005:   PetscFunctionReturn(PETSC_SUCCESS);
    1006: }

    1008: /*MC
    1009:   MATSCHURCOMPLEMENT -  "schurcomplement" - Matrix type that behaves like the Schur complement of a matrix.

    1011:   Level: intermediate

    1013: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatType`, `MatCreateSchurComplement()`, `MatSchurComplementComputeExplicitOperator()`,
    1014:           `MatSchurComplementGetSubMatrices()`, `MatSchurComplementGetKSP()`
    1015: M*/
    1016: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
    1017: {
    1018:   Mat_SchurComplement *Na;

    1020:   PetscFunctionBegin;
    1021:   PetscCall(PetscNew(&Na));
    1022:   N->data = (void *)Na;

    1024:   N->ops->destroy        = MatDestroy_SchurComplement;
    1025:   N->ops->getvecs        = MatCreateVecs_SchurComplement;
    1026:   N->ops->view           = MatView_SchurComplement;
    1027:   N->ops->mult           = MatMult_SchurComplement;
    1028:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
    1029:   N->ops->multadd        = MatMultAdd_SchurComplement;
    1030:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
    1031:   N->assembled           = PETSC_FALSE;
    1032:   N->preallocated        = PETSC_FALSE;

    1034:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)N), &Na->ksp));
    1035:   PetscCall(PetscObjectChangeTypeName((PetscObject)N, MATSCHURCOMPLEMENT));
    1036:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", MatProductSetFromOptions_SchurComplement_Dense));
    1037:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", MatProductSetFromOptions_SchurComplement_Dense));
    1038:   PetscFunctionReturn(PETSC_SUCCESS);
    1039: }