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);CHKERRQ(ierr); 27 28 /* Determine file from which we read the matrix A */ 29 ierr = PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 30 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 31 32 /* Load matrix A */ 33 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr); 34 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 35 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 36 ierr = MatLoad(A,fd);CHKERRQ(ierr); 37 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 38 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 39 if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%D, %D)", m, n); 40 41 /* if A is symmetric, set its flag -- required by MatGetInertia() */ 42 ierr = MatIsSymmetric(A,0.0,&flg);CHKERRQ(ierr); 43 44 ierr = MatViewFromOptions(A,NULL,"-A_view");CHKERRQ(ierr); 45 46 /* Create dense matrix C and X; C holds true solution with identical colums */ 47 nrhs = 2; 48 ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr); 49 ierr = PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %D\n",nrhs);CHKERRQ(ierr); 50 ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr); 51 ierr = MatSetOptionsPrefix(C,"rhs_");CHKERRQ(ierr); 52 ierr = MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);CHKERRQ(ierr); 53 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr); 54 ierr = MatSetFromOptions(C);CHKERRQ(ierr); 55 ierr = MatSetUp(C);CHKERRQ(ierr); 56 57 ierr = PetscOptionsGetBool(NULL,NULL,"-view_factor",&view,NULL);CHKERRQ(ierr); 58 ierr = PetscOptionsGetBool(NULL,NULL,"-test_matmatsolve",&testMatMatSolve,NULL);CHKERRQ(ierr); 59 ierr = PetscOptionsGetBool(NULL,NULL,"-cholesky",&chol,NULL);CHKERRQ(ierr); 60 #if defined(PETSC_HAVE_MUMPS) 61 ierr = PetscOptionsGetBool(NULL,NULL,"-test_mumps_opts",&test_mumps_opts,NULL);CHKERRQ(ierr); 62 #endif 63 64 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 65 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 66 ierr = MatSetRandom(C,rand);CHKERRQ(ierr); 67 ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);CHKERRQ(ierr); 68 69 /* Create vectors */ 70 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 71 ierr = VecSetSizes(x,n,PETSC_DECIDE);CHKERRQ(ierr); 72 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 73 ierr = VecDuplicate(x,&b);CHKERRQ(ierr); 74 ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* save the true solution */ 75 76 /* Test Factorization */ 77 ierr = MatGetOrdering(A,MATORDERINGND,&perm,&iperm);CHKERRQ(ierr); 78 79 ierr = PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);CHKERRQ(ierr); 80 switch (ipack) { 81 #if defined(PETSC_HAVE_SUPERLU) 82 case 0: 83 if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!"); 84 ierr = PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");CHKERRQ(ierr); 85 ierr = MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 86 matsolvexx = PETSC_TRUE; 87 break; 88 #endif 89 #if defined(PETSC_HAVE_SUPERLU_DIST) 90 case 1: 91 if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!"); 92 ierr = PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");CHKERRQ(ierr); 93 ierr = MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 94 matsolvexx = PETSC_TRUE; 95 break; 96 #endif 97 #if defined(PETSC_HAVE_MUMPS) 98 case 2: 99 if (chol) { 100 ierr = PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n");CHKERRQ(ierr); 101 ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 102 } else { 103 ierr = PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");CHKERRQ(ierr); 104 ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 105 } 106 matsolvexx = PETSC_TRUE; 107 if (test_mumps_opts) { 108 /* test mumps options */ 109 PetscInt icntl; 110 PetscReal cntl; 111 112 icntl = 2; /* sequential matrix ordering */ 113 ierr = MatMumpsSetIcntl(F,7,icntl);CHKERRQ(ierr); 114 115 cntl = 1.e-6; /* threshold for row pivot detection */ 116 ierr = MatMumpsSetIcntl(F,24,1);CHKERRQ(ierr); 117 ierr = MatMumpsSetCntl(F,3,cntl);CHKERRQ(ierr); 118 } 119 break; 120 #endif 121 #if defined(PETSC_HAVE_MKL_PARDISO) 122 case 3: 123 if (chol) { 124 ierr = PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n");CHKERRQ(ierr); 125 ierr = MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 126 } else { 127 ierr = PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");CHKERRQ(ierr); 128 ierr = MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 129 } 130 break; 131 #endif 132 default: 133 if (chol) { 134 ierr = PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n");CHKERRQ(ierr); 135 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 136 } else { 137 ierr = PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");CHKERRQ(ierr); 138 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 139 } 140 matsolvexx = PETSC_TRUE; 141 } 142 143 ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr); 144 info.fill = 5.0; 145 info.shifttype = (PetscReal) MAT_SHIFT_NONE; 146 if (chol) { 147 ierr = MatCholeskyFactorSymbolic(F,A,perm,&info);CHKERRQ(ierr); 148 } else { 149 ierr = MatLUFactorSymbolic(F,A,perm,iperm,&info);CHKERRQ(ierr); 150 } 151 152 for (nfact = 0; nfact < 2; nfact++) { 153 if (chol) { 154 ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the CHOLESKY numfactorization \n",nfact);CHKERRQ(ierr); 155 ierr = MatCholeskyFactorNumeric(F,A,&info);CHKERRQ(ierr); 156 } else { 157 ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the LU numfactorization \n",nfact);CHKERRQ(ierr); 158 ierr = MatLUFactorNumeric(F,A,&info);CHKERRQ(ierr); 159 } 160 if (view) { 161 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 162 ierr = MatView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 163 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 164 view = PETSC_FALSE; 165 } 166 167 #if defined(PETSC_HAVE_SUPERLU_DIST) 168 if (ipack == 1) { /* Test MatSuperluDistGetDiagU() 169 -- input: matrix factor F; output: main diagonal of matrix U on all processes */ 170 PetscInt M; 171 PetscScalar *diag; 172 #if !defined(PETSC_USE_COMPLEX) 173 PetscInt nneg,nzero,npos; 174 #endif 175 176 ierr = MatGetSize(F,&M,NULL);CHKERRQ(ierr); 177 ierr = PetscMalloc1(M,&diag);CHKERRQ(ierr); 178 ierr = MatSuperluDistGetDiagU(F,diag);CHKERRQ(ierr); 179 ierr = PetscFree(diag);CHKERRQ(ierr); 180 181 #if !defined(PETSC_USE_COMPLEX) 182 /* Test MatGetInertia() */ 183 ierr = MatGetInertia(F,&nneg,&nzero,&npos);CHKERRQ(ierr); 184 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);CHKERRQ(ierr); 185 #endif 186 } 187 #endif 188 189 /* Test MatMatSolve() */ 190 if (testMatMatSolve) { 191 if (!nfact) { 192 ierr = MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);CHKERRQ(ierr); 193 } else { 194 ierr = MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr); 195 } 196 for (nsolve = 0; nsolve < 2; nsolve++) { 197 ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the MatMatSolve \n",nsolve);CHKERRQ(ierr); 198 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr); 199 200 /* Check the error */ 201 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 202 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 203 if (norm > tol) { 204 ierr = PetscPrintf(PETSC_COMM_WORLD,"%D-the MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);CHKERRQ(ierr); 205 } 206 } 207 if (matsolvexx) { 208 /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */ 209 ierr = MatCopy(RHS,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 210 ierr = MatMatSolve(F,X,X);CHKERRQ(ierr); 211 /* Check the error */ 212 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 213 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 214 if (norm > tol) { 215 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm);CHKERRQ(ierr); 216 } 217 } 218 219 if (ipack == 2 && size == 1) { 220 Mat spRHS,spRHST,RHST; 221 222 ierr = MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);CHKERRQ(ierr); 223 ierr = MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);CHKERRQ(ierr); 224 ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr); 225 for (nsolve = 0; nsolve < 2; nsolve++) { 226 ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the sparse MatMatSolve \n",nsolve);CHKERRQ(ierr); 227 ierr = MatMatSolve(F,spRHS,X);CHKERRQ(ierr); 228 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,"%D-the sparse MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);CHKERRQ(ierr); 234 } 235 } 236 ierr = MatDestroy(&spRHST);CHKERRQ(ierr); 237 ierr = MatDestroy(&spRHS);CHKERRQ(ierr); 238 ierr = MatDestroy(&RHST);CHKERRQ(ierr); 239 } 240 } 241 242 /* Test MatSolve() */ 243 if (testMatSolve) { 244 for (nsolve = 0; nsolve < 2; nsolve++) { 245 ierr = VecSetRandom(x,rand);CHKERRQ(ierr); 246 ierr = VecCopy(x,u);CHKERRQ(ierr); 247 ierr = MatMult(A,x,b);CHKERRQ(ierr); 248 249 ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the MatSolve \n",nsolve);CHKERRQ(ierr); 250 ierr = MatSolve(F,b,x);CHKERRQ(ierr); 251 252 /* Check the error */ 253 ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr); /* u <- (-1.0)x + u */ 254 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 255 if (norm > tol) { 256 PetscReal resi; 257 ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */ 258 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); /* u <- (-1.0)b + u */ 259 ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr); 260 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %D\n",(double)norm,(double)resi,nfact);CHKERRQ(ierr); 261 } 262 } 263 } 264 } 265 266 /* Free data structures */ 267 ierr = MatDestroy(&A);CHKERRQ(ierr); 268 ierr = MatDestroy(&C);CHKERRQ(ierr); 269 ierr = MatDestroy(&F);CHKERRQ(ierr); 270 ierr = MatDestroy(&X);CHKERRQ(ierr); 271 if (testMatMatSolve) { 272 ierr = MatDestroy(&RHS);CHKERRQ(ierr); 273 } 274 275 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 276 ierr = ISDestroy(&perm);CHKERRQ(ierr); 277 ierr = ISDestroy(&iperm);CHKERRQ(ierr); 278 ierr = VecDestroy(&x);CHKERRQ(ierr); 279 ierr = VecDestroy(&b);CHKERRQ(ierr); 280 ierr = VecDestroy(&u);CHKERRQ(ierr); 281 ierr = PetscFinalize(); 282 return ierr; 283 } 284 285 286 /*TEST 287 288 test: 289 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 290 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10 291 output_file: output/ex125.out 292 293 test: 294 suffix: mkl_pardiso 295 requires: mkl_pardiso datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 296 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3 297 298 test: 299 suffix: mumps 300 requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 301 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2 302 output_file: output/ex125_mumps_seq.out 303 304 test: 305 suffix: mumps_2 306 nsize: 3 307 requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 308 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2 309 output_file: output/ex125_mumps_par.out 310 311 test: 312 suffix: superlu_dist 313 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist 314 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM 315 316 test: 317 suffix: superlu_dist_2 318 nsize: 3 319 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist 320 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM 321 output_file: output/ex125_superlu_dist.out 322 323 test: 324 suffix: superlu_dist_complex 325 nsize: 3 326 requires: datafilespath superlu_dist complex double !define(PETSC_USE_64BIT_INDICES) 327 args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1 328 output_file: output/ex125_superlu_dist_complex.out 329 330 TEST*/ 331