1 static char help[] = "Tests CG, MINRES and SYMMLQ on symmetric matrices with SBAIJ format. The preconditioner ICC only works on sequential SBAIJ format. \n\n"; 2 3 #include <petscksp.h> 4 5 int main(int argc, char **args) 6 { 7 Mat C; 8 PetscScalar v, none = -1.0; 9 PetscInt i, j, Ii, J, Istart, Iend, N, m = 4, n = 4, its, k; 10 PetscMPIInt size, rank; 11 PetscReal err_norm, res_norm; 12 Vec x, b, u, u_tmp; 13 PC pc; 14 KSP ksp; 15 16 PetscFunctionBeginUser; 17 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 18 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 22 N = m * n; 23 24 /* Generate matrix */ 25 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 26 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N)); 27 PetscCall(MatSetFromOptions(C)); 28 PetscCall(MatSetUp(C)); 29 PetscCall(MatGetOwnershipRange(C, &Istart, &Iend)); 30 for (Ii = Istart; Ii < Iend; Ii++) { 31 v = -1.0; 32 i = Ii / n; 33 j = Ii - i * n; 34 if (i > 0) { 35 J = Ii - n; 36 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES)); 37 } 38 if (i < m - 1) { 39 J = Ii + n; 40 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES)); 41 } 42 if (j > 0) { 43 J = Ii - 1; 44 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES)); 45 } 46 if (j < n - 1) { 47 J = Ii + 1; 48 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES)); 49 } 50 v = 4.0; 51 PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, ADD_VALUES)); 52 } 53 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 54 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 55 56 /* a shift can make C indefinite. Preconditioners LU, ILU (for BAIJ format) and ICC may fail */ 57 /* PetscCall(MatShift(C,alpha)); */ 58 /* PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); */ 59 60 /* Setup and solve for system */ 61 /* Create vectors. */ 62 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 63 PetscCall(VecSetSizes(x, PETSC_DECIDE, N)); 64 PetscCall(VecSetFromOptions(x)); 65 PetscCall(VecDuplicate(x, &b)); 66 PetscCall(VecDuplicate(x, &u)); 67 PetscCall(VecDuplicate(x, &u_tmp)); 68 /* Set exact solution u; then compute right-hand-side vector b. */ 69 PetscCall(VecSet(u, 1.0)); 70 PetscCall(MatMult(C, u, b)); 71 72 for (k = 0; k < 3; k++) { 73 if (k == 0) { /* CG */ 74 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 75 PetscCall(KSPSetOperators(ksp, C, C)); 76 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n CG: \n")); 77 PetscCall(KSPSetType(ksp, KSPCG)); 78 } else if (k == 1) { /* MINRES */ 79 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 80 PetscCall(KSPSetOperators(ksp, C, C)); 81 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MINRES: \n")); 82 PetscCall(KSPSetType(ksp, KSPMINRES)); 83 } else { /* SYMMLQ */ 84 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 85 PetscCall(KSPSetOperators(ksp, C, C)); 86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n SYMMLQ: \n")); 87 PetscCall(KSPSetType(ksp, KSPSYMMLQ)); 88 } 89 PetscCall(KSPGetPC(ksp, &pc)); 90 /* PetscCall(PCSetType(pc,PCICC)); */ 91 PetscCall(PCSetType(pc, PCJACOBI)); 92 PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT)); 93 94 /* 95 Set runtime options, e.g., 96 -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> 97 These options will override those specified above as long as 98 KSPSetFromOptions() is called _after_ any other customization 99 routines. 100 */ 101 PetscCall(KSPSetFromOptions(ksp)); 102 103 /* Solve linear system; */ 104 PetscCall(KSPSetUp(ksp)); 105 PetscCall(KSPSolve(ksp, b, x)); 106 107 PetscCall(KSPGetIterationNumber(ksp, &its)); 108 /* Check error */ 109 PetscCall(VecCopy(u, u_tmp)); 110 PetscCall(VecAXPY(u_tmp, none, x)); 111 PetscCall(VecNorm(u_tmp, NORM_2, &err_norm)); 112 PetscCall(MatMult(C, x, u_tmp)); 113 PetscCall(VecAXPY(u_tmp, none, b)); 114 PetscCall(VecNorm(u_tmp, NORM_2, &res_norm)); 115 116 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its)); 117 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g;", (double)res_norm)); 118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm %g.\n", (double)err_norm)); 119 PetscCall(KSPDestroy(&ksp)); 120 } 121 122 /* 123 Free work space. All PETSc objects should be destroyed when they 124 are no longer needed. 125 */ 126 PetscCall(VecDestroy(&b)); 127 PetscCall(VecDestroy(&u)); 128 PetscCall(VecDestroy(&x)); 129 PetscCall(VecDestroy(&u_tmp)); 130 PetscCall(MatDestroy(&C)); 131 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 138 test: 139 args: -pc_type icc -mat_type seqsbaij -mat_ignore_lower_triangular 140 141 test: 142 suffix: 2 143 args: -pc_type icc -pc_factor_levels 2 -mat_type seqsbaij -mat_ignore_lower_triangular 144 145 test: 146 suffix: 3 147 nsize: 2 148 args: -pc_type bjacobi -sub_pc_type icc -mat_type mpisbaij -mat_ignore_lower_triangular -ksp_max_it 8 149 150 test: 151 suffix: 4 152 nsize: 2 153 args: -pc_type bjacobi -sub_pc_type icc -sub_pc_factor_levels 1 -mat_type mpisbaij -mat_ignore_lower_triangular 154 155 TEST*/ 156