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;

 11:   MatNestGetSubMats(A, &nr, &nc, NULL);
 12:   PetscMalloc1(nr, &r);
 13:   PetscMalloc1(nc, &c);
 14:   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:       MatCreateSubMatrix(A, r[i], c[j], MAT_INITIAL_MATRIX, &sA);
 21:       MatDuplicate(sA, MAT_COPY_VALUES, &nA);
 22:       MatGetOptionsPrefix(sA, &prefix);
 23:       MatSetOptionsPrefix(nA, prefix);
 24:       MatNestSetSubMat(A, i, j, nA);
 25:       MatDestroy(&nA);
 26:       MatDestroy(&sA);
 27:     }
 28:   }
 29:   PetscFree(r);
 30:   PetscFree(c);
 31:   return 0;
 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;

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

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

 67:   MatDestroy(&M);

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

 81:   KSPDestroy(&ksp);
 82:   VecDestroy(&x);
 83:   VecDestroy(&b);
 84:   MatDestroy(&A);
 85:   MatDestroy(&P);
 86:   for (i = 0; i < 2; i++) {
 87:     ISDestroy(&f[i]);
 88:     for (j = 0; j < 2; j++) {
 89:       MatDestroy(&sA[i][j]);
 90:       MatDestroy(&sP[i][j]);
 91:     }
 92:   }
 93:   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*/