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,setfromoptions_first=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 = PetscOptionsGetBool(NULL,NULL,"-setfromoptions_first",&setfromoptions_first,NULL);CHKERRQ(ierr); 96 if (setfromoptions_first) { 97 /* code path for changing from KSPLSQR to KSPREONLY */ 98 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 99 } 100 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 101 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 102 ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr); 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 ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);CHKERRQ(ierr); 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 ierr = PetscOptionsInsertString(NULL,"-mat_mumps_icntl_13 1");CHKERRQ(ierr); 113 #else 114 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 115 if (size>1) SETERRQ(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 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 121 } 122 123 /* get inertia */ 124 ierr = PetscOptionsGetBool(NULL,NULL,"-solve",&solve,NULL);CHKERRQ(ierr); 125 ierr = PetscOptionsGetBool(NULL,NULL,"-sameA",&sameA,NULL);CHKERRQ(ierr); 126 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 127 ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr); 128 ierr = MatGetInertia(F,&in,NULL,NULL);CHKERRQ(ierr); 129 ierr = PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);CHKERRQ(ierr); 130 if (solve) { 131 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solving the intermediate KSP\n");CHKERRQ(ierr); 132 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 133 } else {ierr = PetscPrintf(PETSC_COMM_WORLD,"NOT Solving the intermediate KSP\n");CHKERRQ(ierr);} 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Solve the linear system 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 139 if (sameA) { 140 ierr = PetscPrintf(PETSC_COMM_WORLD,"Setting A\n");CHKERRQ(ierr); 141 ierr = MatAXPY(A,1.1,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 142 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 143 } else { 144 ierr = PetscPrintf(PETSC_COMM_WORLD,"Setting B\n");CHKERRQ(ierr); 145 ierr = MatAXPY(B,1.1,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 146 ierr = KSPSetOperators(ksp,B,B);CHKERRQ(ierr); 147 } 148 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 149 ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr); 150 ierr = MatGetInertia(F,&in,NULL,NULL);CHKERRQ(ierr); 151 ierr = PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);CHKERRQ(ierr); 152 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 153 ierr = MatDestroy(&B);CHKERRQ(ierr); 154 155 /* Free work space.*/ 156 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 157 ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); 158 ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); 159 160 ierr = PetscFinalize(); 161 return ierr; 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