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