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 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 19 PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL)); 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 21 22 PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp)); 23 PetscCall(KSPSetFromOptions(ksp)); 24 PetscCall(KSPGetOperators(ksp,&A,NULL)); 25 26 PetscCall(MatSetSizes(A,m,m,PETSC_DETERMINE,PETSC_DETERMINE)); 27 PetscCall(MatSetFromOptions(A)); 28 PetscCall(MatSetUp(A)); 29 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 30 PetscCall(MatGetSize(A,&M,NULL)); 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 PetscCall(MatSetValues(A,1,&Ii,3,j,v,INSERT_VALUES)); 36 } 37 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 38 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 39 PetscCall(MatCreateVecs(A,&x,&b)); 40 41 PetscCall(VecSetFromOptions(b)); 42 PetscCall(VecSet(b,1.0)); 43 PetscCall(VecSetFromOptions(x)); 44 PetscCall(VecSet(x,2.0)); 45 46 PetscCall(KSPSolve(ksp,b,x)); 47 48 PetscCall(VecDestroy(&b)); 49 PetscCall(VecDestroy(&x)); 50 PetscCall(KSPDestroy(&ksp)); 51 52 PetscCall(PetscFinalize()); 53 return 0; 54 } 55