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