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