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