1 static char help[] = "Tests MatSolve() and MatMatSolve() with MUMPS or MKL_PARDISO sequential solvers in Schur complement mode.\n\ 2 Example: mpiexec -n 1 ./ex192 -f <matrix binary file> -nrhs 4 -symmetric_solve -hermitian_solve -schur_ratio 0.3\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A,RHS,C,F,X,S; 9 Vec u,x,b; 10 Vec xschur,bschur,uschur; 11 IS is_schur; 12 PetscErrorCode ierr; 13 PetscMPIInt size; 14 PetscInt isolver=0,size_schur,m,n,nfact,nsolve,nrhs; 15 PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON; 16 PetscRandom rand; 17 PetscBool data_provided,herm,symm,use_lu,cuda = PETSC_FALSE; 18 PetscReal sratio = 5.1/12.; 19 PetscViewer fd; /* viewer */ 20 char solver[256]; 21 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 22 23 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 24 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 25 PetscCheckFalse(size > 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor test"); 26 /* Determine which type of solver we want to test for */ 27 herm = PETSC_FALSE; 28 symm = PETSC_FALSE; 29 ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);CHKERRQ(ierr); 31 if (herm) symm = PETSC_TRUE; 32 ierr = PetscOptionsGetBool(NULL,NULL,"-cuda_solve",&cuda,NULL);CHKERRQ(ierr); 33 ierr = PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);CHKERRQ(ierr); 34 35 /* Determine file from which we read the matrix A */ 36 ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&data_provided);CHKERRQ(ierr); 37 if (!data_provided) { /* get matrices from PETSc distribution */ 38 ierr = PetscStrncpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/",sizeof(file));CHKERRQ(ierr); 39 if (symm) { 40 #if defined (PETSC_USE_COMPLEX) 41 ierr = PetscStrlcat(file,"hpd-complex-",sizeof(file));CHKERRQ(ierr); 42 #else 43 ierr = PetscStrlcat(file,"spd-real-",sizeof(file));CHKERRQ(ierr); 44 #endif 45 } else { 46 #if defined (PETSC_USE_COMPLEX) 47 ierr = PetscStrlcat(file,"nh-complex-",sizeof(file));CHKERRQ(ierr); 48 #else 49 ierr = PetscStrlcat(file,"ns-real-",sizeof(file));CHKERRQ(ierr); 50 #endif 51 } 52 #if defined(PETSC_USE_64BIT_INDICES) 53 ierr = PetscStrlcat(file,"int64-",sizeof(file));CHKERRQ(ierr); 54 #else 55 ierr = PetscStrlcat(file,"int32-",sizeof(file));CHKERRQ(ierr); 56 #endif 57 #if defined (PETSC_USE_REAL_SINGLE) 58 ierr = PetscStrlcat(file,"float32",sizeof(file));CHKERRQ(ierr); 59 #else 60 ierr = PetscStrlcat(file,"float64",sizeof(file));CHKERRQ(ierr); 61 #endif 62 } 63 /* Load matrix A */ 64 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr); 65 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 66 ierr = MatLoad(A,fd);CHKERRQ(ierr); 67 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 68 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 69 PetscCheckFalse(m != n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 70 71 /* Create dense matrix C and X; C holds true solution with identical columns */ 72 nrhs = 2; 73 ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr); 74 ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr); 75 ierr = MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);CHKERRQ(ierr); 76 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr); 77 ierr = MatSetFromOptions(C);CHKERRQ(ierr); 78 ierr = MatSetUp(C);CHKERRQ(ierr); 79 80 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 81 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 82 ierr = MatSetRandom(C,rand);CHKERRQ(ierr); 83 ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);CHKERRQ(ierr); 84 85 /* Create vectors */ 86 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 87 ierr = VecSetSizes(x,n,PETSC_DECIDE);CHKERRQ(ierr); 88 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 89 ierr = VecDuplicate(x,&b);CHKERRQ(ierr); 90 ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* save the true solution */ 91 92 ierr = PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);CHKERRQ(ierr); 93 switch (isolver) { 94 #if defined(PETSC_HAVE_MUMPS) 95 case 0: 96 ierr = PetscStrcpy(solver,MATSOLVERMUMPS);CHKERRQ(ierr); 97 break; 98 #endif 99 #if defined(PETSC_HAVE_MKL_PARDISO) 100 case 1: 101 ierr = PetscStrcpy(solver,MATSOLVERMKL_PARDISO);CHKERRQ(ierr); 102 break; 103 #endif 104 default: 105 ierr = PetscStrcpy(solver,MATSOLVERPETSC);CHKERRQ(ierr); 106 break; 107 } 108 109 #if defined (PETSC_USE_COMPLEX) 110 if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */ 111 PetscScalar im = PetscSqrtScalar((PetscScalar)-1.); 112 PetscScalar val = -1.0; 113 val = val + im; 114 ierr = MatSetValue(A,1,0,val,INSERT_VALUES);CHKERRQ(ierr); 115 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 116 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 117 } 118 #endif 119 120 ierr = PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);CHKERRQ(ierr); 121 PetscCheckFalse(sratio < 0. || sratio > 1.,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %g", (double)sratio); 122 size_schur = (PetscInt)(sratio*m); 123 124 ierr = PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %" PetscInt_FMT ", sym %d, herm %d, size schur %" PetscInt_FMT ", size mat %" PetscInt_FMT "\n",solver,nrhs,symm,herm,size_schur,m);CHKERRQ(ierr); 125 126 /* Test LU/Cholesky Factorization */ 127 use_lu = PETSC_FALSE; 128 if (!symm) use_lu = PETSC_TRUE; 129 #if defined (PETSC_USE_COMPLEX) 130 if (isolver == 1) use_lu = PETSC_TRUE; 131 #endif 132 if (cuda && symm && !herm) use_lu = PETSC_TRUE; 133 134 if (herm && !use_lu) { /* test also conversion routines inside the solver packages */ 135 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 136 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr); 137 } 138 139 if (use_lu) { 140 ierr = MatGetFactor(A,solver,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 141 } else { 142 if (herm) { 143 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 144 ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); 145 } else { 146 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 147 ierr = MatSetOption(A,MAT_SPD,PETSC_FALSE);CHKERRQ(ierr); 148 } 149 ierr = MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 150 } 151 ierr = ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);CHKERRQ(ierr); 152 ierr = MatFactorSetSchurIS(F,is_schur);CHKERRQ(ierr); 153 154 ierr = ISDestroy(&is_schur);CHKERRQ(ierr); 155 if (use_lu) { 156 ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr); 157 } else { 158 ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr); 159 } 160 161 for (nfact = 0; nfact < 3; nfact++) { 162 Mat AD; 163 164 if (!nfact) { 165 ierr = VecSetRandom(x,rand);CHKERRQ(ierr); 166 if (symm && herm) { 167 ierr = VecAbs(x);CHKERRQ(ierr); 168 } 169 ierr = MatDiagonalSet(A,x,ADD_VALUES);CHKERRQ(ierr); 170 } 171 if (use_lu) { 172 ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr); 173 } else { 174 ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr); 175 } 176 if (cuda) { 177 ierr = MatFactorGetSchurComplement(F,&S,NULL);CHKERRQ(ierr); 178 ierr = MatSetType(S,MATSEQDENSECUDA);CHKERRQ(ierr); 179 ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr); 180 ierr = MatFactorRestoreSchurComplement(F,&S,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr); 181 } 182 ierr = MatFactorCreateSchurComplement(F,&S,NULL);CHKERRQ(ierr); 183 if (!cuda) { 184 ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr); 185 } 186 ierr = VecDuplicate(xschur,&uschur);CHKERRQ(ierr); 187 if (nfact == 1 && (!cuda || (herm && symm))) { 188 ierr = MatFactorInvertSchurComplement(F);CHKERRQ(ierr); 189 } 190 for (nsolve = 0; nsolve < 2; nsolve++) { 191 ierr = VecSetRandom(x,rand);CHKERRQ(ierr); 192 ierr = VecCopy(x,u);CHKERRQ(ierr); 193 194 if (nsolve) { 195 ierr = MatMult(A,x,b);CHKERRQ(ierr); 196 ierr = MatSolve(F,b,x);CHKERRQ(ierr); 197 } else { 198 ierr = MatMultTranspose(A,x,b);CHKERRQ(ierr); 199 ierr = MatSolveTranspose(F,b,x);CHKERRQ(ierr); 200 } 201 /* Check the error */ 202 ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr); /* u <- (-1.0)x + u */ 203 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 204 if (norm > tol) { 205 PetscReal resi; 206 if (nsolve) { 207 ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */ 208 } else { 209 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); /* u = A*x */ 210 } 211 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); /* u <- (-1.0)b + u */ 212 ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr); 213 if (nsolve) { 214 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %" PetscInt_FMT ", s %" PetscInt_FMT ") MatSolve error: Norm of error %g, residual %g\n",nfact,nsolve,(double)norm,(double)resi);CHKERRQ(ierr); 215 } else { 216 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %" PetscInt_FMT ", s %" PetscInt_FMT ") MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,(double)norm,(double)resi);CHKERRQ(ierr); 217 } 218 } 219 ierr = VecSetRandom(xschur,rand);CHKERRQ(ierr); 220 ierr = VecCopy(xschur,uschur);CHKERRQ(ierr); 221 if (nsolve) { 222 ierr = MatMult(S,xschur,bschur);CHKERRQ(ierr); 223 ierr = MatFactorSolveSchurComplement(F,bschur,xschur);CHKERRQ(ierr); 224 } else { 225 ierr = MatMultTranspose(S,xschur,bschur);CHKERRQ(ierr); 226 ierr = MatFactorSolveSchurComplementTranspose(F,bschur,xschur);CHKERRQ(ierr); 227 } 228 /* Check the error */ 229 ierr = VecAXPY(uschur,-1.0,xschur);CHKERRQ(ierr); /* u <- (-1.0)x + u */ 230 ierr = VecNorm(uschur,NORM_2,&norm);CHKERRQ(ierr); 231 if (norm > tol) { 232 PetscReal resi; 233 if (nsolve) { 234 ierr = MatMult(S,xschur,uschur);CHKERRQ(ierr); /* u = A*x */ 235 } else { 236 ierr = MatMultTranspose(S,xschur,uschur);CHKERRQ(ierr); /* u = A*x */ 237 } 238 ierr = VecAXPY(uschur,-1.0,bschur);CHKERRQ(ierr); /* u <- (-1.0)b + u */ 239 ierr = VecNorm(uschur,NORM_2,&resi);CHKERRQ(ierr); 240 if (nsolve) { 241 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %" PetscInt_FMT ", s %" PetscInt_FMT ") MatFactorSolveSchurComplement error: Norm of error %g, residual %g\n",nfact,nsolve,(double)norm,(double)resi);CHKERRQ(ierr); 242 } else { 243 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %" PetscInt_FMT ", s %" PetscInt_FMT ") MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,(double)norm,(double)resi);CHKERRQ(ierr); 244 } 245 } 246 } 247 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);CHKERRQ(ierr); 248 if (!nfact) { 249 ierr = MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);CHKERRQ(ierr); 250 } else { 251 ierr = MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr); 252 } 253 ierr = MatDestroy(&AD);CHKERRQ(ierr); 254 for (nsolve = 0; nsolve < 2; nsolve++) { 255 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr); 256 257 /* Check the error */ 258 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 259 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 260 if (norm > tol) { 261 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %" PetscInt_FMT ", s %" PetscInt_FMT ") MatMatSolve: Norm of error %g\n",nfact,nsolve,(double)norm);CHKERRQ(ierr); 262 } 263 } 264 if (isolver == 0) { 265 Mat spRHS,spRHST,RHST; 266 267 ierr = MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);CHKERRQ(ierr); 268 ierr = MatConvert(RHST,MATSEQAIJ,MAT_INITIAL_MATRIX,&spRHST);CHKERRQ(ierr); 269 ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr); 270 for (nsolve = 0; nsolve < 2; nsolve++) { 271 ierr = MatMatSolve(F,spRHS,X);CHKERRQ(ierr); 272 273 /* Check the error */ 274 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 275 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 276 if (norm > tol) { 277 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %" PetscInt_FMT ", s %" PetscInt_FMT ") sparse MatMatSolve: Norm of error %g\n",nfact,nsolve,(double)norm);CHKERRQ(ierr); 278 } 279 } 280 ierr = MatDestroy(&spRHST);CHKERRQ(ierr); 281 ierr = MatDestroy(&spRHS);CHKERRQ(ierr); 282 ierr = MatDestroy(&RHST);CHKERRQ(ierr); 283 } 284 ierr = MatDestroy(&S);CHKERRQ(ierr); 285 ierr = VecDestroy(&xschur);CHKERRQ(ierr); 286 ierr = VecDestroy(&bschur);CHKERRQ(ierr); 287 ierr = VecDestroy(&uschur);CHKERRQ(ierr); 288 } 289 /* Free data structures */ 290 ierr = MatDestroy(&A);CHKERRQ(ierr); 291 ierr = MatDestroy(&C);CHKERRQ(ierr); 292 ierr = MatDestroy(&F);CHKERRQ(ierr); 293 ierr = MatDestroy(&X);CHKERRQ(ierr); 294 ierr = MatDestroy(&RHS);CHKERRQ(ierr); 295 ierr = PetscRandomDestroy(&rand);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 testset: 306 requires: mkl_pardiso double !complex 307 args: -solver 1 308 309 test: 310 suffix: mkl_pardiso 311 test: 312 requires: cuda 313 suffix: mkl_pardiso_cuda 314 args: -cuda_solve 315 output_file: output/ex192_mkl_pardiso.out 316 test: 317 suffix: mkl_pardiso_1 318 args: -symmetric_solve 319 output_file: output/ex192_mkl_pardiso_1.out 320 test: 321 requires: cuda 322 suffix: mkl_pardiso_cuda_1 323 args: -symmetric_solve -cuda_solve 324 output_file: output/ex192_mkl_pardiso_1.out 325 test: 326 suffix: mkl_pardiso_3 327 args: -symmetric_solve -hermitian_solve 328 output_file: output/ex192_mkl_pardiso_3.out 329 test: 330 requires: cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI) 331 suffix: mkl_pardiso_cuda_3 332 args: -symmetric_solve -hermitian_solve -cuda_solve 333 output_file: output/ex192_mkl_pardiso_3.out 334 335 testset: 336 requires: mumps double !complex 337 args: -solver 0 338 339 test: 340 suffix: mumps 341 test: 342 requires: cuda 343 suffix: mumps_cuda 344 args: -cuda_solve 345 output_file: output/ex192_mumps.out 346 test: 347 suffix: mumps_2 348 args: -symmetric_solve 349 output_file: output/ex192_mumps_2.out 350 test: 351 requires: cuda 352 suffix: mumps_cuda_2 353 args: -symmetric_solve -cuda_solve 354 output_file: output/ex192_mumps_2.out 355 test: 356 suffix: mumps_3 357 args: -symmetric_solve -hermitian_solve 358 output_file: output/ex192_mumps_3.out 359 test: 360 requires: cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI) 361 suffix: mumps_cuda_3 362 args: -symmetric_solve -hermitian_solve -cuda_solve 363 output_file: output/ex192_mumps_3.out 364 365 TEST*/ 366