1 2 static char help[] = "Solves a linear system in parallel with KSP. \n\ 3 Contributed by Jose E. Roman, SLEPc developer, for testing repeated call of KSPSetOperators(), 2014 \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 PetscRandom rctx; /* random number generator context */ 12 PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7; 13 PetscErrorCode ierr; 14 PetscBool flg = PETSC_FALSE; 15 PetscScalar v; 16 PC pc; 17 PetscInt in; 18 Mat F,B; 19 PetscBool solve=PETSC_FALSE,sameA=PETSC_FALSE; 20 #if defined(PETSC_USE_LOG) 21 PetscLogStage stage; 22 #endif 23 #if !defined(PETSC_HAVE_MUMPS) 24 PetscMPIInt size; 25 #endif 26 27 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 28 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 30 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 31 Compute the matrix and right-hand-side vector that define 32 the linear system, Ax = b. 33 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 34 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 35 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); 36 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 37 ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr); 38 ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr); 39 ierr = MatSetUp(A);CHKERRQ(ierr); 40 41 ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); 42 43 ierr = PetscLogStageRegister("Assembly", &stage);CHKERRQ(ierr); 44 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 45 for (Ii=Istart; Ii<Iend; Ii++) { 46 v = -1.0; i = Ii/n; j = Ii - i*n; 47 if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 48 if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 49 if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 50 if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 51 v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr); 52 } 53 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 54 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 55 ierr = PetscLogStagePop();CHKERRQ(ierr); 56 57 /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */ 58 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 59 60 /* Create parallel vectors. */ 61 ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr); 62 ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr); 63 ierr = VecSetFromOptions(u);CHKERRQ(ierr); 64 ierr = VecDuplicate(u,&b);CHKERRQ(ierr); 65 ierr = VecDuplicate(b,&x);CHKERRQ(ierr); 66 67 /* 68 Set exact solution; then compute right-hand-side vector. 69 By default we use an exact solution of a vector with all 70 elements of 1.0; Alternatively, using the runtime option 71 -random_sol forms a solution vector with random components. 72 */ 73 ierr = PetscOptionsGetBool(NULL,NULL,"-random_exact_sol",&flg,NULL);CHKERRQ(ierr); 74 if (flg) { 75 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr); 76 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 77 ierr = VecSetRandom(u,rctx);CHKERRQ(ierr); 78 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 79 } else { 80 ierr = VecSet(u,1.0);CHKERRQ(ierr); 81 } 82 ierr = MatMult(A,u,b);CHKERRQ(ierr); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Create the linear solver and set various options 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 /* Create linear solver context */ 88 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 89 90 /* Set operators. */ 91 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 92 93 ierr = KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 94 95 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 96 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 97 ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr); 98 #if defined(PETSC_HAVE_MUMPS) 99 #if defined(PETSC_USE_COMPLEX) 100 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Spectrum slicing with MUMPS is not available for complex scalars"); 101 #endif 102 ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);CHKERRQ(ierr); 103 /* 104 must use runtime option '-mat_mumps_icntl_13 1' (turn off scaLAPACK for 105 matrix inertia), currently there is no better way of setting this in program 106 */ 107 ierr = PetscOptionsInsertString(NULL,"-mat_mumps_icntl_13 1");CHKERRQ(ierr); 108 #else 109 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 110 if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Configure with MUMPS if you want to run this example in parallel"); 111 #endif 112 113 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 114 115 /* get inertia */ 116 ierr = PetscOptionsGetBool(NULL,NULL,"-solve",&solve,NULL);CHKERRQ(ierr); 117 ierr = PetscOptionsGetBool(NULL,NULL,"-sameA",&sameA,NULL);CHKERRQ(ierr); 118 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 119 ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr); 120 ierr = MatGetInertia(F,&in,NULL,NULL);CHKERRQ(ierr); 121 ierr = PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);CHKERRQ(ierr); 122 if (solve) { 123 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solving the intermediate KSP\n");CHKERRQ(ierr); 124 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 125 } else {ierr = PetscPrintf(PETSC_COMM_WORLD,"NOT Solving the intermediate KSP\n");CHKERRQ(ierr);} 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Solve the linear system 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 131 if (sameA) { 132 ierr = PetscPrintf(PETSC_COMM_WORLD,"Seting A\n");CHKERRQ(ierr); 133 ierr = MatAXPY(A,1.1,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 134 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 135 } else { 136 ierr = PetscPrintf(PETSC_COMM_WORLD,"Seting B\n");CHKERRQ(ierr); 137 ierr = MatAXPY(B,1.1,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 138 ierr = KSPSetOperators(ksp,B,B);CHKERRQ(ierr); 139 } 140 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 141 ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr); 142 ierr = MatGetInertia(F,&in,NULL,NULL);CHKERRQ(ierr); 143 ierr = PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);CHKERRQ(ierr); 144 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 145 ierr = MatDestroy(&B);CHKERRQ(ierr); 146 147 /* Free work space.*/ 148 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 149 ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); 150 ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); 151 152 ierr = PetscFinalize(); 153 return ierr; 154 } 155 156 /*TEST 157 158 build: 159 requires: !complex 160 161 test: 162 args: 163 164 test: 165 suffix: 2 166 args: -sameA 167 168 TEST*/ 169