Actual source code: ex8.c

  1: static char help[] = "Solves a linear system in parallel with KSP. \n\
  2: Contributed by Jose E. Roman, SLEPc developer, for testing repeated call of KSPSetOperators(), 2014 \n\n";

  4: #include <petscksp.h>
  5: int main(int argc, char **args)
  6: {
  7:   Vec           x, b, u; /* approx solution, RHS, exact solution */
  8:   Mat           A;       /* linear system matrix */
  9:   KSP           ksp;     /* linear solver context */
 10:   PetscRandom   rctx;    /* random number generator context */
 11:   PetscInt      i, j, Ii, J, Istart, Iend, m = 8, n = 7;
 12:   PetscBool     flg = PETSC_FALSE;
 13:   PetscScalar   v;
 14:   PC            pc;
 15:   PetscInt      in;
 16:   Mat           F, B;
 17:   PetscBool     solve = PETSC_FALSE, sameA = PETSC_FALSE, setfromoptions_first = PETSC_FALSE;
 18:   PetscLogStage stage;
 19: #if !defined(PETSC_HAVE_MUMPS)
 20:   PetscMPIInt size;
 21: #endif

 23:   PetscFunctionBeginUser;
 24:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 27:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 28:          Compute the matrix and right-hand-side vector that define
 29:          the linear system, Ax = b.
 30:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 31:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 32:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
 33:   PetscCall(MatSetFromOptions(A));
 34:   PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
 35:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
 36:   PetscCall(MatSetUp(A));

 38:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

 40:   PetscCall(PetscLogStageRegister("Assembly", &stage));
 41:   PetscCall(PetscLogStagePush(stage));
 42:   for (Ii = Istart; Ii < Iend; Ii++) {
 43:     v = -1.0;
 44:     i = Ii / n;
 45:     j = Ii - i * n;
 46:     if (i > 0) {
 47:       J = Ii - n;
 48:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 49:     }
 50:     if (i < m - 1) {
 51:       J = Ii + n;
 52:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 53:     }
 54:     if (j > 0) {
 55:       J = Ii - 1;
 56:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 57:     }
 58:     if (j < n - 1) {
 59:       J = Ii + 1;
 60:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 61:     }
 62:     v = 4.0;
 63:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
 64:   }
 65:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 66:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 67:   PetscCall(PetscLogStagePop());

 69:   /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
 70:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));

 72:   /* Create parallel vectors. */
 73:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 74:   PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n));
 75:   PetscCall(VecSetFromOptions(u));
 76:   PetscCall(VecDuplicate(u, &b));
 77:   PetscCall(VecDuplicate(b, &x));

 79:   /*
 80:      Set exact solution; then compute right-hand-side vector.
 81:      By default we use an exact solution of a vector with all
 82:      elements of 1.0;  Alternatively, using the runtime option
 83:      -random_sol forms a solution vector with random components.
 84:   */
 85:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL));
 86:   if (flg) {
 87:     PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
 88:     PetscCall(PetscRandomSetFromOptions(rctx));
 89:     PetscCall(VecSetRandom(u, rctx));
 90:     PetscCall(PetscRandomDestroy(&rctx));
 91:   } else {
 92:     PetscCall(VecSet(u, 1.0));
 93:   }
 94:   PetscCall(MatMult(A, u, b));

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:                 Create the linear solver and set various options
 98:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 99:   /* Create linear solver context */
100:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

102:   /* Set operators. */
103:   PetscCall(KSPSetOperators(ksp, A, A));

105:   PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));

107:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-setfromoptions_first", &setfromoptions_first, NULL));
108:   if (setfromoptions_first) {
109:     /* code path for changing from KSPLSQR to KSPREONLY */
110:     PetscCall(KSPSetFromOptions(ksp));
111:   }
112:   PetscCall(KSPSetType(ksp, KSPPREONLY));
113:   PetscCall(KSPGetPC(ksp, &pc));
114:   PetscCall(PCSetType(pc, PCCHOLESKY));
115: #if defined(PETSC_HAVE_MUMPS)
116:   #if defined(PETSC_USE_COMPLEX)
117:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Spectrum slicing with MUMPS is not available for complex scalars");
118:   #endif
119:   PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
120:   /*
121:      must use runtime option '-mat_mumps_icntl_13 1' (turn off ScaLAPACK for
122:      matrix inertia), currently there is no better way of setting this in program
123:   */
124:   PetscCall(PetscOptionsInsertString(NULL, "-mat_mumps_icntl_13 1"));
125: #else
126:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
127:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Configure with MUMPS if you want to run this example in parallel");
128: #endif

130:   if (!setfromoptions_first) {
131:     /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
132:     PetscCall(KSPSetFromOptions(ksp));
133:   }

135:   /* get inertia */
136:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-solve", &solve, NULL));
137:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-sameA", &sameA, NULL));
138:   PetscCall(KSPSetUp(ksp));
139:   PetscCall(PCFactorGetMatrix(pc, &F));
140:   PetscCall(MatGetInertia(F, &in, NULL, NULL));
141:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in));
142:   if (solve) {
143:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving the intermediate KSP\n"));
144:     PetscCall(KSPSolve(ksp, b, x));
145:   } else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NOT Solving the intermediate KSP\n"));

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:                       Solve the linear system
149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
151:   if (sameA) {
152:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting A\n"));
153:     PetscCall(MatAXPY(A, 1.1, B, DIFFERENT_NONZERO_PATTERN));
154:     PetscCall(KSPSetOperators(ksp, A, A));
155:   } else {
156:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting B\n"));
157:     PetscCall(MatAXPY(B, 1.1, A, DIFFERENT_NONZERO_PATTERN));
158:     PetscCall(KSPSetOperators(ksp, B, B));
159:   }
160:   PetscCall(KSPSetUp(ksp));
161:   PetscCall(PCFactorGetMatrix(pc, &F));
162:   PetscCall(MatGetInertia(F, &in, NULL, NULL));
163:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in));
164:   PetscCall(KSPSolve(ksp, b, x));
165:   PetscCall(MatDestroy(&B));

167:   /* Free work space.*/
168:   PetscCall(KSPDestroy(&ksp));
169:   PetscCall(VecDestroy(&u));
170:   PetscCall(VecDestroy(&x));
171:   PetscCall(VecDestroy(&b));
172:   PetscCall(MatDestroy(&A));

174:   PetscCall(PetscFinalize());
175:   return 0;
176: }

178: /*TEST

180:     build:
181:       requires: !complex

183:     test:
184:       args:

186:     test:
187:       suffix: 2
188:       args: -sameA

190:     test:
191:       suffix: 3
192:       args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}

194: TEST*/