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