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:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
    135:   PetscFunctionReturn(PETSC_SUCCESS);
    136: }

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

    141:   Collective

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

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

    153:   Level: intermediate

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

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

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

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

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

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

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

    188:   Collective

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

    198:   Level: intermediate

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

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

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

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

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

    222:   PetscFunctionBegin;
    223:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    224:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    225:   PetscCheck(!S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementUpdateSubMatrices() for already used matrix");
    230:   PetscCheckSameComm(A00, 2, Ap00, 3);
    231:   PetscCheckSameComm(A00, 2, A01, 4);
    232:   PetscCheckSameComm(A00, 2, A10, 5);
    233:   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);
    234:   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);
    235:   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);
    236:   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);
    237:   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);
    238:   if (A11) {
    240:     PetscCheckSameComm(A00, 2, A11, 6);
    241:     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);
    242:   }

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

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

    264:   Not Collective

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

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

    272:   Level: intermediate

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

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

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

    294:   Not Collective

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

    300:   Level: developer

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

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

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

    326: /*@
    327:   MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

    329:   Collective

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

    339:   Level: intermediate

    341:   Notes:
    342:   All four matrices must have the same MPI communicator

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

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

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

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

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

    361:   PetscFunctionBegin;
    363:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    364:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    365:   PetscCheck(S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementSetSubMatrices() for a new matrix");
    370:   PetscCheckSameComm(A00, 2, Ap00, 3);
    371:   PetscCheckSameComm(A00, 2, A01, 4);
    372:   PetscCheckSameComm(A00, 2, A10, 5);
    373:   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);
    374:   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);
    375:   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);
    376:   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);
    377:   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);
    378:   if (A11) {
    380:     PetscCheckSameComm(A00, 2, A11, 6);
    381:     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);
    382:   }

    384:   PetscCall(PetscObjectReference((PetscObject)A00));
    385:   PetscCall(PetscObjectReference((PetscObject)Ap00));
    386:   PetscCall(PetscObjectReference((PetscObject)A01));
    387:   PetscCall(PetscObjectReference((PetscObject)A10));
    388:   PetscCall(PetscObjectReference((PetscObject)A11));

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

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

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

    406: /*@
    407:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

    409:   Collective

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

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

    421:   Level: intermediate

    423:   Note:
    424:   Use `NULL` for any unneeded output argument.

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

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

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

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

    449: /*@
    450:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

    452:   Collective

    454:   Input Parameter:
    455: . A - the matrix obtained with `MatCreateSchurComplement()`

    457:   Output Parameter:
    458: . S - the Schur complement matrix

    460:   Level: advanced

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

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

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

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

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

    486:   PetscFunctionBegin;
    487:   PetscAssertPointer(S, 2);
    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, PETSC_DECIDE, 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:     PetscMemType type;

    525:     PetscCall(MatDenseGetArrayWriteAndMemType(AinvBd, &v, &type)); // no easy way to resize a Mat, so create a new one with the same data pointer
    526:     PetscCall(MatCreateDenseWithMemType(PetscObjectComm((PetscObject)A), type, AinvBd->rmap->n, A->cmap->n, AinvBd->rmap->N, A->cmap->N, PETSC_DECIDE, v, &AinvB));
    527:     PetscCall(MatDenseReplaceArrayWithMemType(AinvB, type, v)); // let MatDestroy() free the data pointer
    528:     PetscCall(MatDenseRestoreArrayWriteAndMemType(AinvBd, &v));
    529:     PetscCall(MatHeaderReplace(AinvBd, &AinvB)); // replace the input composed Mat with just A00^-1 A01 (trailing columns are removed)
    530:   }
    531:   PetscCall(MatDestroy(&Bd));
    532:   if (!set) PetscCall(MatFilter(AinvBd, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
    533:   if (D && !*S) {
    534:     PetscCall(MatGetLocalSize(D, &m, &n));
    535:     PetscCall(MatGetSize(D, &M, &N));
    536:     PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)A), vtype, m, n, M, N, PETSC_DECIDE, NULL, S));
    537:   } else if (*S) {
    538:     PetscCall(MatGetType(AinvBd, &mtype));
    539:     PetscCall(MatSetType(*S, mtype));
    540:   }
    541:   PetscCall(MatMatMult(C, AinvBd, *S ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DETERMINE, S));
    542:   if (!set) PetscCall(MatDestroy(&AinvBd));
    543:   else {
    544:     PetscCall(MatScale(AinvBd, -1.0));
    545:     PetscCall(MatFilter(AinvBd, PETSC_MACHINE_EPSILON, PETSC_FALSE, PETSC_FALSE));
    546:     PetscCall(MatFilter(*S, PETSC_MACHINE_EPSILON, PETSC_FALSE, PETSC_FALSE));
    547:   }
    548:   if (D) {
    549:     PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
    550:     if (flg) {
    551:       PetscCall(MatIsSymmetricKnown(A, &set, &symm));
    552:       if (!set || !symm) PetscCall(MatConvert(D, MATBAIJ, MAT_INITIAL_MATRIX, &E)); /* convert the (1,1) block to nonsymmetric storage for MatAXPY() */
    553:     }
    554:     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 */
    555:     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 */
    556:   }
    557:   PetscCall(MatDestroy(&E));
    558:   PetscCall(MatScale(*S, -1.0));
    559:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sexplicit_operator_", ((PetscObject)A)->prefix ? ((PetscObject)A)->prefix : ""));
    560:   PetscCall(MatSetOptionsPrefix(*S, prefix));
    561:   PetscObjectOptionsBegin((PetscObject)*S);
    562:   PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, !E && flg ? MATSBAIJ : mtype, type, 256, &set));
    563:   if (set) PetscCall(MatConvert(*S, type, MAT_INPLACE_MATRIX, S));
    564:   flg = PETSC_FALSE;
    565:   PetscCall(PetscOptionsBool("-mat_symmetric", "Sets the MAT_SYMMETRIC option", "MatSetOption", flg, &flg, &set));
    566:   if (set) PetscCall(MatSetOption(*S, MAT_SYMMETRIC, flg));
    567:   if (PetscDefined(USE_COMPLEX)) {
    568:     flg = PETSC_FALSE;
    569:     PetscCall(PetscOptionsBool("-mat_hermitian", "Sets the MAT_HERMITIAN option", "MatSetOption", flg, &flg, &set));
    570:     if (set) PetscCall(MatSetOption(*S, MAT_HERMITIAN, flg));
    571:   }
    572:   PetscOptionsEnd();
    573:   PetscFunctionReturn(PETSC_SUCCESS);
    574: }

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

    583:   PetscFunctionBegin;
    593:   if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);

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

    599:   reuse = MAT_INITIAL_MATRIX;
    600:   if (mreuse == MAT_REUSE_MATRIX) {
    601:     PetscCall(MatSchurComplementGetSubMatrices(*S, &A, &Ap, &B, &C, &D));
    602:     PetscCheck(A && Ap && B && C, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Attempting to reuse matrix but Schur complement matrices unset");
    603:     PetscCheck(A == Ap, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Matrix for constructing the preconditioner does not match operator");
    604:     PetscCall(MatDestroy(&Ap)); /* get rid of extra reference */
    605:     reuse = MAT_REUSE_MATRIX;
    606:   }
    607:   PetscCall(MatCreateSubMatrix(mat, isrow0, iscol0, reuse, &A));
    608:   PetscCall(MatCreateSubMatrix(mat, isrow0, iscol1, reuse, &B));
    609:   PetscCall(MatCreateSubMatrix(mat, isrow1, iscol0, reuse, &C));
    610:   PetscCall(MatCreateSubMatrix(mat, isrow1, iscol1, reuse, &D));
    611:   switch (mreuse) {
    612:   case MAT_INITIAL_MATRIX:
    613:     PetscCall(MatCreateSchurComplement(A, A, B, C, D, S));
    614:     break;
    615:   case MAT_REUSE_MATRIX:
    616:     PetscCall(MatSchurComplementUpdateSubMatrices(*S, A, A, B, C, D));
    617:     break;
    618:   default:
    619:     PetscCheck(mreuse == MAT_IGNORE_MATRIX, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unrecognized value of mreuse %d", (int)mreuse);
    620:   }
    621:   if (preuse != MAT_IGNORE_MATRIX) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, ainvtype, preuse, Sp));
    622:   PetscCall(MatDestroy(&A));
    623:   PetscCall(MatDestroy(&B));
    624:   PetscCall(MatDestroy(&C));
    625:   PetscCall(MatDestroy(&D));
    626:   PetscFunctionReturn(PETSC_SUCCESS);
    627: }

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

    632:   Collective

    634:   Input Parameters:
    635: + A        - matrix in which the complement is to be taken
    636: . isrow0   - rows to eliminate
    637: . iscol0   - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
    638: . isrow1   - rows in which the Schur complement is formed
    639: . iscol1   - columns in which the Schur complement is formed
    640: . mreuse   - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in `S`
    641: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming `Sp`:
    642:              `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
    643: - preuse   - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in `Sp`

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

    649:   Level: advanced

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

    655:   Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
    656:   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
    657:   "MatGetSchurComplement_C" to their function.  If their function needs to fall back to the default implementation, it
    658:   should call `MatGetSchurComplement_Basic()`.

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

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

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

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

    671: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatCreateSchurComplement()`, `MatSchurComplementAinvType`
    672: @*/
    673: PetscErrorCode MatGetSchurComplement(Mat A, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
    674: {
    675:   PetscErrorCode (*f)(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatReuse, Mat *) = NULL;

    677:   PetscFunctionBegin;
    689:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
    690:   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. */
    691:     PetscCall(PetscObjectQueryFunction((PetscObject)*S, "MatGetSchurComplement_C", &f));
    692:   }
    693:   if (f) PetscCall((*f)(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, preuse, Sp));
    694:   else PetscCall(MatGetSchurComplement_Basic(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, ainvtype, preuse, Sp));
    695:   PetscFunctionReturn(PETSC_SUCCESS);
    696: }

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

    701:   Not Collective

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

    708:   Options Database Key:
    709: . -mat_schur_complement_ainv_type (diag|lump|blockdiag|full) - set Schur complement type

    711:   Level: advanced

    713: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementGetAinvType()`
    714: @*/
    715: PetscErrorCode MatSchurComplementSetAinvType(Mat S, MatSchurComplementAinvType ainvtype)
    716: {
    717:   PetscBool            isschur;
    718:   Mat_SchurComplement *schur;

    720:   PetscFunctionBegin;
    722:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    723:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
    725:   schur = (Mat_SchurComplement *)S->data;
    726:   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);
    727:   schur->ainvtype = ainvtype;
    728:   PetscFunctionReturn(PETSC_SUCCESS);
    729: }

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

    734:   Not Collective

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

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

    743:   Level: advanced

    745: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementSetAinvType()`
    746: @*/
    747: PetscErrorCode MatSchurComplementGetAinvType(Mat S, MatSchurComplementAinvType *ainvtype)
    748: {
    749:   PetscBool            isschur;
    750:   Mat_SchurComplement *schur;

    752:   PetscFunctionBegin;
    754:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    755:   PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
    756:   schur = (Mat_SchurComplement *)S->data;
    757:   if (ainvtype) *ainvtype = schur->ainvtype;
    758:   PetscFunctionReturn(PETSC_SUCCESS);
    759: }

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

    765:   Collective

    767:   Input Parameters:
    768: + A00      - the upper-left part of the original matrix $A = [A00 A01; A10 A11]$
    769: . A01      - (optional) the upper-right part of the original matrix $A = [A00 A01; A10 A11]$
    770: . A10      - (optional) the lower-left part of the original matrix $A = [A00 A01; A10 A11]$
    771: . A11      - (optional) the lower-right part of the original matrix $A = [A00 A01; A10 A11]$
    772: . ainvtype - type of approximation for DIAGFORM(A00) used when forming $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$. See `MatSchurComplementAinvType`.
    773: - 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`

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

    778:   Level: advanced

    780: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementAinvType`
    781: @*/
    782: PetscErrorCode MatCreateSchurComplementPmat(Mat A00, Mat A01, Mat A10, Mat A11, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
    783: {
    784:   PetscInt N00;

    786:   PetscFunctionBegin;
    787:   /* 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. */
    788:   /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
    790:   if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);

    792:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
    793:   PetscCall(MatGetSize(A00, &N00, NULL));
    794:   if (!A01 || !A10 || !N00) {
    795:     if (preuse == MAT_INITIAL_MATRIX) {
    796:       PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, Sp));
    797:     } else { /* MAT_REUSE_MATRIX */
    798:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
    799:       PetscCall(MatCopy(A11, *Sp, DIFFERENT_NONZERO_PATTERN));
    800:     }
    801:   } else {
    802:     Mat       AdB, T;
    803:     Vec       diag;
    804:     PetscBool flg;

    806:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
    807:       PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
    808:       if (flg) {
    809:         PetscCall(MatTransposeGetMat(A01, &T));
    810:         PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &AdB));
    811:       } else {
    812:         PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
    813:         if (flg) {
    814:           PetscCall(MatHermitianTransposeGetMat(A01, &T));
    815:           PetscCall(MatHermitianTranspose(T, MAT_INITIAL_MATRIX, &AdB));
    816:         }
    817:       }
    818:       if (!flg) PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &AdB));
    819:       else {
    820:         PetscScalar shift, scale;

    822:         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));
    823:         PetscCall(MatShift(AdB, shift));
    824:         PetscCall(MatScale(AdB, scale));
    825:       }
    826:       PetscCall(MatCreateVecs(A00, &diag, NULL));
    827:       if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) PetscCall(MatGetRowSum(A00, diag));
    828:       else PetscCall(MatGetDiagonal(A00, diag));
    829:       PetscCall(VecReciprocal(diag));
    830:       PetscCall(MatDiagonalScale(AdB, diag, NULL));
    831:       PetscCall(VecDestroy(&diag));
    832:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
    833:       Mat      A00_inv;
    834:       MatType  type;
    835:       MPI_Comm comm;

    837:       PetscCall(PetscObjectGetComm((PetscObject)A00, &comm));
    838:       PetscCall(MatGetType(A00, &type));
    839:       PetscCall(MatCreate(comm, &A00_inv));
    840:       PetscCall(MatSetType(A00_inv, type));
    841:       PetscCall(MatInvertBlockDiagonalMat(A00, A00_inv));
    842:       PetscCall(MatMatMult(A00_inv, A01, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AdB));
    843:       PetscCall(MatDestroy(&A00_inv));
    844:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", ainvtype);
    845:     /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
    846:          MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult()  */
    847:     if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
    848:     PetscCall(MatMatMult(A10, AdB, MAT_INITIAL_MATRIX, PETSC_DETERMINE, Sp));
    849:     PetscCall(MatScale(*Sp, -1.0));
    850:     if (A11) { /* TODO: when can we pass SAME_NONZERO_PATTERN? */
    851:       PetscCall(MatAXPY(*Sp, 1.0, A11, DIFFERENT_NONZERO_PATTERN));
    852:     }
    853:     PetscCall(MatDestroy(&AdB));
    854:   }
    855:   PetscFunctionReturn(PETSC_SUCCESS);
    856: }

    858: static PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S, MatReuse preuse, Mat *Sp)
    859: {
    860:   Mat                  A, B, C, D;
    861:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
    862:   MatNullSpace         sp;

    864:   PetscFunctionBegin;
    865:   if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
    866:   PetscCall(MatSchurComplementGetSubMatrices(S, &A, NULL, &B, &C, &D));
    867:   PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Schur complement component matrices unset");
    868:   if (schur->ainvtype != MAT_SCHUR_COMPLEMENT_AINV_FULL) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, schur->ainvtype, preuse, Sp));
    869:   else {
    870:     if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
    871:     PetscCall(MatSchurComplementComputeExplicitOperator(S, Sp));
    872:   }
    873:   /* If the Schur complement has a nullspace, then Sp nullspace contains it, independently of the ainv type */
    874:   PetscCall(MatGetNullSpace(S, &sp));
    875:   if (sp) PetscCall(MatSetNullSpace(*Sp, sp));
    876:   PetscCall(MatGetTransposeNullSpace(S, &sp));
    877:   if (sp) PetscCall(MatSetTransposeNullSpace(*Sp, sp));
    878:   PetscFunctionReturn(PETSC_SUCCESS);
    879: }

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

    884:   Collective

    886:   Input Parameters:
    887: + S      - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of $A11 - A10 ksp(A00,Ap00) A01$
    888: - 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`

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

    893:   Level: advanced

    895:   Notes:
    896:   The approximation of `Sp` depends on the argument passed to `MatSchurComplementSetAinvType()`
    897:   `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
    898:   -mat_schur_complement_ainv_type <diag,lump,blockdiag,full>

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

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

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

    911: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetAinvType()`
    912: @*/
    913: PetscErrorCode MatSchurComplementGetPmat(Mat S, MatReuse preuse, Mat *Sp)
    914: {
    915:   PetscErrorCode (*f)(Mat, MatReuse, Mat *);

    917:   PetscFunctionBegin;
    921:   if (preuse != MAT_IGNORE_MATRIX) {
    922:     PetscAssertPointer(Sp, 3);
    923:     if (preuse == MAT_INITIAL_MATRIX) *Sp = NULL;
    925:   }
    926:   PetscCheck(!S->factortype, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

    928:   PetscCall(PetscObjectQueryFunction((PetscObject)S, "MatSchurComplementGetPmat_C", &f));
    929:   if (f) PetscCall((*f)(S, preuse, Sp));
    930:   else PetscCall(MatSchurComplementGetPmat_Basic(S, preuse, Sp));
    931:   PetscFunctionReturn(PETSC_SUCCESS);
    932: }

    934: static PetscErrorCode MatProductNumeric_SchurComplement_Dense(Mat C)
    935: {
    936:   Mat_Product         *product = C->product;
    937:   Mat_SchurComplement *Na      = (Mat_SchurComplement *)product->A->data;
    938:   Mat                  work1, work2;
    939:   PetscScalar         *v;
    940:   PetscInt             lda;

    942:   PetscFunctionBegin;
    943:   PetscCall(MatMatMult(Na->B, product->B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &work1));
    944:   PetscCall(MatDuplicate(work1, MAT_DO_NOT_COPY_VALUES, &work2));
    945:   PetscCall(KSPMatSolve(Na->ksp, work1, work2));
    946:   PetscCall(MatDestroy(&work1));
    947:   PetscCall(MatDenseGetArrayWrite(C, &v));
    948:   PetscCall(MatDenseGetLDA(C, &lda));
    949:   PetscCall(MatCreateDense(PetscObjectComm((PetscObject)C), C->rmap->n, C->cmap->n, C->rmap->N, C->cmap->N, v, &work1));
    950:   PetscCall(MatDenseSetLDA(work1, lda));
    951:   PetscCall(MatMatMult(Na->C, work2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &work1));
    952:   PetscCall(MatDenseRestoreArrayWrite(C, &v));
    953:   PetscCall(MatDestroy(&work2));
    954:   PetscCall(MatDestroy(&work1));
    955:   if (Na->D) {
    956:     PetscCall(MatMatMult(Na->D, product->B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &work1));
    957:     PetscCall(MatAYPX(C, -1.0, work1, SAME_NONZERO_PATTERN));
    958:     PetscCall(MatDestroy(&work1));
    959:   } else PetscCall(MatScale(C, -1.0));
    960:   PetscFunctionReturn(PETSC_SUCCESS);
    961: }

    963: static PetscErrorCode MatProductSymbolic_SchurComplement_Dense(Mat C)
    964: {
    965:   Mat_Product *product = C->product;
    966:   Mat          A = product->A, B = product->B;
    967:   PetscInt     m = A->rmap->n, n = B->cmap->n, M = A->rmap->N, N = B->cmap->N;
    968:   PetscBool    flg;

    970:   PetscFunctionBegin;
    971:   PetscCall(MatSetSizes(C, m, n, M, N));
    972:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &flg, MATSEQDENSE, MATMPIDENSE, ""));
    973:   if (!flg) {
    974:     PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
    975:     C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
    976:   }
    977:   PetscCall(MatSetUp(C));
    978:   C->ops->productnumeric = MatProductNumeric_SchurComplement_Dense;
    979:   PetscFunctionReturn(PETSC_SUCCESS);
    980: }

    982: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Dense(Mat C)
    983: {
    984:   Mat_Product *product = C->product;

    986:   PetscFunctionBegin;
    987:   if (product->type != MATPRODUCT_AB) PetscFunctionReturn(PETSC_SUCCESS);
    988:   C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
    989:   PetscFunctionReturn(PETSC_SUCCESS);
    990: }

    992: static PetscErrorCode MatProductNumeric_SchurComplement_Any(Mat C)
    993: {
    994:   Mat_SchurComplement *Na = (Mat_SchurComplement *)C->data;

    996:   PetscFunctionBegin;
    997:   if (Na->D && Na->D->product) PetscCall(MatProductNumeric(Na->D));
    998:   if (Na->B->product) PetscCall(MatProductNumeric(Na->B));
    999:   if (Na->C->product) PetscCall(MatProductNumeric(Na->C));
    1000:   C->assembled = PETSC_TRUE;
    1001:   PetscFunctionReturn(PETSC_SUCCESS);
    1002: }

    1004: static PetscErrorCode MatProductSymbolic_SchurComplement_Any(Mat C)
    1005: {
    1006:   Mat_SchurComplement *Na = (Mat_SchurComplement *)C->data;

    1008:   PetscFunctionBegin;
    1009:   if (Na->D && Na->D->product) PetscCall(MatProductSymbolic(Na->D));
    1010:   if (Na->B->product) PetscCall(MatProductSymbolic(Na->B));
    1011:   if (Na->C->product) PetscCall(MatProductSymbolic(Na->C));
    1012:   C->ops->productnumeric = MatProductNumeric_SchurComplement_Any;
    1013:   C->preallocated        = PETSC_TRUE;
    1014:   C->assembled           = PETSC_FALSE;
    1015:   PetscFunctionReturn(PETSC_SUCCESS);
    1016: }

    1018: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Any(Mat C)
    1019: {
    1020:   Mat_Product         *product = C->product;
    1021:   Mat_SchurComplement *Na, *Ca;
    1022:   Mat                  B = product->B, S = product->A, pB = NULL, pC = NULL, pD = NULL;
    1023:   KSP                  ksp;
    1024:   PetscInt             m = PETSC_DECIDE, n = PETSC_DECIDE, M = PETSC_DECIDE, N = PETSC_DECIDE;
    1025:   MatProductType       pbtype = MATPRODUCT_UNSPECIFIED, pctype = MATPRODUCT_UNSPECIFIED;
    1026:   PetscBool            isschur;

    1028:   PetscFunctionBegin;
    1029:   if (product->type == MATPRODUCT_ABC || product->type == MATPRODUCT_AtB) PetscFunctionReturn(PETSC_SUCCESS);
    1030:   /* A * S not yet supported (should be easy though) */
    1031:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
    1032:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);

    1034:   Na = (Mat_SchurComplement *)S->data;
    1035:   if (Na->D) {
    1036:     PetscCall(MatProductCreate(Na->D, B, NULL, &pD));
    1037:     PetscCall(MatProductSetType(pD, product->type));
    1038:     PetscCall(MatProductSetFromOptions(pD));
    1039:   }
    1040:   if (pD && !pD->ops->productsymbolic) {
    1041:     PetscCall(MatDestroy(&pD));
    1042:     PetscFunctionReturn(PETSC_SUCCESS);
    1043:   }

    1045:   /* S = A11 - A10 M A01 */
    1046:   switch (product->type) {
    1047:   case MATPRODUCT_AB: /* A11 B - A10 * M * A01 * B */
    1048:     pbtype = product->type;
    1049:     PetscCall(PetscObjectReference((PetscObject)Na->C));
    1050:     pC = Na->C;
    1051:     m  = S->rmap->n;
    1052:     M  = S->rmap->N;
    1053:     n  = B->cmap->n;
    1054:     N  = B->cmap->N;
    1055:     break;
    1056:   case MATPRODUCT_ABt: /* A11 B^t - A10 * M * A01 * B^t */
    1057:     pbtype = product->type;
    1058:     PetscCall(PetscObjectReference((PetscObject)Na->C));
    1059:     pC = Na->C;
    1060:     m  = S->rmap->n;
    1061:     M  = S->rmap->N;
    1062:     n  = B->rmap->n;
    1063:     N  = B->rmap->N;
    1064:     break;
    1065:   case MATPRODUCT_PtAP: /* Pt A11 P - Pt * A10 * M * A01 * P */
    1066:     pbtype = MATPRODUCT_AB;
    1067:     pctype = MATPRODUCT_AtB;
    1068:     m      = B->cmap->n;
    1069:     M      = B->cmap->N;
    1070:     n      = B->cmap->n;
    1071:     N      = B->cmap->N;
    1072:     break;
    1073:   case MATPRODUCT_RARt: /* R A11 Rt - R * A10 * M * A01 * Rt */
    1074:     pbtype = MATPRODUCT_ABt;
    1075:     pctype = MATPRODUCT_AB;
    1076:     m      = B->rmap->n;
    1077:     M      = B->rmap->N;
    1078:     n      = B->rmap->n;
    1079:     N      = B->rmap->N;
    1080:     break;
    1081:   default:
    1082:     break;
    1083:   }
    1084:   PetscCall(MatProductCreate(Na->B, B, NULL, &pB));
    1085:   PetscCall(MatProductSetType(pB, pbtype));
    1086:   PetscCall(MatProductSetFromOptions(pB));
    1087:   if (!pB->ops->productsymbolic) {
    1088:     PetscCall(MatDestroy(&pB));
    1089:     PetscFunctionReturn(PETSC_SUCCESS);
    1090:   }
    1091:   if (pC == NULL) { /* Some work can in principle be saved here if we recognize symmetry */
    1092:     PetscCall(MatProductCreate(B, Na->C, NULL, &pC));
    1093:     PetscCall(MatProductSetType(pC, pctype));
    1094:     PetscCall(MatProductSetFromOptions(pC));
    1095:     if (!pC->ops->productsymbolic) {
    1096:       PetscCall(MatDestroy(&pC));
    1097:       PetscFunctionReturn(PETSC_SUCCESS);
    1098:     }
    1099:   }
    1100:   PetscCall(MatSetType(C, MATSCHURCOMPLEMENT));
    1101:   PetscCall(MatSetSizes(C, m, n, M, N));
    1102:   PetscCall(PetscLayoutSetUp(C->rmap));
    1103:   PetscCall(PetscLayoutSetUp(C->cmap));
    1104:   PetscCall(PetscObjectReference((PetscObject)Na->A));
    1105:   PetscCall(PetscObjectReference((PetscObject)Na->Ap));
    1106:   Ca                      = (Mat_SchurComplement *)C->data;
    1107:   Ca->A                   = Na->A;
    1108:   Ca->Ap                  = Na->Ap;
    1109:   Ca->B                   = pB;
    1110:   Ca->C                   = pC;
    1111:   Ca->D                   = pD;
    1112:   C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Any;
    1113:   PetscCall(MatSchurComplementGetKSP(S, &ksp));
    1114:   PetscCall(MatSchurComplementSetKSP(C, ksp));
    1115:   PetscFunctionReturn(PETSC_SUCCESS);
    1116: }

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

    1121:   Level: intermediate

    1123: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatType`, `MatCreateSchurComplement()`, `MatSchurComplementComputeExplicitOperator()`,
    1124:           `MatSchurComplementGetSubMatrices()`, `MatSchurComplementGetKSP()`
    1125: M*/
    1126: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
    1127: {
    1128:   Mat_SchurComplement *Na;

    1130:   PetscFunctionBegin;
    1131:   PetscCall(PetscNew(&Na));
    1132:   N->data = (void *)Na;

    1134:   N->ops->destroy        = MatDestroy_SchurComplement;
    1135:   N->ops->getvecs        = MatCreateVecs_SchurComplement;
    1136:   N->ops->view           = MatView_SchurComplement;
    1137:   N->ops->mult           = MatMult_SchurComplement;
    1138:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
    1139:   N->ops->multadd        = MatMultAdd_SchurComplement;
    1140:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
    1141:   N->assembled           = PETSC_FALSE;
    1142:   N->preallocated        = PETSC_FALSE;

    1144:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)N), &Na->ksp));
    1145:   PetscCall(PetscObjectChangeTypeName((PetscObject)N, MATSCHURCOMPLEMENT));
    1146:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", MatProductSetFromOptions_SchurComplement_Dense));
    1147:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", MatProductSetFromOptions_SchurComplement_Dense));
    1148:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_SchurComplement_Any));
    1149:   PetscFunctionReturn(PETSC_SUCCESS);
    1150: }