1 2 static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\ 3 This version first preloads and solves a small system, then loads \n\ 4 another (larger) system and solves it as well. This example illustrates\n\ 5 preloading of instructions with the smaller system so that more accurate\n\ 6 performance monitoring can be done with the larger one (that actually\n\ 7 is the system of interest). See the 'Performance Hints' chapter of the\n\ 8 users manual for a discussion of preloading. Input parameters include\n\ 9 -f0 <input_file> : first file to load (small system)\n\ 10 -f1 <input_file> : second file to load (larger system)\n\n\ 11 -nearnulldim <0> : number of vectors in the near-null space immediately following matrix\n\n\ 12 -trans : solve transpose system instead\n\n"; 13 /* 14 This code can be used to test PETSc interface to other packages.\n\ 15 Examples of command line options: \n\ 16 ./ex72 -f0 <datafile> -ksp_type preonly \n\ 17 -help -ksp_view \n\ 18 -num_numfac <num_numfac> -num_rhs <num_rhs> \n\ 19 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu or superlu_dist or mumps \n\ 20 -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_type mumps \n\ 21 mpiexec -n <np> ./ex72 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij 22 \n\n"; 23 */ 24 /*T 25 Concepts: KSP^solving a linear system 26 Processors: n 27 T*/ 28 29 /* 30 Include "petscksp.h" so that we can use KSP solvers. Note that this file 31 automatically includes: 32 petscsys.h - base PETSc routines petscvec.h - vectors 33 petscmat.h - matrices 34 petscis.h - index sets petscksp.h - Krylov subspace methods 35 petscviewer.h - viewers petscpc.h - preconditioners 36 */ 37 #include <petscksp.h> 38 39 int main(int argc,char **args) 40 { 41 KSP ksp; /* linear solver context */ 42 Mat A; /* matrix */ 43 Vec x,b,u; /* approx solution, RHS, exact solution */ 44 PetscViewer viewer; /* viewer */ 45 char file[4][PETSC_MAX_PATH_LEN]; /* input file name */ 46 PetscBool table =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE; 47 PetscBool outputSoln=PETSC_FALSE,constantnullspace = PETSC_FALSE; 48 PetscErrorCode ierr; 49 PetscInt its,num_numfac,m,n,M,nearnulldim = 0; 50 PetscReal norm; 51 PetscBool preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE; 52 PetscMPIInt rank; 53 char initialguessfilename[PETSC_MAX_PATH_LEN]; 54 55 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 56 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 57 ierr = PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);CHKERRQ(ierr); 58 ierr = PetscOptionsGetBool(NULL,NULL,"-constantnullspace",&constantnullspace,NULL);CHKERRQ(ierr); 59 ierr = PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsGetBool(NULL,NULL,"-initialguess",&initialguess,NULL);CHKERRQ(ierr); 61 ierr = PetscOptionsGetBool(NULL,NULL,"-output_solution",&outputSoln,NULL);CHKERRQ(ierr); 62 ierr = PetscOptionsGetString(NULL,NULL,"-initialguessfilename",initialguessfilename,sizeof(initialguessfilename),&initialguessfile);CHKERRQ(ierr); 63 ierr = PetscOptionsGetInt(NULL,NULL,"-nearnulldim",&nearnulldim,NULL);CHKERRQ(ierr); 64 65 /* 66 Determine files from which we read the two linear systems 67 (matrix and right-hand-side vector). 68 */ 69 ierr = PetscOptionsGetString(NULL,NULL,"-f",file[0],sizeof(file[0]),&flg);CHKERRQ(ierr); 70 if (flg) { 71 ierr = PetscStrcpy(file[1],file[0]);CHKERRQ(ierr); 72 preload = PETSC_FALSE; 73 } else { 74 ierr = PetscOptionsGetString(NULL,NULL,"-f0",file[0],sizeof(file[0]),&flg);CHKERRQ(ierr); 75 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate binary file with the -f0 or -f option"); 76 ierr = PetscOptionsGetString(NULL,NULL,"-f1",file[1],sizeof(file[1]),&flg);CHKERRQ(ierr); 77 if (!flg) preload = PETSC_FALSE; /* don't bother with second system */ 78 } 79 80 /* ----------------------------------------------------------- 81 Beginning of linear solver loop 82 ----------------------------------------------------------- */ 83 /* 84 Loop through the linear solve 2 times. 85 - The intention here is to preload and solve a small system; 86 then load another (larger) system and solve it as well. 87 This process preloads the instructions with the smaller 88 system so that more accurate performance monitoring (via 89 -log_view) can be done with the larger one (that actually 90 is the system of interest). 91 */ 92 PetscPreLoadBegin(preload,"Load system"); 93 94 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 95 Load system 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 98 /* 99 Open binary file. Note that we use FILE_MODE_READ to indicate 100 reading from this file. 101 */ 102 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&viewer);CHKERRQ(ierr); 103 104 /* 105 Load the matrix and vector; then destroy the viewer. 106 */ 107 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 108 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 109 ierr = MatLoad(A,viewer);CHKERRQ(ierr); 110 if (nearnulldim) { 111 MatNullSpace nullsp; 112 Vec *nullvecs; 113 PetscInt i; 114 ierr = PetscMalloc1(nearnulldim,&nullvecs);CHKERRQ(ierr); 115 for (i=0; i<nearnulldim; i++) { 116 ierr = VecCreate(PETSC_COMM_WORLD,&nullvecs[i]);CHKERRQ(ierr); 117 ierr = VecLoad(nullvecs[i],viewer);CHKERRQ(ierr); 118 } 119 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nearnulldim,nullvecs,&nullsp);CHKERRQ(ierr); 120 ierr = MatSetNearNullSpace(A,nullsp);CHKERRQ(ierr); 121 for (i=0; i<nearnulldim; i++) {ierr = VecDestroy(&nullvecs[i]);CHKERRQ(ierr);} 122 ierr = PetscFree(nullvecs);CHKERRQ(ierr); 123 ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr); 124 } 125 if (constantnullspace) { 126 MatNullSpace constant; 127 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constant);CHKERRQ(ierr); 128 ierr = MatSetNullSpace(A,constant);CHKERRQ(ierr); 129 ierr = MatNullSpaceDestroy(&constant);CHKERRQ(ierr); 130 } 131 flg = PETSC_FALSE; 132 ierr = PetscOptionsGetString(NULL,NULL,"-rhs",file[2],sizeof(file[2]),&flg);CHKERRQ(ierr); 133 ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr); 134 if (flg) { /* rhs is stored in a separate file */ 135 if (file[2][0] == '0' || file[2][0] == 0) { 136 PetscInt m; 137 PetscScalar one = 1.0; 138 ierr = PetscInfo(0,"Using vector of ones for RHS\n");CHKERRQ(ierr); 139 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 140 ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr); 141 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 142 ierr = VecSet(b,one);CHKERRQ(ierr); 143 } else { 144 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 145 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&viewer);CHKERRQ(ierr); 146 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 147 ierr = VecLoad(b,viewer);CHKERRQ(ierr); 148 } 149 } else { /* rhs is stored in the same file as matrix */ 150 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 151 ierr = VecLoad(b,viewer);CHKERRQ(ierr); 152 } 153 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 154 155 /* Make A singular for testing zero-pivot of ilu factorization */ 156 /* Example: ./ex72 -f0 <datafile> -test_zeropivot -pc_factor_shift_type <shift_type> */ 157 flg = PETSC_FALSE; 158 ierr = PetscOptionsGetBool(NULL,NULL, "-test_zeropivot", &flg,NULL);CHKERRQ(ierr); 159 if (flg) { /* set a row as zeros */ 160 PetscInt row=0; 161 ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);CHKERRQ(ierr); 162 ierr = MatZeroRows(A,1,&row,0.0,NULL,NULL);CHKERRQ(ierr); 163 } 164 165 /* Check whether A is symmetric, then set A->symmetric option */ 166 flg = PETSC_FALSE; 167 ierr = PetscOptionsGetBool(NULL,NULL, "-check_symmetry", &flg,NULL);CHKERRQ(ierr); 168 if (flg) { 169 ierr = MatIsSymmetric(A,0.0,&isSymmetric);CHKERRQ(ierr); 170 if (!isSymmetric) { 171 ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");CHKERRQ(ierr); 172 } 173 } 174 175 /* 176 If the loaded matrix is larger than the vector (due to being padded 177 to match the block size of the system), then create a new padded vector. 178 */ 179 180 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 181 /* if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/ 182 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 183 ierr = VecGetSize(b,&m);CHKERRQ(ierr); 184 if (M != m) { /* Create a new vector b by padding the old one */ 185 PetscInt j,mvec,start,end,indx; 186 Vec tmp; 187 PetscScalar *bold; 188 189 ierr = VecCreate(PETSC_COMM_WORLD,&tmp);CHKERRQ(ierr); 190 ierr = VecSetSizes(tmp,n,PETSC_DECIDE);CHKERRQ(ierr); 191 ierr = VecSetFromOptions(tmp);CHKERRQ(ierr); 192 ierr = VecGetOwnershipRange(b,&start,&end);CHKERRQ(ierr); 193 ierr = VecGetLocalSize(b,&mvec);CHKERRQ(ierr); 194 ierr = VecGetArray(b,&bold);CHKERRQ(ierr); 195 for (j=0; j<mvec; j++) { 196 indx = start+j; 197 ierr = VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);CHKERRQ(ierr); 198 } 199 ierr = VecRestoreArray(b,&bold);CHKERRQ(ierr); 200 ierr = VecDestroy(&b);CHKERRQ(ierr); 201 ierr = VecAssemblyBegin(tmp);CHKERRQ(ierr); 202 ierr = VecAssemblyEnd(tmp);CHKERRQ(ierr); 203 b = tmp; 204 } 205 206 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 207 ierr = VecDuplicate(b,&u);CHKERRQ(ierr); 208 if (initialguessfile) { 209 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 210 ierr = VecLoad(x,viewer);CHKERRQ(ierr); 211 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 212 initialguess = PETSC_TRUE; 213 } else if (initialguess) { 214 ierr = VecSet(x,1.0);CHKERRQ(ierr); 215 } else { 216 ierr = VecSet(x,0.0);CHKERRQ(ierr); 217 } 218 219 /* Check scaling in A */ 220 flg = PETSC_FALSE; 221 ierr = PetscOptionsGetBool(NULL,NULL, "-check_scaling", &flg,NULL);CHKERRQ(ierr); 222 if (flg) { 223 Vec max, min; 224 PetscInt idx; 225 PetscReal val; 226 227 ierr = VecDuplicate(x, &max);CHKERRQ(ierr); 228 ierr = VecDuplicate(x, &min);CHKERRQ(ierr); 229 ierr = MatGetRowMaxAbs(A, max, NULL);CHKERRQ(ierr); 230 ierr = MatGetRowMinAbs(A, min, NULL);CHKERRQ(ierr); 231 { 232 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);CHKERRQ(ierr); 233 ierr = VecView(max, viewer);CHKERRQ(ierr); 234 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 235 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);CHKERRQ(ierr); 236 ierr = VecView(min, viewer);CHKERRQ(ierr); 237 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 238 } 239 ierr = VecView(max, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 240 ierr = VecMax(max, &idx, &val);CHKERRQ(ierr); 241 ierr = PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 242 ierr = VecView(min, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 243 ierr = VecMin(min, &idx, &val);CHKERRQ(ierr); 244 ierr = PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 245 ierr = VecMin(max, &idx, &val);CHKERRQ(ierr); 246 ierr = PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 247 ierr = VecPointwiseDivide(max, max, min);CHKERRQ(ierr); 248 ierr = VecMax(max, &idx, &val);CHKERRQ(ierr); 249 ierr = PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 250 ierr = VecView(max, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 251 ierr = VecDestroy(&max);CHKERRQ(ierr); 252 ierr = VecDestroy(&min);CHKERRQ(ierr); 253 } 254 255 /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 256 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 257 Setup solve for system 258 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 259 /* 260 Conclude profiling last stage; begin profiling next stage. 261 */ 262 PetscPreLoadStage("KSPSetUpSolve"); 263 264 /* 265 Create linear solver; set operators; set runtime options. 266 */ 267 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 268 ierr = KSPSetInitialGuessNonzero(ksp,initialguess);CHKERRQ(ierr); 269 num_numfac = 1; 270 ierr = PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);CHKERRQ(ierr); 271 while (num_numfac--) { 272 PC pc; 273 PetscBool lsqr,isbddc,ismatis; 274 char str[32]; 275 276 ierr = PetscOptionsGetString(NULL,NULL,"-ksp_type",str,sizeof(str),&lsqr);CHKERRQ(ierr); 277 if (lsqr) { 278 ierr = PetscStrcmp("lsqr",str,&lsqr);CHKERRQ(ierr); 279 } 280 if (lsqr) { 281 Mat BtB; 282 ierr = MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);CHKERRQ(ierr); 283 ierr = KSPSetOperators(ksp,A,BtB);CHKERRQ(ierr); 284 ierr = MatDestroy(&BtB);CHKERRQ(ierr); 285 } else { 286 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 287 } 288 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 289 290 /* if we test BDDC, make sure pmat is of type MATIS */ 291 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 292 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr); 293 ierr = PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);CHKERRQ(ierr); 294 if (isbddc && !ismatis) { 295 Mat J; 296 297 ierr = MatConvert(A,MATIS,MAT_INITIAL_MATRIX,&J);CHKERRQ(ierr); 298 ierr = KSPSetOperators(ksp,A,J);CHKERRQ(ierr); 299 ierr = MatDestroy(&J);CHKERRQ(ierr); 300 } 301 302 /* 303 Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to 304 enable more precise profiling of setting up the preconditioner. 305 These calls are optional, since both will be called within 306 KSPSolve() if they haven't been called already. 307 */ 308 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 309 ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr); 310 311 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 312 Solve system 313 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 314 315 /* 316 Solve linear system; 317 */ 318 if (trans) { 319 ierr = KSPSolveTranspose(ksp,b,x);CHKERRQ(ierr); 320 ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 321 } else { 322 PetscInt num_rhs=1; 323 ierr = PetscOptionsGetInt(NULL,NULL,"-num_rhs",&num_rhs,NULL);CHKERRQ(ierr); 324 cknorm = PETSC_FALSE; 325 ierr = PetscOptionsGetBool(NULL,NULL,"-cknorm",&cknorm,NULL);CHKERRQ(ierr); 326 while (num_rhs--) { 327 if (num_rhs == 1) VecSet(x,0.0); 328 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 329 } 330 ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 331 if (cknorm) { /* Check error for each rhs */ 332 if (trans) { 333 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); 334 } else { 335 ierr = MatMult(A,x,u);CHKERRQ(ierr); 336 } 337 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); 338 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 339 ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);CHKERRQ(ierr); 340 if (!PetscIsNanScalar(norm)) { 341 if (norm < 1.e-12) { 342 ierr = PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");CHKERRQ(ierr); 343 } else { 344 ierr = PetscPrintf(PETSC_COMM_WORLD," Residual norm %g\n",(double)norm);CHKERRQ(ierr); 345 } 346 } 347 } 348 } /* while (num_rhs--) */ 349 350 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 351 Check error, print output, free data structures. 352 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 353 354 /* 355 Check error 356 */ 357 if (trans) { 358 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); 359 } else { 360 ierr = MatMult(A,x,u);CHKERRQ(ierr); 361 } 362 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); 363 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 364 /* 365 Write output (optinally using table for solver details). 366 - PetscPrintf() handles output for multiprocessor jobs 367 by printing from only one processor in the communicator. 368 - KSPView() prints information about the linear solver. 369 */ 370 if (table) { 371 char *matrixname,kspinfo[120]; 372 373 /* 374 Open a string viewer; then write info to it. 375 */ 376 ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);CHKERRQ(ierr); 377 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 378 ierr = PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);CHKERRQ(ierr); 379 ierr = PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);CHKERRQ(ierr); 380 381 /* 382 Destroy the viewer 383 */ 384 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 385 } else { 386 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);CHKERRQ(ierr); 387 if (!PetscIsNanScalar(norm)) { 388 if (norm < 1.e-12 && !PetscIsNanScalar((PetscScalar)norm)) { 389 ierr = PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");CHKERRQ(ierr); 390 } else { 391 ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);CHKERRQ(ierr); 392 } 393 } 394 } 395 ierr = PetscOptionsGetString(NULL,NULL,"-solution",file[3],sizeof(file[3]),&flg);CHKERRQ(ierr); 396 if (flg) { 397 Vec xstar; 398 PetscReal norm; 399 400 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);CHKERRQ(ierr); 401 ierr = VecCreate(PETSC_COMM_WORLD,&xstar);CHKERRQ(ierr); 402 ierr = VecLoad(xstar,viewer);CHKERRQ(ierr); 403 ierr = VecAXPY(xstar, -1.0, x);CHKERRQ(ierr); 404 ierr = VecNorm(xstar, NORM_2, &norm);CHKERRQ(ierr); 405 ierr = PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)norm);CHKERRQ(ierr); 406 ierr = VecDestroy(&xstar);CHKERRQ(ierr); 407 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 408 } 409 if (outputSoln) { 410 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 411 ierr = VecView(x, viewer);CHKERRQ(ierr); 412 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 413 } 414 415 flg = PETSC_FALSE; 416 ierr = PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);CHKERRQ(ierr); 417 if (flg) { 418 KSPConvergedReason reason; 419 ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr); 420 ierr = PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);CHKERRQ(ierr); 421 } 422 423 } /* while (num_numfac--) */ 424 425 /* 426 Free work space. All PETSc objects should be destroyed when they 427 are no longer needed. 428 */ 429 ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); 430 ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); 431 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 432 PetscPreLoadEnd(); 433 /* ----------------------------------------------------------- 434 End of linear solver loop 435 ----------------------------------------------------------- */ 436 437 ierr = PetscFinalize(); 438 return ierr; 439 } 440 441 /*TEST 442 443 build: 444 requires: !complex 445 446 testset: 447 suffix: 1 448 nsize: 2 449 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int@PETSC_INDEX_SIZE@-float@PETSC_SCALAR_SIZE@ 450 requires: !__float128 451 452 testset: 453 suffix: 1a 454 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int@PETSC_INDEX_SIZE@-float@PETSC_SCALAR_SIZE@ 455 requires: !__float128 456 457 testset: 458 nsize: 2 459 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 460 args: -f0 ${DATAFILESPATH}/matrices/medium 461 args: -ksp_type bicg 462 test: 463 suffix: 2 464 465 testset: 466 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 467 args: -f0 ${DATAFILESPATH}/matrices/medium 468 args: -ksp_type bicg 469 test: 470 suffix: 4 471 args: -pc_type lu 472 test: 473 suffix: 5 474 475 testset: 476 suffix: 6 477 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 478 args: -f0 ${DATAFILESPATH}/matrices/fem1 479 args: -pc_factor_levels 2 -pc_factor_fill 1.73 -ksp_gmres_cgs_refinement_type refine_always 480 481 testset: 482 TODO: Matrix row/column sizes are not compatible with block size 483 suffix: 7 484 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 485 args: -f0 ${DATAFILESPATH}/matrices/medium 486 args: -viewer_binary_skip_info -mat_type seqbaij 487 args: -matload_block_size {{2 3 4 5 6 7 8}separate output} 488 args: -ksp_max_it 100 -ksp_gmres_cgs_refinement_type refine_always 489 args: -ksp_rtol 1.0e-15 -ksp_monitor_short 490 test: 491 suffix: a 492 test: 493 suffix: b 494 args: -pc_factor_mat_ordering_type nd 495 test: 496 suffix: c 497 args: -pc_factor_levels 1 498 499 testset: 500 TODO: Matrix row/column sizes are not compatible with block size 501 suffix: 7_d 502 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 503 args: -f0 ${DATAFILESPATH}/matrices/medium 504 args: -viewer_binary_skip_info -mat_type seqbaij 505 args: -matload_block_size {{2 3 4 5 6 7 8}shared output} 506 args: -ksp_type preonly -pc_type lu 507 508 testset: 509 suffix: 8 510 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 511 args: -f0 ${DATAFILESPATH}/matrices/medium 512 args: -ksp_diagonal_scale -pc_type eisenstat -ksp_monitor_short -ksp_diagonal_scale_fix -ksp_gmres_cgs_refinement_type refine_always -mat_no_inode 513 514 testset: 515 TODO: Matrix row/column sizes are not compatible with block size 516 suffix: 9 517 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 518 args: -f0 ${DATAFILESPATH}/matrices/medium 519 args: -viewer_binary_skip_info -matload_block_size {{1 2 3 4 5 6 7}separate output} -ksp_max_it 100 -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol 1.0e-15 -ksp_monitor_short 520 test: 521 suffix: a 522 args: -mat_type seqbaij 523 test: 524 suffix: b 525 args: -mat_type seqbaij -trans 526 test: 527 suffix: c 528 nsize: 2 529 args: -mat_type mpibaij 530 test: 531 suffix: d 532 nsize: 2 533 args: -mat_type mpibaij -trans 534 test: 535 suffix: e 536 nsize: 3 537 args: -mat_type mpibaij 538 test: 539 suffix: f 540 nsize: 3 541 args: -mat_type mpibaij -trans 542 543 testset: 544 suffix: 10 545 nsize: 2 546 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 547 args: -ksp_type fgmres -pc_type ksp -f0 ${DATAFILESPATH}/matrices/medium -ksp_fgmres_modifypcksp -ksp_monitor_short 548 549 testset: 550 suffix: 12 551 requires: matlab 552 args: -pc_type lu -pc_factor_mat_solver_type matlab -f0 ${DATAFILESPATH}/matrices/arco1 553 554 testset: 555 suffix: 13 556 requires: lusol 557 args: -f0 ${DATAFILESPATH}/matrices/arco1 558 args: -mat_type lusol -pc_type lu 559 560 testset: 561 nsize: 3 562 args: -f0 ${DATAFILESPATH}/matrices/medium 563 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 564 test: 565 suffix: 14 566 requires: spai 567 args: -pc_type spai 568 test: 569 suffix: 15 570 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 571 args: -pc_type hypre -pc_hypre_type pilut 572 test: 573 suffix: 16 574 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 575 args: -pc_type hypre -pc_hypre_type parasails 576 test: 577 suffix: 17 578 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 579 args: -pc_type hypre -pc_hypre_type boomeramg 580 test: 581 suffix: 18 582 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 583 args: -pc_type hypre -pc_hypre_type euclid 584 585 testset: 586 suffix: 19 587 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 588 args: -f0 ${DATAFILESPATH}/matrices/poisson1 589 args: -ksp_type cg -pc_type icc 590 args: -pc_factor_levels {{0 2 4}separate output} 591 test: 592 test: 593 args: -mat_type seqsbaij 594 595 testset: 596 suffix: ILU 597 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 598 args: -f0 ${DATAFILESPATH}/matrices/small 599 args: -pc_factor_levels 1 600 test: 601 test: 602 # This is tested against regular ILU (used to be denoted ILUBAIJ) 603 args: -mat_type baij 604 605 testset: 606 suffix: aijcusparse 607 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) cuda 608 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_monitor_short -ksp_view -mat_view ascii::ascii_info -mat_type aijcusparse -pc_factor_mat_solver_type cusparse -pc_type ilu -vec_type cuda 609 610 testset: 611 TODO: No output file. Need to determine if deprecated 612 suffix: asm_viennacl 613 nsize: 2 614 requires: viennacl 615 args: -pc_type asm -pc_asm_sub_mat_type aijviennacl -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int${PETSC_INDEX_SIZE}-float${PETSC_SCALAR_SIZE} 616 617 testset: 618 nsize: 2 619 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 620 args: -f0 ${DATAFILESPATH}/matrices/poisson2.gz -ksp_monitor_short -ksp_rtol 1.E-9 -pc_type hypre -pc_hypre_type boomeramg 621 test: 622 suffix: boomeramg_euclid 623 args: -pc_hypre_boomeramg_smooth_type Euclid -pc_hypre_boomeramg_smooth_num_levels 2 -pc_hypre_boomeramg_eu_level 1 -pc_hypre_boomeramg_eu_droptolerance 0.01 624 TODO: Need to determine if deprecated 625 test: 626 suffix: boomeramg_euclid_bj 627 args: -pc_hypre_boomeramg_smooth_type Euclid -pc_hypre_boomeramg_smooth_num_levels 2 -pc_hypre_boomeramg_eu_level 1 -pc_hypre_boomeramg_eu_droptolerance 0.01 -pc_hypre_boomeramg_eu_bj 628 TODO: Need to determine if deprecated 629 test: 630 suffix: boomeramg_parasails 631 args: -pc_hypre_boomeramg_smooth_type ParaSails -pc_hypre_boomeramg_smooth_num_levels 2 632 test: 633 suffix: boomeramg_pilut 634 args: -pc_hypre_boomeramg_smooth_type Pilut -pc_hypre_boomeramg_smooth_num_levels 2 635 test: 636 suffix: boomeramg_schwarz 637 args: -pc_hypre_boomeramg_smooth_type Schwarz-smoothers 638 639 testset: 640 suffix: cg_singlereduction 641 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 642 args: -f0 ${DATAFILESPATH}/matrices/small 643 args: -mat_type mpisbaij -ksp_type cg -pc_type eisenstat -ksp_monitor_short -ksp_converged_reason 644 test: 645 test: 646 args: -ksp_cg_single_reduction 647 648 testset: 649 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 650 args: -f0 ${DATAFILESPATH}/matrices/poisson2.gz 651 args: -ksp_monitor_short -pc_type icc 652 test: 653 suffix: cr 654 args: -ksp_type cr 655 test: 656 suffix: lcd 657 args: -ksp_type lcd 658 659 testset: 660 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 661 args: -f0 ${DATAFILESPATH}/matrices/small 662 args: -ksp_monitor_short -ksp_view -mat_view ascii::ascii_info 663 test: 664 suffix: seqaijcrl 665 args: -mat_type seqaijcrl 666 test: 667 suffix: seqaijperm 668 args: -mat_type seqaijperm 669 670 testset: 671 nsize: 2 672 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 673 args: -f0 ${DATAFILESPATH}/matrices/small 674 args: -ksp_monitor_short -ksp_view 675 # Different output files 676 test: 677 suffix: mpiaijcrl 678 args: -mat_type mpiaijcrl 679 test: 680 suffix: mpiaijperm 681 args: -mat_type mpiaijperm 682 683 testset: 684 nsize: 4 685 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 686 args: -ksp_monitor_short -ksp_view 687 test: 688 suffix: xxt 689 args: -f0 ${DATAFILESPATH}/matrices/poisson1 -check_symmetry -ksp_type cg -pc_type tfs 690 test: 691 suffix: xyt 692 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type gmres -pc_type tfs 693 694 testset: 695 # The output file here is the same as mumps 696 suffix: mumps_cholesky 697 output_file: output/ex72_mumps.out 698 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 699 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 700 nsize: {{1 2}} 701 test: 702 args: -mat_type sbaij -mat_ignore_lower_triangular 703 test: 704 args: -mat_type aij 705 test: 706 args: -mat_type aij -matload_spd 707 708 testset: 709 # The output file here is the same as mumps 710 suffix: mumps_lu 711 output_file: output/ex72_mumps.out 712 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 713 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 714 test: 715 args: -mat_type seqaij 716 test: 717 nsize: 2 718 args: -mat_type mpiaij 719 test: 720 args: -mat_type seqbaij -matload_block_size 2 721 test: 722 nsize: 2 723 args: -mat_type mpibaij -matload_block_size 2 724 test: 725 args: -mat_type aij -mat_mumps_icntl_7 5 726 TODO: Need to determine if deprecated 727 728 test: 729 suffix: mumps_lu_parmetis 730 output_file: output/ex72_mumps.out 731 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps parmetis 732 nsize: 2 733 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 -mat_type mpiaij -mat_mumps_icntl_28 2 -mat_mumps_icntl_29 2 734 735 test: 736 suffix: mumps_lu_ptscotch 737 output_file: output/ex72_mumps.out 738 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps ptscotch 739 nsize: 2 740 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 -mat_type mpiaij -mat_mumps_icntl_28 2 -mat_mumps_icntl_29 1 741 742 testset: 743 # The output file here is the same as mumps 744 suffix: mumps_redundant 745 output_file: output/ex72_mumps_redundant.out 746 nsize: 8 747 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 748 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type preonly -pc_type redundant -pc_redundant_number {{8 7 6 5 4 3 2 1}} -redundant_pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 749 750 testset: 751 suffix: pastix_cholesky 752 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) pastix 753 output_file: output/ex72_mumps.out 754 nsize: {{1 2}} 755 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_factor_mat_solver_type pastix -num_numfac 2 -num_rhs 2 -pc_type cholesky -mat_type sbaij -mat_ignore_lower_triangular 756 757 testset: 758 suffix: pastix_lu 759 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) pastix 760 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type pastix -num_numfac 2 -num_rhs 2 761 output_file: output/ex72_mumps.out 762 test: 763 args: -mat_type seqaij 764 test: 765 nsize: 2 766 args: -mat_type mpiaij 767 768 testset: 769 suffix: pastix_redundant 770 output_file: output/ex72_mumps_redundant.out 771 nsize: 8 772 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) pastix 773 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type preonly -pc_type redundant -pc_redundant_number {{8 7 6 5 4 3 2 1}} -redundant_pc_factor_mat_solver_type pastix -num_numfac 2 -num_rhs 2 774 775 testset: 776 suffix: superlu_dist_lu 777 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) superlu_dist 778 output_file: output/ex72_mumps.out 779 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu_dist -num_numfac 2 -num_rhs 2 780 nsize: {{1 2}} 781 782 testset: 783 suffix: superlu_dist_redundant 784 nsize: 8 785 output_file: output/ex72_mumps_redundant.out 786 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) superlu_dist 787 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type preonly -pc_type redundant -pc_redundant_number {{8 7 6 5 4 3 2 1}} -redundant_pc_factor_mat_solver_type superlu_dist -num_numfac 2 -num_rhs 2 788 789 testset: 790 suffix: superlu_lu 791 output_file: output/ex72_mumps.out 792 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) superlu 793 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu -num_numfac 2 -num_rhs 2 794 795 testset: 796 suffix: umfpack 797 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) suitesparse 798 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -mat_type seqaij -pc_factor_mat_solver_type umfpack -num_numfac 2 -num_rhs 2 799 800 testset: 801 suffix: zeropivot 802 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 803 args: -f0 ${DATAFILESPATH}/matrices/small -test_zeropivot -ksp_converged_reason -ksp_type fgmres -pc_type ksp 804 test: 805 nsize: 3 806 args: -ksp_pc_type bjacobi 807 test: 808 nsize: 2 809 args: -ksp_ksp_type cg -ksp_pc_type bjacobi -ksp_pc_bjacobi_blocks 1 810 #test: 811 #nsize: 3 812 #args: -ksp_ksp_converged_reason -ksp_pc_type bjacobi -ksp_sub_ksp_converged_reason 813 #TODO: Need to determine if deprecated 814 815 testset: 816 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 817 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type fgmres 818 test: 819 suffix: bddc_seq 820 nsize: 1 821 args: -pc_type bddc 822 test: 823 suffix: bddc_par 824 nsize: 2 825 args: -pc_type bddc 826 test: 827 requires: parmetis 828 suffix: bddc_par_nd_parmetis 829 filter: sed -e "s/Number of iterations = [0-9]/Number of iterations = 9/g" 830 nsize: 4 831 args: -ksp_error_if_not_converged -pc_type bddc -mat_is_disassemble_l2g_type nd -mat_partitioning_type parmetis 832 test: 833 requires: ptscotch defined(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND) 834 suffix: bddc_par_nd_ptscotch 835 filter: sed -e "s/Number of iterations = [0-9]/Number of iterations = 9/g" 836 nsize: 4 837 args: -ksp_error_if_not_converged -pc_type bddc -mat_is_disassemble_l2g_type nd -mat_partitioning_type ptscotch 838 839 test: 840 suffix: hpddm 841 output_file: output/ex72_bddc_seq.out 842 filter: sed -e "s/Number of iterations = 2/Number of iterations = 1/g" 843 nsize: 2 844 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int@PETSC_INDEX_SIZE@-float@PETSC_SCALAR_SIZE@ -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 5 -pc_hpddm_levels_1_st_pc_type mat 845 requires: !__float128 hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES) 846 TEST*/ 847