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 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 #if defined(PETSC_USE_LOG) 19 PetscLogStage stage; 20 #endif 21 #if !defined(PETSC_HAVE_MUMPS) 22 PetscMPIInt size; 23 #endif 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 29 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 30 Compute the matrix and right-hand-side vector that define 31 the linear system, Ax = b. 32 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 33 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 34 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n)); 35 PetscCall(MatSetFromOptions(A)); 36 PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL)); 37 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL)); 38 PetscCall(MatSetUp(A)); 39 40 PetscCall(MatGetOwnershipRange(A, &Istart, &Iend)); 41 42 PetscCall(PetscLogStageRegister("Assembly", &stage)); 43 PetscCall(PetscLogStagePush(stage)); 44 for (Ii = Istart; Ii < Iend; Ii++) { 45 v = -1.0; 46 i = Ii / n; 47 j = Ii - i * n; 48 if (i > 0) { 49 J = Ii - n; 50 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 51 } 52 if (i < m - 1) { 53 J = Ii + n; 54 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 55 } 56 if (j > 0) { 57 J = Ii - 1; 58 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 59 } 60 if (j < n - 1) { 61 J = Ii + 1; 62 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 63 } 64 v = 4.0; 65 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 66 } 67 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 68 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 69 PetscCall(PetscLogStagePop()); 70 71 /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */ 72 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 73 74 /* Create parallel vectors. */ 75 PetscCall(VecCreate(PETSC_COMM_WORLD, &u)); 76 PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n)); 77 PetscCall(VecSetFromOptions(u)); 78 PetscCall(VecDuplicate(u, &b)); 79 PetscCall(VecDuplicate(b, &x)); 80 81 /* 82 Set exact solution; then compute right-hand-side vector. 83 By default we use an exact solution of a vector with all 84 elements of 1.0; Alternatively, using the runtime option 85 -random_sol forms a solution vector with random components. 86 */ 87 PetscCall(PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL)); 88 if (flg) { 89 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 90 PetscCall(PetscRandomSetFromOptions(rctx)); 91 PetscCall(VecSetRandom(u, rctx)); 92 PetscCall(PetscRandomDestroy(&rctx)); 93 } else { 94 PetscCall(VecSet(u, 1.0)); 95 } 96 PetscCall(MatMult(A, u, b)); 97 98 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 Create the linear solver and set various options 100 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 101 /* Create linear solver context */ 102 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 103 104 /* Set operators. */ 105 PetscCall(KSPSetOperators(ksp, A, A)); 106 107 PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 108 109 PetscCall(PetscOptionsGetBool(NULL, NULL, "-setfromoptions_first", &setfromoptions_first, NULL)); 110 if (setfromoptions_first) { 111 /* code path for changing from KSPLSQR to KSPREONLY */ 112 PetscCall(KSPSetFromOptions(ksp)); 113 } 114 PetscCall(KSPSetType(ksp, KSPPREONLY)); 115 PetscCall(KSPGetPC(ksp, &pc)); 116 PetscCall(PCSetType(pc, PCCHOLESKY)); 117 #if defined(PETSC_HAVE_MUMPS) 118 #if defined(PETSC_USE_COMPLEX) 119 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Spectrum slicing with MUMPS is not available for complex scalars"); 120 #endif 121 PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 122 /* 123 must use runtime option '-mat_mumps_icntl_13 1' (turn off ScaLAPACK for 124 matrix inertia), currently there is no better way of setting this in program 125 */ 126 PetscCall(PetscOptionsInsertString(NULL, "-mat_mumps_icntl_13 1")); 127 #else 128 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 129 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Configure with MUMPS if you want to run this example in parallel"); 130 #endif 131 132 if (!setfromoptions_first) { 133 /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */ 134 PetscCall(KSPSetFromOptions(ksp)); 135 } 136 137 /* get inertia */ 138 PetscCall(PetscOptionsGetBool(NULL, NULL, "-solve", &solve, NULL)); 139 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sameA", &sameA, NULL)); 140 PetscCall(KSPSetUp(ksp)); 141 PetscCall(PCFactorGetMatrix(pc, &F)); 142 PetscCall(MatGetInertia(F, &in, NULL, NULL)); 143 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in)); 144 if (solve) { 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving the intermediate KSP\n")); 146 PetscCall(KSPSolve(ksp, b, x)); 147 } else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NOT Solving the intermediate KSP\n")); 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Solve the linear system 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 153 if (sameA) { 154 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting A\n")); 155 PetscCall(MatAXPY(A, 1.1, B, DIFFERENT_NONZERO_PATTERN)); 156 PetscCall(KSPSetOperators(ksp, A, A)); 157 } else { 158 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting B\n")); 159 PetscCall(MatAXPY(B, 1.1, A, DIFFERENT_NONZERO_PATTERN)); 160 PetscCall(KSPSetOperators(ksp, B, B)); 161 } 162 PetscCall(KSPSetUp(ksp)); 163 PetscCall(PCFactorGetMatrix(pc, &F)); 164 PetscCall(MatGetInertia(F, &in, NULL, NULL)); 165 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in)); 166 PetscCall(KSPSolve(ksp, b, x)); 167 PetscCall(MatDestroy(&B)); 168 169 /* Free work space.*/ 170 PetscCall(KSPDestroy(&ksp)); 171 PetscCall(VecDestroy(&u)); 172 PetscCall(VecDestroy(&x)); 173 PetscCall(VecDestroy(&b)); 174 PetscCall(MatDestroy(&A)); 175 176 PetscCall(PetscFinalize()); 177 return 0; 178 } 179 180 /*TEST 181 182 build: 183 requires: !complex 184 185 test: 186 args: 187 188 test: 189 suffix: 2 190 args: -sameA 191 192 test: 193 suffix: 3 194 args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output} 195 196 TEST*/ 197