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

 27:   PetscInitialize(&argc, &args, (char *)0, help);
 28:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 29:   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:   MatCreate(PETSC_COMM_WORLD, &A);
 35:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
 36:   MatSetFromOptions(A);
 37:   MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
 38:   MatSeqAIJSetPreallocation(A, 5, NULL);
 39:   MatSetUp(A);

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

 43:   PetscLogStageRegister("Assembly", &stage);
 44:   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:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 52:     }
 53:     if (i < m - 1) {
 54:       J = Ii + n;
 55:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 56:     }
 57:     if (j > 0) {
 58:       J = Ii - 1;
 59:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 60:     }
 61:     if (j < n - 1) {
 62:       J = Ii + 1;
 63:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 64:     }
 65:     v = 4.0;
 66:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
 67:   }
 68:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 70:   PetscLogStagePop();

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

 75:   /* Create parallel vectors. */
 76:   VecCreate(PETSC_COMM_WORLD, &u);
 77:   VecSetSizes(u, PETSC_DECIDE, m * n);
 78:   VecSetFromOptions(u);
 79:   VecDuplicate(u, &b);
 80:   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:   PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL);
 89:   if (flg) {
 90:     PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
 91:     PetscRandomSetFromOptions(rctx);
 92:     VecSetRandom(u, rctx);
 93:     PetscRandomDestroy(&rctx);
 94:   } else {
 95:     VecSet(u, 1.0);
 96:   }
 97:   MatMult(A, u, b);

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

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

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

110:   PetscOptionsGetBool(NULL, NULL, "-setfromoptions_first", &setfromoptions_first, NULL);
111:   if (setfromoptions_first) {
112:     /* code path for changing from KSPLSQR to KSPREONLY */
113:     KSPSetFromOptions(ksp);
114:   }
115:   KSPSetType(ksp, KSPPREONLY);
116:   KSPGetPC(ksp, &pc);
117:   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:   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:   PetscOptionsInsertString(NULL, "-mat_mumps_icntl_13 1");
128: #else
129:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
131: #endif

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

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

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

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

177:   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*/