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