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