Actual source code: ex9.c


  2: static char help[] = "Tests repeated setups and solves of PCFIELDSPLIT.\n\n";
  3: #include <petscksp.h>

  5: static PetscErrorCode replace_submats(Mat A)
  6: {
  7:   IS      *r, *c;
  8:   PetscInt i, j, nr, nc;

 10:   PetscFunctionBeginUser;
 11:   PetscCall(MatNestGetSubMats(A, &nr, &nc, NULL));
 12:   PetscCall(PetscMalloc1(nr, &r));
 13:   PetscCall(PetscMalloc1(nc, &c));
 14:   PetscCall(MatNestGetISs(A, r, c));
 15:   for (i = 0; i < nr; i++) {
 16:     for (j = 0; j < nc; j++) {
 17:       Mat         sA, nA;
 18:       const char *prefix;

 20:       PetscCall(MatCreateSubMatrix(A, r[i], c[j], MAT_INITIAL_MATRIX, &sA));
 21:       PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &nA));
 22:       PetscCall(MatGetOptionsPrefix(sA, &prefix));
 23:       PetscCall(MatSetOptionsPrefix(nA, prefix));
 24:       PetscCall(MatNestSetSubMat(A, i, j, nA));
 25:       PetscCall(MatDestroy(&nA));
 26:       PetscCall(MatDestroy(&sA));
 27:     }
 28:   }
 29:   PetscCall(PetscFree(r));
 30:   PetscCall(PetscFree(c));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: int main(int argc, char *argv[])
 35: {
 36:   KSP       ksp;
 37:   PC        pc;
 38:   Mat       M, A, P, sA[2][2], sP[2][2];
 39:   Vec       x, b;
 40:   IS        f[2];
 41:   PetscInt  i, j, rstart, rend;
 42:   PetscBool missA, missM;

 44:   PetscFunctionBeginUser;
 45:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 46:   PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &M));
 47:   PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
 48:   PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
 49:   PetscCall(MatShift(M, 1.));
 50:   PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
 51:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)M), 7, rstart, 1, &f[0]));
 52:   PetscCall(ISComplement(f[0], rstart, rend, &f[1]));
 53:   for (i = 0; i < 2; i++) {
 54:     for (j = 0; j < 2; j++) {
 55:       PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sA[i][j]));
 56:       PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sP[i][j]));
 57:     }
 58:   }
 59:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sA[0][0], &A));
 60:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sP[0][0], &P));

 62:   /* Tests MatMissingDiagonal_Nest */
 63:   PetscCall(MatMissingDiagonal(M, &missM, NULL));
 64:   PetscCall(MatMissingDiagonal(A, &missA, NULL));
 65:   if (missM != missA) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Unexpected %s != %s\n", missM ? "true" : "false", missA ? "true" : "false"));

 67:   PetscCall(MatDestroy(&M));

 69:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
 70:   PetscCall(KSPSetOperators(ksp, A, P));
 71:   PetscCall(KSPGetPC(ksp, &pc));
 72:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
 73:   PetscCall(KSPSetFromOptions(ksp));
 74:   PetscCall(MatCreateVecs(A, &x, &b));
 75:   PetscCall(VecSetRandom(b, NULL));
 76:   PetscCall(KSPSolve(ksp, b, x));
 77:   PetscCall(replace_submats(A));
 78:   PetscCall(replace_submats(P));
 79:   PetscCall(KSPSolve(ksp, b, x));

 81:   PetscCall(KSPDestroy(&ksp));
 82:   PetscCall(VecDestroy(&x));
 83:   PetscCall(VecDestroy(&b));
 84:   PetscCall(MatDestroy(&A));
 85:   PetscCall(MatDestroy(&P));
 86:   for (i = 0; i < 2; i++) {
 87:     PetscCall(ISDestroy(&f[i]));
 88:     for (j = 0; j < 2; j++) {
 89:       PetscCall(MatDestroy(&sA[i][j]));
 90:       PetscCall(MatDestroy(&sP[i][j]));
 91:     }
 92:   }
 93:   PetscCall(PetscFinalize());
 94:   return 0;
 95: }

 97: /*TEST

 99:    test:
100:      nsize: 1
101:      filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
102:      args: -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_type {{additive multiplicative}} -ksp_converged_reason -ksp_error_if_not_converged

104:    test:
105:      suffix: schur
106:      nsize: 1
107:      filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
108:      args: -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_type schur -pc_fieldsplit_schur_scale 1.0 -pc_fieldsplit_schur_fact_type {{diag lower upper full}} -ksp_converged_reason -ksp_error_if_not_converged

110: TEST*/