1 static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\ 2 Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A,RHS,C,F,X; 9 Vec u,x,b; 10 PetscMPIInt size; 11 PetscInt m,n,nfact,nsolve,nrhs,ipack=0; 12 PetscReal norm,tol=1.e-10; 13 IS perm,iperm; 14 MatFactorInfo info; 15 PetscRandom rand; 16 PetscBool flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE; 17 PetscBool chol=PETSC_FALSE,view=PETSC_FALSE,matsolvexx = PETSC_FALSE; 18 #if defined(PETSC_HAVE_MUMPS) 19 PetscBool test_mumps_opts=PETSC_FALSE; 20 #endif 21 PetscViewer fd; /* viewer */ 22 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 23 24 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 25 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 26 27 /* Determine file from which we read the matrix A */ 28 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg)); 29 if (flg) { /* Load matrix A */ 30 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 31 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 32 CHKERRQ(MatSetFromOptions(A)); 33 CHKERRQ(MatLoad(A,fd)); 34 CHKERRQ(PetscViewerDestroy(&fd)); 35 } else { 36 n = 13; 37 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 38 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 39 CHKERRQ(MatSetType(A,MATAIJ)); 40 CHKERRQ(MatSetFromOptions(A)); 41 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n)); 42 CHKERRQ(MatSetUp(A)); 43 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 44 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 45 CHKERRQ(MatShift(A,1.0)); 46 } 47 CHKERRQ(MatGetLocalSize(A,&m,&n)); 48 PetscCheckFalse(m != n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 49 50 /* if A is symmetric, set its flag -- required by MatGetInertia() */ 51 CHKERRQ(MatIsSymmetric(A,0.0,&flg)); 52 53 CHKERRQ(MatViewFromOptions(A,NULL,"-A_view")); 54 55 /* Create dense matrix C and X; C holds true solution with identical columns */ 56 nrhs = 2; 57 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL)); 58 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %" PetscInt_FMT "\n",nrhs)); 59 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&C)); 60 CHKERRQ(MatSetOptionsPrefix(C,"rhs_")); 61 CHKERRQ(MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs)); 62 CHKERRQ(MatSetType(C,MATDENSE)); 63 CHKERRQ(MatSetFromOptions(C)); 64 CHKERRQ(MatSetUp(C)); 65 66 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-view_factor",&view,NULL)); 67 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_matmatsolve",&testMatMatSolve,NULL)); 68 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-cholesky",&chol,NULL)); 69 #if defined(PETSC_HAVE_MUMPS) 70 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_mumps_opts",&test_mumps_opts,NULL)); 71 #endif 72 73 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 74 CHKERRQ(PetscRandomSetFromOptions(rand)); 75 CHKERRQ(MatSetRandom(C,rand)); 76 CHKERRQ(MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X)); 77 78 /* Create vectors */ 79 CHKERRQ(MatCreateVecs(A,&x,&b)); 80 CHKERRQ(VecDuplicate(x,&u)); /* save the true solution */ 81 82 /* Test Factorization */ 83 CHKERRQ(MatGetOrdering(A,MATORDERINGND,&perm,&iperm)); 84 85 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL)); 86 switch (ipack) { 87 #if defined(PETSC_HAVE_SUPERLU) 88 case 0: 89 PetscCheck(!chol,PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!"); 90 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n")); 91 CHKERRQ(MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F)); 92 matsolvexx = PETSC_TRUE; 93 break; 94 #endif 95 #if defined(PETSC_HAVE_SUPERLU_DIST) 96 case 1: 97 PetscCheck(!chol,PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!"); 98 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n")); 99 CHKERRQ(MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F)); 100 matsolvexx = PETSC_TRUE; 101 break; 102 #endif 103 #if defined(PETSC_HAVE_MUMPS) 104 case 2: 105 if (chol) { 106 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n")); 107 CHKERRQ(MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F)); 108 } else { 109 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n")); 110 CHKERRQ(MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F)); 111 } 112 matsolvexx = PETSC_TRUE; 113 if (test_mumps_opts) { 114 /* test mumps options */ 115 PetscInt icntl; 116 PetscReal cntl; 117 118 icntl = 2; /* sequential matrix ordering */ 119 CHKERRQ(MatMumpsSetIcntl(F,7,icntl)); 120 121 cntl = 1.e-6; /* threshold for row pivot detection */ 122 CHKERRQ(MatMumpsSetIcntl(F,24,1)); 123 CHKERRQ(MatMumpsSetCntl(F,3,cntl)); 124 } 125 break; 126 #endif 127 #if defined(PETSC_HAVE_MKL_PARDISO) 128 case 3: 129 if (chol) { 130 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n")); 131 CHKERRQ(MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F)); 132 } else { 133 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n")); 134 CHKERRQ(MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F)); 135 } 136 break; 137 #endif 138 #if defined(PETSC_HAVE_CUDA) 139 case 4: 140 if (chol) { 141 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," CUSPARSE CHOLESKY:\n")); 142 CHKERRQ(MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_CHOLESKY,&F)); 143 } else { 144 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," CUSPARSE LU:\n")); 145 CHKERRQ(MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_LU,&F)); 146 } 147 break; 148 #endif 149 default: 150 if (chol) { 151 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n")); 152 CHKERRQ(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F)); 153 } else { 154 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n")); 155 CHKERRQ(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F)); 156 } 157 matsolvexx = PETSC_TRUE; 158 } 159 160 CHKERRQ(MatFactorInfoInitialize(&info)); 161 info.fill = 5.0; 162 info.shifttype = (PetscReal) MAT_SHIFT_NONE; 163 if (chol) { 164 CHKERRQ(MatCholeskyFactorSymbolic(F,A,perm,&info)); 165 } else { 166 CHKERRQ(MatLUFactorSymbolic(F,A,perm,iperm,&info)); 167 } 168 169 for (nfact = 0; nfact < 2; nfact++) { 170 if (chol) { 171 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the CHOLESKY numfactorization \n",nfact)); 172 CHKERRQ(MatCholeskyFactorNumeric(F,A,&info)); 173 } else { 174 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the LU numfactorization \n",nfact)); 175 CHKERRQ(MatLUFactorNumeric(F,A,&info)); 176 } 177 if (view) { 178 CHKERRQ(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO)); 179 CHKERRQ(MatView(F,PETSC_VIEWER_STDOUT_WORLD)); 180 CHKERRQ(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 181 view = PETSC_FALSE; 182 } 183 184 #if defined(PETSC_HAVE_SUPERLU_DIST) 185 if (ipack == 1) { /* Test MatSuperluDistGetDiagU() 186 -- input: matrix factor F; output: main diagonal of matrix U on all processes */ 187 PetscInt M; 188 PetscScalar *diag; 189 #if !defined(PETSC_USE_COMPLEX) 190 PetscInt nneg,nzero,npos; 191 #endif 192 193 CHKERRQ(MatGetSize(F,&M,NULL)); 194 CHKERRQ(PetscMalloc1(M,&diag)); 195 CHKERRQ(MatSuperluDistGetDiagU(F,diag)); 196 CHKERRQ(PetscFree(diag)); 197 198 #if !defined(PETSC_USE_COMPLEX) 199 /* Test MatGetInertia() */ 200 CHKERRQ(MatGetInertia(F,&nneg,&nzero,&npos)); 201 CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos)); 202 #endif 203 } 204 #endif 205 206 #if defined(PETSC_HAVE_MUMPS) 207 /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */ 208 if (ipack == 2) { 209 if (chol) { 210 CHKERRQ(MatCholeskyFactorSymbolic(F,A,perm,&info)); 211 CHKERRQ(MatCholeskyFactorNumeric(F,A,&info)); 212 } else { 213 CHKERRQ(MatLUFactorSymbolic(F,A,perm,iperm,&info)); 214 CHKERRQ(MatLUFactorNumeric(F,A,&info)); 215 } 216 } 217 #endif 218 219 /* Test MatMatSolve() */ 220 if (testMatMatSolve) { 221 if (!nfact) { 222 CHKERRQ(MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS)); 223 } else { 224 CHKERRQ(MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS)); 225 } 226 for (nsolve = 0; nsolve < 2; nsolve++) { 227 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the MatMatSolve \n",nsolve)); 228 CHKERRQ(MatMatSolve(F,RHS,X)); 229 230 /* Check the error */ 231 CHKERRQ(MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN)); 232 CHKERRQ(MatNorm(X,NORM_FROBENIUS,&norm)); 233 if (norm > tol) { 234 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n",nsolve,(double)norm,nsolve)); 235 } 236 } 237 if (matsolvexx) { 238 /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */ 239 CHKERRQ(MatCopy(RHS,X,SAME_NONZERO_PATTERN)); 240 CHKERRQ(MatMatSolve(F,X,X)); 241 /* Check the error */ 242 CHKERRQ(MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN)); 243 CHKERRQ(MatNorm(X,NORM_FROBENIUS,&norm)); 244 if (norm > tol) { 245 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm)); 246 } 247 } 248 249 if (ipack == 2 && size == 1) { 250 Mat spRHS,spRHST,RHST; 251 252 CHKERRQ(MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST)); 253 CHKERRQ(MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST)); 254 CHKERRQ(MatCreateTranspose(spRHST,&spRHS)); 255 for (nsolve = 0; nsolve < 2; nsolve++) { 256 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the sparse MatMatSolve \n",nsolve)); 257 CHKERRQ(MatMatSolve(F,spRHS,X)); 258 259 /* Check the error */ 260 CHKERRQ(MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN)); 261 CHKERRQ(MatNorm(X,NORM_FROBENIUS,&norm)); 262 if (norm > tol) { 263 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n",nsolve,(double)norm,nsolve)); 264 } 265 } 266 CHKERRQ(MatDestroy(&spRHST)); 267 CHKERRQ(MatDestroy(&spRHS)); 268 CHKERRQ(MatDestroy(&RHST)); 269 } 270 } 271 272 /* Test MatSolve() */ 273 if (testMatSolve) { 274 for (nsolve = 0; nsolve < 2; nsolve++) { 275 CHKERRQ(VecSetRandom(x,rand)); 276 CHKERRQ(VecCopy(x,u)); 277 CHKERRQ(MatMult(A,x,b)); 278 279 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the MatSolve \n",nsolve)); 280 CHKERRQ(MatSolve(F,b,x)); 281 282 /* Check the error */ 283 CHKERRQ(VecAXPY(u,-1.0,x)); /* u <- (-1.0)x + u */ 284 CHKERRQ(VecNorm(u,NORM_2,&norm)); 285 if (norm > tol) { 286 PetscReal resi; 287 CHKERRQ(MatMult(A,x,u)); /* u = A*x */ 288 CHKERRQ(VecAXPY(u,-1.0,b)); /* u <- (-1.0)b + u */ 289 CHKERRQ(VecNorm(u,NORM_2,&resi)); 290 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n",(double)norm,(double)resi,nfact)); 291 } 292 } 293 } 294 } 295 296 /* Free data structures */ 297 CHKERRQ(MatDestroy(&A)); 298 CHKERRQ(MatDestroy(&C)); 299 CHKERRQ(MatDestroy(&F)); 300 CHKERRQ(MatDestroy(&X)); 301 if (testMatMatSolve) { 302 CHKERRQ(MatDestroy(&RHS)); 303 } 304 305 CHKERRQ(PetscRandomDestroy(&rand)); 306 CHKERRQ(ISDestroy(&perm)); 307 CHKERRQ(ISDestroy(&iperm)); 308 CHKERRQ(VecDestroy(&x)); 309 CHKERRQ(VecDestroy(&b)); 310 CHKERRQ(VecDestroy(&u)); 311 CHKERRQ(PetscFinalize()); 312 return 0; 313 } 314 315 /*TEST 316 317 test: 318 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 319 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10 320 output_file: output/ex125.out 321 322 test: 323 suffix: 2 324 args: -mat_solver_type 10 325 output_file: output/ex125.out 326 327 test: 328 suffix: mkl_pardiso 329 requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 330 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3 331 332 test: 333 suffix: mkl_pardiso_2 334 requires: mkl_pardiso 335 args: -mat_solver_type 3 336 output_file: output/ex125_mkl_pardiso.out 337 338 test: 339 suffix: mumps 340 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 341 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2 342 output_file: output/ex125_mumps_seq.out 343 344 test: 345 suffix: mumps_2 346 nsize: 3 347 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 348 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2 349 output_file: output/ex125_mumps_par.out 350 351 test: 352 suffix: mumps_3 353 requires: mumps 354 args: -mat_solver_type 2 355 output_file: output/ex125_mumps_seq.out 356 357 test: 358 suffix: mumps_4 359 nsize: 3 360 requires: mumps 361 args: -mat_solver_type 2 362 output_file: output/ex125_mumps_par.out 363 364 test: 365 suffix: mumps_5 366 nsize: 3 367 requires: mumps 368 args: -mat_solver_type 2 -cholesky 369 output_file: output/ex125_mumps_par_cholesky.out 370 371 test: 372 suffix: superlu_dist 373 nsize: {{1 3}} 374 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) superlu_dist 375 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM 376 377 test: 378 suffix: superlu_dist_2 379 nsize: {{1 3}} 380 requires: superlu_dist !complex 381 args: -n 36 -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM 382 output_file: output/ex125_superlu_dist.out 383 384 test: 385 suffix: superlu_dist_complex 386 nsize: 3 387 requires: datafilespath superlu_dist complex double !defined(PETSC_USE_64BIT_INDICES) 388 args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1 389 output_file: output/ex125_superlu_dist_complex.out 390 391 test: 392 suffix: superlu_dist_complex_2 393 nsize: 3 394 requires: superlu_dist complex 395 args: -mat_solver_type 1 396 output_file: output/ex125_superlu_dist_complex.out 397 398 test: 399 suffix: cusparse 400 requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 401 args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output} 402 403 test: 404 suffix: cusparse_2 405 requires: cuda 406 args: -mat_type aijcusparse -mat_solver_type 4 -cholesky {{0 1}separate output} 407 408 TEST*/ 409