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