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