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