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,p,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,NULL,&n);CHKERRQ(ierr); 181 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 182 ierr = VecGetSize(b,&m);CHKERRQ(ierr); 183 ierr = VecGetLocalSize(b,&p);CHKERRQ(ierr); 184 preload = (PetscBool)(M != m || p != n); /* Global or local dimension mismatch */ 185 ierr = MPIU_Allreduce(&preload,&flg,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRMPI(ierr); 186 if (flg) { /* Create a new vector b by padding the old one */ 187 PetscInt j,mvec,start,end,indx; 188 Vec tmp; 189 PetscScalar *bold; 190 191 ierr = VecCreate(PETSC_COMM_WORLD,&tmp);CHKERRQ(ierr); 192 ierr = VecSetSizes(tmp,n,PETSC_DECIDE);CHKERRQ(ierr); 193 ierr = VecSetFromOptions(tmp);CHKERRQ(ierr); 194 ierr = VecGetOwnershipRange(b,&start,&end);CHKERRQ(ierr); 195 ierr = VecGetLocalSize(b,&mvec);CHKERRQ(ierr); 196 ierr = VecGetArray(b,&bold);CHKERRQ(ierr); 197 for (j=0; j<mvec; j++) { 198 indx = start+j; 199 ierr = VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);CHKERRQ(ierr); 200 } 201 ierr = VecRestoreArray(b,&bold);CHKERRQ(ierr); 202 ierr = VecDestroy(&b);CHKERRQ(ierr); 203 ierr = VecAssemblyBegin(tmp);CHKERRQ(ierr); 204 ierr = VecAssemblyEnd(tmp);CHKERRQ(ierr); 205 b = tmp; 206 } 207 208 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 209 ierr = VecDuplicate(b,&u);CHKERRQ(ierr); 210 if (initialguessfile) { 211 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 212 ierr = VecLoad(x,viewer);CHKERRQ(ierr); 213 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 214 initialguess = PETSC_TRUE; 215 } else if (initialguess) { 216 ierr = VecSet(x,1.0);CHKERRQ(ierr); 217 } else { 218 ierr = VecSet(x,0.0);CHKERRQ(ierr); 219 } 220 221 /* Check scaling in A */ 222 flg = PETSC_FALSE; 223 ierr = PetscOptionsGetBool(NULL,NULL, "-check_scaling", &flg,NULL);CHKERRQ(ierr); 224 if (flg) { 225 Vec max, min; 226 PetscInt idx; 227 PetscReal val; 228 229 ierr = VecDuplicate(x, &max);CHKERRQ(ierr); 230 ierr = VecDuplicate(x, &min);CHKERRQ(ierr); 231 ierr = MatGetRowMaxAbs(A, max, NULL);CHKERRQ(ierr); 232 ierr = MatGetRowMinAbs(A, min, NULL);CHKERRQ(ierr); 233 { 234 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);CHKERRQ(ierr); 235 ierr = VecView(max, viewer);CHKERRQ(ierr); 236 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 237 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);CHKERRQ(ierr); 238 ierr = VecView(min, viewer);CHKERRQ(ierr); 239 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 240 } 241 ierr = VecView(max, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 242 ierr = VecMax(max, &idx, &val);CHKERRQ(ierr); 243 ierr = PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 244 ierr = VecView(min, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 245 ierr = VecMin(min, &idx, &val);CHKERRQ(ierr); 246 ierr = PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 247 ierr = VecMin(max, &idx, &val);CHKERRQ(ierr); 248 ierr = PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 249 ierr = VecPointwiseDivide(max, max, min);CHKERRQ(ierr); 250 ierr = VecMax(max, &idx, &val);CHKERRQ(ierr); 251 ierr = PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %g at row %D\n", (double)val, idx);CHKERRQ(ierr); 252 ierr = VecView(max, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 253 ierr = VecDestroy(&max);CHKERRQ(ierr); 254 ierr = VecDestroy(&min);CHKERRQ(ierr); 255 } 256 257 /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 258 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 259 Setup solve for system 260 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 261 /* 262 Conclude profiling last stage; begin profiling next stage. 263 */ 264 PetscPreLoadStage("KSPSetUpSolve"); 265 266 /* 267 Create linear solver; set operators; set runtime options. 268 */ 269 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 270 ierr = KSPSetInitialGuessNonzero(ksp,initialguess);CHKERRQ(ierr); 271 num_numfac = 1; 272 ierr = PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);CHKERRQ(ierr); 273 while (num_numfac--) { 274 PC pc; 275 PetscBool lsqr,isbddc,ismatis; 276 char str[32]; 277 278 ierr = PetscOptionsGetString(NULL,NULL,"-ksp_type",str,sizeof(str),&lsqr);CHKERRQ(ierr); 279 if (lsqr) { 280 ierr = PetscStrcmp("lsqr",str,&lsqr);CHKERRQ(ierr); 281 } 282 if (lsqr) { 283 Mat BtB; 284 ierr = MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);CHKERRQ(ierr); 285 ierr = KSPSetOperators(ksp,A,BtB);CHKERRQ(ierr); 286 ierr = MatDestroy(&BtB);CHKERRQ(ierr); 287 } else { 288 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 289 } 290 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 291 292 /* if we test BDDC, make sure pmat is of type MATIS */ 293 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 294 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr); 295 ierr = PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);CHKERRQ(ierr); 296 if (isbddc && !ismatis) { 297 Mat J; 298 299 ierr = MatConvert(A,MATIS,MAT_INITIAL_MATRIX,&J);CHKERRQ(ierr); 300 ierr = KSPSetOperators(ksp,A,J);CHKERRQ(ierr); 301 ierr = MatDestroy(&J);CHKERRQ(ierr); 302 } 303 304 /* 305 Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to 306 enable more precise profiling of setting up the preconditioner. 307 These calls are optional, since both will be called within 308 KSPSolve() if they haven't been called already. 309 */ 310 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 311 ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr); 312 313 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 314 Solve system 315 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 316 317 /* 318 Solve linear system; 319 */ 320 if (trans) { 321 ierr = KSPSolveTranspose(ksp,b,x);CHKERRQ(ierr); 322 ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 323 } else { 324 PetscInt num_rhs=1; 325 ierr = PetscOptionsGetInt(NULL,NULL,"-num_rhs",&num_rhs,NULL);CHKERRQ(ierr); 326 cknorm = PETSC_FALSE; 327 ierr = PetscOptionsGetBool(NULL,NULL,"-cknorm",&cknorm,NULL);CHKERRQ(ierr); 328 while (num_rhs--) { 329 if (num_rhs == 1) VecSet(x,0.0); 330 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 331 } 332 ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 333 if (cknorm) { /* Check error for each rhs */ 334 if (trans) { 335 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); 336 } else { 337 ierr = MatMult(A,x,u);CHKERRQ(ierr); 338 } 339 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); 340 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 341 ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);CHKERRQ(ierr); 342 if (!PetscIsNanScalar(norm)) { 343 if (norm < 1.e-12) { 344 ierr = PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");CHKERRQ(ierr); 345 } else { 346 ierr = PetscPrintf(PETSC_COMM_WORLD," Residual norm %g\n",(double)norm);CHKERRQ(ierr); 347 } 348 } 349 } 350 } /* while (num_rhs--) */ 351 352 /* - - - - - - - - - - - New Stage - - - - - - - - - - - - - 353 Check error, print output, free data structures. 354 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 355 356 /* 357 Check error 358 */ 359 if (trans) { 360 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); 361 } else { 362 ierr = MatMult(A,x,u);CHKERRQ(ierr); 363 } 364 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); 365 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 366 /* 367 Write output (optinally using table for solver details). 368 - PetscPrintf() handles output for multiprocessor jobs 369 by printing from only one processor in the communicator. 370 - KSPView() prints information about the linear solver. 371 */ 372 if (table) { 373 char *matrixname,kspinfo[120]; 374 375 /* 376 Open a string viewer; then write info to it. 377 */ 378 ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);CHKERRQ(ierr); 379 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 380 ierr = PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);CHKERRQ(ierr); 381 ierr = PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);CHKERRQ(ierr); 382 383 /* 384 Destroy the viewer 385 */ 386 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 387 } else { 388 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);CHKERRQ(ierr); 389 if (!PetscIsNanScalar(norm)) { 390 if (norm < 1.e-12 && !PetscIsNanScalar((PetscScalar)norm)) { 391 ierr = PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");CHKERRQ(ierr); 392 } else { 393 ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);CHKERRQ(ierr); 394 } 395 } 396 } 397 ierr = PetscOptionsGetString(NULL,NULL,"-solution",file[3],sizeof(file[3]),&flg);CHKERRQ(ierr); 398 if (flg) { 399 Vec xstar; 400 PetscReal norm; 401 402 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);CHKERRQ(ierr); 403 ierr = VecCreate(PETSC_COMM_WORLD,&xstar);CHKERRQ(ierr); 404 ierr = VecLoad(xstar,viewer);CHKERRQ(ierr); 405 ierr = VecAXPY(xstar, -1.0, x);CHKERRQ(ierr); 406 ierr = VecNorm(xstar, NORM_2, &norm);CHKERRQ(ierr); 407 ierr = PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)norm);CHKERRQ(ierr); 408 ierr = VecDestroy(&xstar);CHKERRQ(ierr); 409 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 410 } 411 if (outputSoln) { 412 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 413 ierr = VecView(x, viewer);CHKERRQ(ierr); 414 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 415 } 416 417 flg = PETSC_FALSE; 418 ierr = PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);CHKERRQ(ierr); 419 if (flg) { 420 KSPConvergedReason reason; 421 ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr); 422 ierr = PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);CHKERRQ(ierr); 423 } 424 425 } /* while (num_numfac--) */ 426 427 /* 428 Free work space. All PETSc objects should be destroyed when they 429 are no longer needed. 430 */ 431 ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); 432 ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); 433 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 434 PetscPreLoadEnd(); 435 /* ----------------------------------------------------------- 436 End of linear solver loop 437 ----------------------------------------------------------- */ 438 439 ierr = PetscFinalize(); 440 return ierr; 441 } 442 443 /*TEST 444 445 build: 446 requires: !complex 447 448 testset: 449 suffix: 1 450 nsize: 2 451 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int@PETSC_INDEX_SIZE@-float@PETSC_SCALAR_SIZE@ 452 requires: !__float128 453 454 testset: 455 suffix: 1a 456 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int@PETSC_INDEX_SIZE@-float@PETSC_SCALAR_SIZE@ 457 requires: !__float128 458 459 testset: 460 nsize: 2 461 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 462 args: -f0 ${DATAFILESPATH}/matrices/medium 463 args: -ksp_type bicg 464 test: 465 suffix: 2 466 467 testset: 468 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 469 args: -f0 ${DATAFILESPATH}/matrices/medium 470 args: -ksp_type bicg 471 test: 472 suffix: 4 473 args: -pc_type lu 474 test: 475 suffix: 5 476 477 testset: 478 suffix: 6 479 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 480 args: -f0 ${DATAFILESPATH}/matrices/fem1 481 args: -pc_factor_levels 2 -pc_factor_fill 1.73 -ksp_gmres_cgs_refinement_type refine_always 482 483 testset: 484 TODO: Matrix row/column sizes are not compatible with block size 485 suffix: 7 486 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 487 args: -f0 ${DATAFILESPATH}/matrices/medium 488 args: -viewer_binary_skip_info -mat_type seqbaij 489 args: -matload_block_size {{2 3 4 5 6 7 8}separate output} 490 args: -ksp_max_it 100 -ksp_gmres_cgs_refinement_type refine_always 491 args: -ksp_rtol 1.0e-15 -ksp_monitor_short 492 test: 493 suffix: a 494 test: 495 suffix: b 496 args: -pc_factor_mat_ordering_type nd 497 test: 498 suffix: c 499 args: -pc_factor_levels 1 500 501 testset: 502 TODO: Matrix row/column sizes are not compatible with block size 503 suffix: 7_d 504 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 505 args: -f0 ${DATAFILESPATH}/matrices/medium 506 args: -viewer_binary_skip_info -mat_type seqbaij 507 args: -matload_block_size {{2 3 4 5 6 7 8}shared output} 508 args: -ksp_type preonly -pc_type lu 509 510 testset: 511 suffix: 8 512 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 513 args: -f0 ${DATAFILESPATH}/matrices/medium 514 args: -ksp_diagonal_scale -pc_type eisenstat -ksp_monitor_short -ksp_diagonal_scale_fix -ksp_gmres_cgs_refinement_type refine_always -mat_no_inode 515 516 testset: 517 TODO: Matrix row/column sizes are not compatible with block size 518 suffix: 9 519 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 520 args: -f0 ${DATAFILESPATH}/matrices/medium 521 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 522 test: 523 suffix: a 524 args: -mat_type seqbaij 525 test: 526 suffix: b 527 args: -mat_type seqbaij -trans 528 test: 529 suffix: c 530 nsize: 2 531 args: -mat_type mpibaij 532 test: 533 suffix: d 534 nsize: 2 535 args: -mat_type mpibaij -trans 536 test: 537 suffix: e 538 nsize: 3 539 args: -mat_type mpibaij 540 test: 541 suffix: f 542 nsize: 3 543 args: -mat_type mpibaij -trans 544 545 testset: 546 suffix: 10 547 nsize: 2 548 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 549 args: -ksp_type fgmres -pc_type ksp -f0 ${DATAFILESPATH}/matrices/medium -ksp_fgmres_modifypcksp -ksp_monitor_short 550 551 testset: 552 suffix: 12 553 requires: matlab 554 args: -pc_type lu -pc_factor_mat_solver_type matlab -f0 ${DATAFILESPATH}/matrices/arco1 555 556 testset: 557 suffix: 13 558 requires: lusol 559 args: -f0 ${DATAFILESPATH}/matrices/arco1 560 args: -mat_type lusol -pc_type lu 561 562 testset: 563 nsize: 3 564 args: -f0 ${DATAFILESPATH}/matrices/medium 565 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 566 test: 567 suffix: 14 568 requires: spai 569 args: -pc_type spai 570 test: 571 suffix: 15 572 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 573 args: -pc_type hypre -pc_hypre_type pilut 574 test: 575 suffix: 16 576 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 577 args: -pc_type hypre -pc_hypre_type parasails 578 test: 579 suffix: 17 580 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 581 args: -pc_type hypre -pc_hypre_type boomeramg 582 test: 583 suffix: 18 584 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 585 args: -pc_type hypre -pc_hypre_type euclid 586 587 testset: 588 suffix: 19 589 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 590 args: -f0 ${DATAFILESPATH}/matrices/poisson1 591 args: -ksp_type cg -pc_type icc 592 args: -pc_factor_levels {{0 2 4}separate output} 593 test: 594 test: 595 args: -mat_type seqsbaij 596 597 testset: 598 suffix: ILU 599 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 600 args: -f0 ${DATAFILESPATH}/matrices/small 601 args: -pc_factor_levels 1 602 test: 603 test: 604 # This is tested against regular ILU (used to be denoted ILUBAIJ) 605 args: -mat_type baij 606 607 testset: 608 suffix: aijcusparse 609 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) cuda 610 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 611 612 testset: 613 TODO: No output file. Need to determine if deprecated 614 suffix: asm_viennacl 615 nsize: 2 616 requires: viennacl 617 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} 618 619 testset: 620 nsize: 2 621 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 622 args: -f0 ${DATAFILESPATH}/matrices/poisson2.gz -ksp_monitor_short -ksp_rtol 1.E-9 -pc_type hypre -pc_hypre_type boomeramg 623 test: 624 suffix: boomeramg_euclid 625 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 626 TODO: Need to determine if deprecated 627 test: 628 suffix: boomeramg_euclid_bj 629 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 630 TODO: Need to determine if deprecated 631 test: 632 suffix: boomeramg_parasails 633 args: -pc_hypre_boomeramg_smooth_type ParaSails -pc_hypre_boomeramg_smooth_num_levels 2 634 test: 635 suffix: boomeramg_pilut 636 args: -pc_hypre_boomeramg_smooth_type Pilut -pc_hypre_boomeramg_smooth_num_levels 2 637 test: 638 suffix: boomeramg_schwarz 639 args: -pc_hypre_boomeramg_smooth_type Schwarz-smoothers 640 641 testset: 642 suffix: cg_singlereduction 643 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 644 args: -f0 ${DATAFILESPATH}/matrices/small 645 args: -mat_type mpisbaij -ksp_type cg -pc_type eisenstat -ksp_monitor_short -ksp_converged_reason 646 test: 647 test: 648 args: -ksp_cg_single_reduction 649 650 testset: 651 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 652 args: -f0 ${DATAFILESPATH}/matrices/poisson2.gz 653 args: -ksp_monitor_short -pc_type icc 654 test: 655 suffix: cr 656 args: -ksp_type cr 657 test: 658 suffix: lcd 659 args: -ksp_type lcd 660 661 testset: 662 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 663 args: -f0 ${DATAFILESPATH}/matrices/small 664 args: -ksp_monitor_short -ksp_view -mat_view ascii::ascii_info 665 test: 666 suffix: seqaijcrl 667 args: -mat_type seqaijcrl 668 test: 669 suffix: seqaijperm 670 args: -mat_type seqaijperm 671 672 testset: 673 nsize: 2 674 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 675 args: -f0 ${DATAFILESPATH}/matrices/small 676 args: -ksp_monitor_short -ksp_view 677 # Different output files 678 test: 679 suffix: mpiaijcrl 680 args: -mat_type mpiaijcrl 681 test: 682 suffix: mpiaijperm 683 args: -mat_type mpiaijperm 684 685 testset: 686 nsize: 4 687 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !defined(PETSC_HAVE_I_MPI_NUMVERSION) 688 args: -ksp_monitor_short -ksp_view 689 test: 690 suffix: xxt 691 args: -f0 ${DATAFILESPATH}/matrices/poisson1 -check_symmetry -ksp_type cg -pc_type tfs 692 test: 693 suffix: xyt 694 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type gmres -pc_type tfs 695 696 testset: 697 # The output file here is the same as mumps 698 suffix: mumps_cholesky 699 output_file: output/ex72_mumps.out 700 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 701 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 702 nsize: {{1 2}} 703 test: 704 args: -mat_type sbaij -mat_ignore_lower_triangular 705 test: 706 args: -mat_type aij 707 test: 708 args: -mat_type aij -matload_spd 709 710 testset: 711 # The output file here is the same as mumps 712 suffix: mumps_lu 713 output_file: output/ex72_mumps.out 714 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 715 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -num_numfac 2 -num_rhs 2 716 test: 717 args: -mat_type seqaij 718 test: 719 nsize: 2 720 args: -mat_type mpiaij 721 test: 722 args: -mat_type seqbaij -matload_block_size 2 723 test: 724 nsize: 2 725 args: -mat_type mpibaij -matload_block_size 2 726 test: 727 args: -mat_type aij -mat_mumps_icntl_7 5 728 TODO: Need to determine if deprecated 729 730 test: 731 suffix: mumps_lu_parmetis 732 output_file: output/ex72_mumps.out 733 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps parmetis 734 nsize: 2 735 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 736 737 test: 738 suffix: mumps_lu_ptscotch 739 output_file: output/ex72_mumps.out 740 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps ptscotch 741 nsize: 2 742 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 743 744 testset: 745 # The output file here is the same as mumps 746 suffix: mumps_redundant 747 output_file: output/ex72_mumps_redundant.out 748 nsize: 8 749 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 750 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 751 752 testset: 753 suffix: pastix_cholesky 754 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) pastix 755 output_file: output/ex72_mumps.out 756 nsize: {{1 2}} 757 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 758 759 testset: 760 suffix: pastix_lu 761 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) pastix 762 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type pastix -num_numfac 2 -num_rhs 2 763 output_file: output/ex72_mumps.out 764 test: 765 args: -mat_type seqaij 766 test: 767 nsize: 2 768 args: -mat_type mpiaij 769 770 testset: 771 suffix: pastix_redundant 772 output_file: output/ex72_mumps_redundant.out 773 nsize: 8 774 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) pastix 775 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 776 777 testset: 778 suffix: superlu_dist_lu 779 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) superlu_dist 780 output_file: output/ex72_mumps.out 781 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu_dist -num_numfac 2 -num_rhs 2 782 nsize: {{1 2}} 783 784 testset: 785 suffix: superlu_dist_redundant 786 nsize: 8 787 output_file: output/ex72_mumps_redundant.out 788 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) superlu_dist 789 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 790 791 testset: 792 suffix: superlu_lu 793 output_file: output/ex72_mumps.out 794 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) superlu 795 args: -f0 ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu -num_numfac 2 -num_rhs 2 796 797 testset: 798 suffix: umfpack 799 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) suitesparse 800 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 801 802 testset: 803 suffix: zeropivot 804 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps 805 args: -f0 ${DATAFILESPATH}/matrices/small -test_zeropivot -ksp_converged_reason -ksp_type fgmres -pc_type ksp 806 test: 807 nsize: 3 808 args: -ksp_pc_type bjacobi 809 test: 810 nsize: 2 811 args: -ksp_ksp_type cg -ksp_pc_type bjacobi -ksp_pc_bjacobi_blocks 1 812 #test: 813 #nsize: 3 814 #args: -ksp_ksp_converged_reason -ksp_pc_type bjacobi -ksp_sub_ksp_converged_reason 815 #TODO: Need to determine if deprecated 816 817 testset: 818 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 819 args: -f0 ${DATAFILESPATH}/matrices/medium -ksp_type fgmres 820 test: 821 suffix: bddc_seq 822 nsize: 1 823 args: -pc_type bddc 824 test: 825 suffix: bddc_par 826 nsize: 2 827 args: -pc_type bddc 828 test: 829 requires: parmetis 830 suffix: bddc_par_nd_parmetis 831 filter: sed -e "s/Number of iterations = [0-9]/Number of iterations = 9/g" 832 nsize: 4 833 args: -ksp_error_if_not_converged -pc_type bddc -mat_is_disassemble_l2g_type nd -mat_partitioning_type parmetis 834 test: 835 requires: ptscotch defined(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND) 836 suffix: bddc_par_nd_ptscotch 837 filter: sed -e "s/Number of iterations = [0-9]/Number of iterations = 9/g" 838 nsize: 4 839 args: -ksp_error_if_not_converged -pc_type bddc -mat_is_disassemble_l2g_type nd -mat_partitioning_type ptscotch 840 841 testset: 842 requires: !__float128 hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES) 843 test: 844 suffix: hpddm_mat 845 output_file: output/ex72_bddc_seq.out 846 filter: sed -e "s/Number of iterations = 2/Number of iterations = 1/g" 847 nsize: 2 848 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 849 test: 850 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) 851 suffix: hpddm_gen_non_hermitian 852 output_file: output/ex72_2.out 853 nsize: 4 854 args: -f0 ${DATAFILESPATH}/matrices/arco1 -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 10 -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_eps_gen_non_hermitian -pc_hpddm_coarse_mat_type baij -pc_hpddm_block_splitting -pc_hpddm_levels_1_eps_threshold 0.7 -pc_hpddm_coarse_pc_type lu -ksp_pc_side right 855 test: 856 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) mumps !defined(PETSCTEST_VALGRIND) 857 suffix: hpddm_gen_non_hermitian_baij 858 output_file: output/ex72_19_pc_factor_levels-0.out 859 nsize: 4 860 timeoutfactor: 2 861 args: -f0 ${DATAFILESPATH}/matrices/arco6 -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 30 -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_eps_gen_non_hermitian -pc_hpddm_coarse_mat_type baij -pc_hpddm_block_splitting -pc_hpddm_levels_1_eps_threshold 0.8 -pc_hpddm_coarse_pc_type lu -ksp_pc_side right -mat_type baij -pc_hpddm_levels_1_sub_pc_factor_mat_solver_type mumps -pc_hpddm_levels_1_eps_tol 1.0e-2 862 TEST*/ 863