Actual source code: ex8.c


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

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

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

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

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

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

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

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

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

105:   /* Set operators. */
106:   PetscCall(KSPSetOperators(ksp, A, A));

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

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

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

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

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

170:   /* Free work space.*/
171:   PetscCall(KSPDestroy(&ksp));
172:   PetscCall(VecDestroy(&u));
173:   PetscCall(VecDestroy(&x));
174:   PetscCall(VecDestroy(&b));
175:   PetscCall(MatDestroy(&A));

177:   PetscCall(PetscFinalize());
178:   return 0;
179: }

181: /*TEST

183:     build:
184:       requires: !complex

186:     test:
187:       args:

189:     test:
190:       suffix: 2
191:       args: -sameA

193:     test:
194:       suffix: 3
195:       args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}

197: TEST*/