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