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 PetscErrorCode ierr; 14 PetscInt i,col[3],its,rstart,rend,N=10,num_numfac; 15 PetscScalar value[3]; 16 17 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 18 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 19 20 /* Create and assemble matrix. */ 21 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 22 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); 23 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 24 ierr = MatSetUp(A);CHKERRQ(ierr); 25 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,2,col+1,value+1,INSERT_VALUES);CHKERRQ(ierr); 32 } else if (i == N-1) { 33 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 34 } else { 35 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 36 } 37 } 38 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 39 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 40 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 41 42 /* Create vectors */ 43 ierr = MatCreateVecs(A,&x,&b);CHKERRQ(ierr); 44 ierr = VecDuplicate(x,&u);CHKERRQ(ierr); 45 46 /* Set exact solution; then compute right-hand-side vector. */ 47 ierr = VecSet(u,1.0);CHKERRQ(ierr); 48 ierr = MatMult(A,u,b);CHKERRQ(ierr); 49 50 /* Create the linear solver and set various options. */ 51 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 52 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 53 ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr); 54 ierr = KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 55 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 56 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 57 58 num_numfac = 1; 59 ierr = PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,1,&i,&one,ADD_VALUES);CHKERRQ(ierr); 65 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 66 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 67 /* Update b */ 68 ierr = MatMult(A,u,b);CHKERRQ(ierr); 69 70 /* Solve the linear system */ 71 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 72 73 /* Check the solution and clean up */ 74 ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr); 75 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 76 ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 77 if (norm > 100*PETSC_MACHINE_EPSILON) { 78 ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr); 79 } 80 } 81 82 /* Free work space. */ 83 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr); 84 ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); 85 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 86 87 ierr = PetscFinalize(); 88 return ierr; 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