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