1 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"; 4 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 25 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)); 40 41 PetscCall(MatGetOwnershipRange(A, &Istart, &Iend)); 42 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()); 71 72 /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */ 73 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 74 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)); 81 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)); 98 99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 100 Create the linear solver and set various options 101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 102 /* Create linear solver context */ 103 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 104 105 /* Set operators. */ 106 PetscCall(KSPSetOperators(ksp, A, A)); 107 108 PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 109 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 132 133 if (!setfromoptions_first) { 134 /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */ 135 PetscCall(KSPSetFromOptions(ksp)); 136 } 137 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")); 149 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)); 169 170 /* Free work space.*/ 171 PetscCall(KSPDestroy(&ksp)); 172 PetscCall(VecDestroy(&u)); 173 PetscCall(VecDestroy(&x)); 174 PetscCall(VecDestroy(&b)); 175 PetscCall(MatDestroy(&A)); 176 177 PetscCall(PetscFinalize()); 178 return 0; 179 } 180 181 /*TEST 182 183 build: 184 requires: !complex 185 186 test: 187 args: 188 189 test: 190 suffix: 2 191 args: -sameA 192 193 test: 194 suffix: 3 195 args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output} 196 197 TEST*/ 198