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