1 2 static char help[] = "Demonstrates the use of fast Richardson for SOR. And\n\ 3 also tests the MatSOR() routines. Input parameters are:\n\ 4 -n <n> : problem dimension\n\n"; 5 6 #include <petscksp.h> 7 #include <petscpc.h> 8 9 int main(int argc, char **args) 10 { 11 Mat mat; /* matrix */ 12 Vec b, ustar, u; /* vectors (RHS, exact solution, approx solution) */ 13 PC pc; /* PC context */ 14 KSP ksp; /* KSP context */ 15 PetscInt n = 10, i, its, col[3]; 16 PetscScalar value[3]; 17 KSPType kspname; 18 PCType pcname; 19 20 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 23 24 /* Create and initialize vectors */ 25 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &b)); 26 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &ustar)); 27 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &u)); 28 PetscCall(VecSet(ustar, 1.0)); 29 PetscCall(VecSet(u, 0.0)); 30 31 /* Create and assemble matrix */ 32 PetscCall(MatCreate(PETSC_COMM_SELF, &mat)); 33 PetscCall(MatSetType(mat, MATSEQAIJ)); 34 PetscCall(MatSetSizes(mat, n, n, n, n)); 35 PetscCall(MatSetFromOptions(mat)); 36 PetscCall(MatSeqAIJSetPreallocation(mat, 3, NULL)); 37 PetscCall(MatSeqBAIJSetPreallocation(mat, 1, 3, NULL)); 38 PetscCall(MatSeqSBAIJSetPreallocation(mat, 1, 3, NULL)); 39 PetscCall(MatSeqSELLSetPreallocation(mat, 3, NULL)); 40 value[0] = -1.0; 41 value[1] = 2.0; 42 value[2] = -1.0; 43 for (i = 1; i < n - 1; i++) { 44 col[0] = i - 1; 45 col[1] = i; 46 col[2] = i + 1; 47 PetscCall(MatSetValues(mat, 1, &i, 3, col, value, INSERT_VALUES)); 48 } 49 i = n - 1; 50 col[0] = n - 2; 51 col[1] = n - 1; 52 PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES)); 53 i = 0; 54 col[0] = 0; 55 col[1] = 1; 56 value[0] = 2.0; 57 value[1] = -1.0; 58 PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES)); 59 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 60 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 61 62 /* Compute right-hand-side vector */ 63 PetscCall(MatMult(mat, ustar, b)); 64 65 /* Create PC context and set up data structures */ 66 PetscCall(PCCreate(PETSC_COMM_WORLD, &pc)); 67 PetscCall(PCSetType(pc, PCNONE)); 68 PetscCall(PCSetFromOptions(pc)); 69 PetscCall(PCSetOperators(pc, mat, mat)); 70 PetscCall(PCSetUp(pc)); 71 72 /* Create KSP context and set up data structures */ 73 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 74 PetscCall(KSPSetType(ksp, KSPRICHARDSON)); 75 PetscCall(KSPSetFromOptions(ksp)); 76 PetscCall(PCSetOperators(pc, mat, mat)); 77 PetscCall(KSPSetPC(ksp, pc)); 78 PetscCall(KSPSetUp(ksp)); 79 80 /* Solve the problem */ 81 PetscCall(KSPGetType(ksp, &kspname)); 82 PetscCall(PCGetType(pc, &pcname)); 83 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Running %s with %s preconditioning\n", kspname, pcname)); 84 PetscCall(KSPSolve(ksp, b, u)); 85 PetscCall(KSPGetIterationNumber(ksp, &its)); 86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations %" PetscInt_FMT "\n", its)); 87 88 /* Free data structures */ 89 PetscCall(KSPDestroy(&ksp)); 90 PetscCall(VecDestroy(&u)); 91 PetscCall(VecDestroy(&ustar)); 92 PetscCall(VecDestroy(&b)); 93 PetscCall(MatDestroy(&mat)); 94 PetscCall(PCDestroy(&pc)); 95 PetscCall(PetscFinalize()); 96 return 0; 97 } 98 99 /*TEST 100 101 testset: 102 args: -ksp_type gmres -ksp_monitor_short -pc_type sor -pc_sor_symmetric 103 output_file: output/ex3_1.out 104 test: 105 suffix: sor_aij 106 test: 107 suffix: sor_seqbaij 108 args: -mat_type seqbaij 109 test: 110 suffix: sor_seqsbaij 111 args: -mat_type seqbaij 112 test: 113 suffix: sor_seqsell 114 args: -mat_type seqsell 115 116 TEST*/ 117