1 2 static char help[] = "Solves a tridiagonal linear system. Designed to compare SOR for different Mat impls.\n\n"; 3 4 #include <petscksp.h> 5 6 int main(int argc,char **args) 7 { 8 KSP ksp; /* linear solver context */ 9 Mat A; /* linear system matrix */ 10 Vec x,b; /* approx solution, RHS */ 11 PetscInt Ii,Istart,Iend; 12 PetscErrorCode ierr; 13 PetscScalar v[3] = {-1./2., 1., -1./2.}; 14 PetscInt j[3]; 15 PetscInt k=15; 16 PetscInt M,m=420; 17 18 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 19 ierr = PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);CHKERRQ(ierr); 20 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 21 22 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 23 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 24 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); 25 26 ierr = MatSetSizes(A,m,m,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 27 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 28 ierr = MatSetUp(A);CHKERRQ(ierr); 29 ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); 30 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 31 for (Ii=Istart; Ii<Iend; Ii++) { 32 j[0] = Ii - k; 33 j[1] = Ii; 34 j[2] = (Ii + k) < M ? (Ii + k) : -1; 35 ierr = MatSetValues(A,1,&Ii,3,j,v,INSERT_VALUES);CHKERRQ(ierr); 36 } 37 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 38 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 39 ierr = MatCreateVecs(A,&x,&b);CHKERRQ(ierr); 40 41 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 42 ierr = VecSet(b,1.0);CHKERRQ(ierr); 43 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 44 ierr = VecSet(x,2.0);CHKERRQ(ierr); 45 46 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 47 48 ierr = VecDestroy(&b);CHKERRQ(ierr); 49 ierr = VecDestroy(&x);CHKERRQ(ierr); 50 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 51 52 ierr = PetscFinalize(); 53 return ierr; 54 } 55