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 PetscScalar v[3] = {-1./2., 1., -1./2.}; 13 PetscInt j[3]; 14 PetscInt k=15; 15 PetscInt M,m=420; 16 17 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 18 PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL)); 19 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 20 21 PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp)); 22 PetscCall(KSPSetFromOptions(ksp)); 23 PetscCall(KSPGetOperators(ksp,&A,NULL)); 24 25 PetscCall(MatSetSizes(A,m,m,PETSC_DETERMINE,PETSC_DETERMINE)); 26 PetscCall(MatSetFromOptions(A)); 27 PetscCall(MatSetUp(A)); 28 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 29 PetscCall(MatGetSize(A,&M,NULL)); 30 for (Ii=Istart; Ii<Iend; Ii++) { 31 j[0] = Ii - k; 32 j[1] = Ii; 33 j[2] = (Ii + k) < M ? (Ii + k) : -1; 34 PetscCall(MatSetValues(A,1,&Ii,3,j,v,INSERT_VALUES)); 35 } 36 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 37 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 38 PetscCall(MatCreateVecs(A,&x,&b)); 39 40 PetscCall(VecSetFromOptions(b)); 41 PetscCall(VecSet(b,1.0)); 42 PetscCall(VecSetFromOptions(x)); 43 PetscCall(VecSet(x,2.0)); 44 45 PetscCall(KSPSolve(ksp,b,x)); 46 47 PetscCall(VecDestroy(&b)); 48 PetscCall(VecDestroy(&x)); 49 PetscCall(KSPDestroy(&ksp)); 50 51 PetscCall(PetscFinalize()); 52 return 0; 53 } 54