1 2 static char help[] = "Solves a tridiagonal linear system with KSP. \n\ 3 It illustrates how to do one symbolic factorization and multiple numeric factorizations using same matrix structure. \n\n"; 4 5 #include <petscksp.h> 6 int main(int argc,char **args) 7 { 8 Vec x, b, u; /* approx solution, RHS, exact solution */ 9 Mat A; /* linear system matrix */ 10 KSP ksp; /* linear solver context */ 11 PC pc; /* preconditioner context */ 12 PetscReal norm; /* norm of solution error */ 13 PetscInt i,col[3],its,rstart,rend,N=10,num_numfac; 14 PetscScalar value[3]; 15 16 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 17 PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 18 19 /* Create and assemble matrix. */ 20 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 21 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N)); 22 PetscCall(MatSetFromOptions(A)); 23 PetscCall(MatSetUp(A)); 24 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 25 26 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 27 for (i=rstart; i<rend; i++) { 28 col[0] = i-1; col[1] = i; col[2] = i+1; 29 if (i == 0) { 30 PetscCall(MatSetValues(A,1,&i,2,col+1,value+1,INSERT_VALUES)); 31 } else if (i == N-1) { 32 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 33 } else { 34 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 35 } 36 } 37 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 38 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 39 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 40 41 /* Create vectors */ 42 PetscCall(MatCreateVecs(A,&x,&b)); 43 PetscCall(VecDuplicate(x,&u)); 44 45 /* Set exact solution; then compute right-hand-side vector. */ 46 PetscCall(VecSet(u,1.0)); 47 PetscCall(MatMult(A,u,b)); 48 49 /* Create the linear solver and set various options. */ 50 PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp)); 51 PetscCall(KSPGetPC(ksp,&pc)); 52 PetscCall(PCSetType(pc,PCJACOBI)); 53 PetscCall(KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT)); 54 PetscCall(KSPSetOperators(ksp,A,A)); 55 PetscCall(KSPSetFromOptions(ksp)); 56 57 num_numfac = 1; 58 PetscCall(PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL)); 59 while (num_numfac--) { 60 /* An example on how to update matrix A for repeated numerical factorization and solve. */ 61 PetscScalar one=1.0; 62 PetscInt i = 0; 63 PetscCall(MatSetValues(A,1,&i,1,&i,&one,ADD_VALUES)); 64 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 65 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 66 /* Update b */ 67 PetscCall(MatMult(A,u,b)); 68 69 /* Solve the linear system */ 70 PetscCall(KSPSolve(ksp,b,x)); 71 72 /* Check the solution and clean up */ 73 PetscCall(VecAXPY(x,-1.0,u)); 74 PetscCall(VecNorm(x,NORM_2,&norm)); 75 PetscCall(KSPGetIterationNumber(ksp,&its)); 76 if (norm > 100*PETSC_MACHINE_EPSILON) { 77 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its)); 78 } 79 } 80 81 /* Free work space. */ 82 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&u)); 83 PetscCall(VecDestroy(&b)); PetscCall(MatDestroy(&A)); 84 PetscCall(KSPDestroy(&ksp)); 85 86 PetscCall(PetscFinalize()); 87 return 0; 88 } 89 90 /*TEST 91 92 test: 93 args: -num_numfac 2 -pc_type lu 94 95 test: 96 suffix: 2 97 args: -num_numfac 2 -pc_type lu -pc_factor_mat_solver_type mumps 98 requires: mumps 99 100 test: 101 suffix: 3 102 nsize: 3 103 args: -num_numfac 2 -pc_type lu -pc_factor_mat_solver_type mumps 104 requires: mumps 105 106 TEST*/ 107