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 if (size > 1) SETERRQ(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 if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", 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 if (sratio < 0. || sratio > 1.) { 122 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio); 123 } 124 size_schur = (PetscInt)(sratio*m); 125 126 ierr = PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %D, sym %d, herm %d, size schur %D, size mat %D\n",solver,nrhs,symm,herm,size_schur,m);CHKERRQ(ierr); 127 128 /* Test LU/Cholesky Factorization */ 129 use_lu = PETSC_FALSE; 130 if (!symm) use_lu = PETSC_TRUE; 131 #if defined (PETSC_USE_COMPLEX) 132 if (isolver == 1) use_lu = PETSC_TRUE; 133 #endif 134 if (cuda && symm && !herm) use_lu = PETSC_TRUE; 135 136 if (herm && !use_lu) { /* test also conversion routines inside the solver packages */ 137 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 138 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr); 139 } 140 141 if (use_lu) { 142 ierr = MatGetFactor(A,solver,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 143 } else { 144 if (herm) { 145 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 146 ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); 147 } else { 148 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 149 ierr = MatSetOption(A,MAT_SPD,PETSC_FALSE);CHKERRQ(ierr); 150 } 151 ierr = MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 152 } 153 ierr = ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);CHKERRQ(ierr); 154 ierr = MatFactorSetSchurIS(F,is_schur);CHKERRQ(ierr); 155 156 ierr = ISDestroy(&is_schur);CHKERRQ(ierr); 157 if (use_lu) { 158 ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr); 159 } else { 160 ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr); 161 } 162 163 for (nfact = 0; nfact < 3; nfact++) { 164 Mat AD; 165 166 if (!nfact) { 167 ierr = VecSetRandom(x,rand);CHKERRQ(ierr); 168 if (symm && herm) { 169 ierr = VecAbs(x);CHKERRQ(ierr); 170 } 171 ierr = MatDiagonalSet(A,x,ADD_VALUES);CHKERRQ(ierr); 172 } 173 if (use_lu) { 174 ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr); 175 } else { 176 ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr); 177 } 178 if (cuda) { 179 ierr = MatFactorGetSchurComplement(F,&S,NULL);CHKERRQ(ierr); 180 ierr = MatSetType(S,MATSEQDENSECUDA);CHKERRQ(ierr); 181 ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr); 182 ierr = MatFactorRestoreSchurComplement(F,&S,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr); 183 } 184 ierr = MatFactorCreateSchurComplement(F,&S,NULL);CHKERRQ(ierr); 185 if (!cuda) { 186 ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr); 187 } 188 ierr = VecDuplicate(xschur,&uschur);CHKERRQ(ierr); 189 if (nfact == 1 && (!cuda || (herm && symm))) { 190 ierr = MatFactorInvertSchurComplement(F);CHKERRQ(ierr); 191 } 192 for (nsolve = 0; nsolve < 2; nsolve++) { 193 ierr = VecSetRandom(x,rand);CHKERRQ(ierr); 194 ierr = VecCopy(x,u);CHKERRQ(ierr); 195 196 if (nsolve) { 197 ierr = MatMult(A,x,b);CHKERRQ(ierr); 198 ierr = MatSolve(F,b,x);CHKERRQ(ierr); 199 } else { 200 ierr = MatMultTranspose(A,x,b);CHKERRQ(ierr); 201 ierr = MatSolveTranspose(F,b,x);CHKERRQ(ierr); 202 } 203 /* Check the error */ 204 ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr); /* u <- (-1.0)x + u */ 205 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 206 if (norm > tol) { 207 PetscReal resi; 208 if (nsolve) { 209 ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */ 210 } else { 211 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); /* u = A*x */ 212 } 213 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); /* u <- (-1.0)b + u */ 214 ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr); 215 if (nsolve) { 216 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr); 217 } else { 218 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr); 219 } 220 } 221 ierr = VecSetRandom(xschur,rand);CHKERRQ(ierr); 222 ierr = VecCopy(xschur,uschur);CHKERRQ(ierr); 223 if (nsolve) { 224 ierr = MatMult(S,xschur,bschur);CHKERRQ(ierr); 225 ierr = MatFactorSolveSchurComplement(F,bschur,xschur);CHKERRQ(ierr); 226 } else { 227 ierr = MatMultTranspose(S,xschur,bschur);CHKERRQ(ierr); 228 ierr = MatFactorSolveSchurComplementTranspose(F,bschur,xschur);CHKERRQ(ierr); 229 } 230 /* Check the error */ 231 ierr = VecAXPY(uschur,-1.0,xschur);CHKERRQ(ierr); /* u <- (-1.0)x + u */ 232 ierr = VecNorm(uschur,NORM_2,&norm);CHKERRQ(ierr); 233 if (norm > tol) { 234 PetscReal resi; 235 if (nsolve) { 236 ierr = MatMult(S,xschur,uschur);CHKERRQ(ierr); /* u = A*x */ 237 } else { 238 ierr = MatMultTranspose(S,xschur,uschur);CHKERRQ(ierr); /* u = A*x */ 239 } 240 ierr = VecAXPY(uschur,-1.0,bschur);CHKERRQ(ierr); /* u <- (-1.0)b + u */ 241 ierr = VecNorm(uschur,NORM_2,&resi);CHKERRQ(ierr); 242 if (nsolve) { 243 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr); 244 } else { 245 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr); 246 } 247 } 248 } 249 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);CHKERRQ(ierr); 250 if (!nfact) { 251 ierr = MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);CHKERRQ(ierr); 252 } else { 253 ierr = MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr); 254 } 255 ierr = MatDestroy(&AD);CHKERRQ(ierr); 256 for (nsolve = 0; nsolve < 2; nsolve++) { 257 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr); 258 259 /* Check the error */ 260 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 261 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 262 if (norm > tol) { 263 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);CHKERRQ(ierr); 264 } 265 } 266 if (isolver == 0) { 267 Mat spRHS,spRHST,RHST; 268 269 ierr = MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);CHKERRQ(ierr); 270 ierr = MatConvert(RHST,MATSEQAIJ,MAT_INITIAL_MATRIX,&spRHST);CHKERRQ(ierr); 271 ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr); 272 for (nsolve = 0; nsolve < 2; nsolve++) { 273 ierr = MatMatSolve(F,spRHS,X);CHKERRQ(ierr); 274 275 /* Check the error */ 276 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 277 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 278 if (norm > tol) { 279 ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) sparse MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);CHKERRQ(ierr); 280 } 281 } 282 ierr = MatDestroy(&spRHST);CHKERRQ(ierr); 283 ierr = MatDestroy(&spRHS);CHKERRQ(ierr); 284 ierr = MatDestroy(&RHST);CHKERRQ(ierr); 285 } 286 ierr = MatDestroy(&S);CHKERRQ(ierr); 287 ierr = VecDestroy(&xschur);CHKERRQ(ierr); 288 ierr = VecDestroy(&bschur);CHKERRQ(ierr); 289 ierr = VecDestroy(&uschur);CHKERRQ(ierr); 290 } 291 /* Free data structures */ 292 ierr = MatDestroy(&A);CHKERRQ(ierr); 293 ierr = MatDestroy(&C);CHKERRQ(ierr); 294 ierr = MatDestroy(&F);CHKERRQ(ierr); 295 ierr = MatDestroy(&X);CHKERRQ(ierr); 296 ierr = MatDestroy(&RHS);CHKERRQ(ierr); 297 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 298 ierr = VecDestroy(&x);CHKERRQ(ierr); 299 ierr = VecDestroy(&b);CHKERRQ(ierr); 300 ierr = VecDestroy(&u);CHKERRQ(ierr); 301 ierr = PetscFinalize(); 302 return ierr; 303 } 304 305 /*TEST 306 307 testset: 308 requires: mkl_pardiso double !complex 309 args: -solver 1 310 311 test: 312 suffix: mkl_pardiso 313 test: 314 requires: cuda 315 suffix: mkl_pardiso_cuda 316 args: -cuda_solve 317 output_file: output/ex192_mkl_pardiso.out 318 test: 319 suffix: mkl_pardiso_1 320 args: -symmetric_solve 321 output_file: output/ex192_mkl_pardiso_1.out 322 test: 323 requires: cuda 324 suffix: mkl_pardiso_cuda_1 325 args: -symmetric_solve -cuda_solve 326 output_file: output/ex192_mkl_pardiso_1.out 327 test: 328 suffix: mkl_pardiso_3 329 args: -symmetric_solve -hermitian_solve 330 output_file: output/ex192_mkl_pardiso_3.out 331 test: 332 requires: cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI) 333 suffix: mkl_pardiso_cuda_3 334 args: -symmetric_solve -hermitian_solve -cuda_solve 335 output_file: output/ex192_mkl_pardiso_3.out 336 337 testset: 338 requires: mumps double !complex 339 args: -solver 0 340 341 test: 342 suffix: mumps 343 test: 344 requires: cuda 345 suffix: mumps_cuda 346 args: -cuda_solve 347 output_file: output/ex192_mumps.out 348 test: 349 suffix: mumps_2 350 args: -symmetric_solve 351 output_file: output/ex192_mumps_2.out 352 test: 353 requires: cuda 354 suffix: mumps_cuda_2 355 args: -symmetric_solve -cuda_solve 356 output_file: output/ex192_mumps_2.out 357 test: 358 suffix: mumps_3 359 args: -symmetric_solve -hermitian_solve 360 output_file: output/ex192_mumps_3.out 361 test: 362 requires: cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI) 363 suffix: mumps_cuda_3 364 args: -symmetric_solve -hermitian_solve -cuda_solve 365 output_file: output/ex192_mumps_3.out 366 367 TEST*/ 368