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