Actual source code: ex192.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      2: static char help[] = "Tests MatSolve() and MatMatSolve() with MUMPS or MKL_PARDISO sequential solvers in Schur complement mode.\n\
      3: Example: mpiexec -n 1 ./ex192 -f <matrix binary file> -nrhs 4 -symmetric_solve -hermitian_solve -schur_ratio 0.3\n\n";

      5:  #include <petscmat.h>

      7: int main(int argc,char **args)
      8: {
      9:   Mat            A,RHS,C,F,X,S;
     10:   Vec            u,x,b;
     11:   Vec            xschur,bschur,uschur;
     12:   IS             is_schur;
     14:   PetscMPIInt    size;
     15:   PetscInt       isolver=0,size_schur,m,n,nfact,nsolve,nrhs;
     16:   PetscReal      norm,tol=PETSC_SQRT_MACHINE_EPSILON;
     17:   PetscRandom    rand;
     18:   PetscBool      data_provided,herm,symm,use_lu,cuda = PETSC_FALSE;
     19:   PetscReal      sratio = 5.1/12.;
     20:   PetscViewer    fd;              /* viewer */
     21:   char           solver[256];
     22:   char           file[PETSC_MAX_PATH_LEN]; /* input file name */

     24:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     25:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
     26:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor test");
     27:   /* Determine which type of solver we want to test for */
     28:   herm = PETSC_FALSE;
     29:   symm = PETSC_FALSE;
     30:   PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);
     31:   PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);
     32:   if (herm) symm = PETSC_TRUE;
     33:   PetscOptionsGetBool(NULL,NULL,"-cuda_solve",&cuda,NULL);
     34:   PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);

     36:   /* Determine file from which we read the matrix A */
     37:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&data_provided);
     38:   if (!data_provided) { /* get matrices from PETSc distribution */
     39:     PetscStrncpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/",sizeof(file));
     40:     if (symm) {
     41: #if defined (PETSC_USE_COMPLEX)
     42:       PetscStrlcat(file,"hpd-complex-",sizeof(file));
     43: #else
     44:       PetscStrlcat(file,"spd-real-",sizeof(file));
     45: #endif
     46:     } else {
     47: #if defined (PETSC_USE_COMPLEX)
     48:       PetscStrlcat(file,"nh-complex-",sizeof(file));
     49: #else
     50:       PetscStrlcat(file,"ns-real-",sizeof(file));
     51: #endif
     52:     }
     53: #if defined(PETSC_USE_64BIT_INDICES)
     54:     PetscStrlcat(file,"int64-",sizeof(file));
     55: #else
     56:     PetscStrlcat(file,"int32-",sizeof(file));
     57: #endif
     58: #if defined (PETSC_USE_REAL_SINGLE)
     59:     PetscStrlcat(file,"float32",sizeof(file));
     60: #else
     61:     PetscStrlcat(file,"float64",sizeof(file));
     62: #endif
     63:   }
     64:   /* Load matrix A */
     65:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
     66:   MatCreate(PETSC_COMM_WORLD,&A);
     67:   MatLoad(A,fd);
     68:   PetscViewerDestroy(&fd);
     69:   MatGetSize(A,&m,&n);
     70:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);

     72:   /* Create dense matrix C and X; C holds true solution with identical colums */
     73:   nrhs = 2;
     74:   PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
     75:   MatCreate(PETSC_COMM_WORLD,&C);
     76:   MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
     77:   MatSetType(C,MATDENSE);
     78:   MatSetFromOptions(C);
     79:   MatSetUp(C);

     81:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
     82:   PetscRandomSetFromOptions(rand);
     83:   MatSetRandom(C,rand);
     84:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);

     86:   /* Create vectors */
     87:   VecCreate(PETSC_COMM_WORLD,&x);
     88:   VecSetSizes(x,n,PETSC_DECIDE);
     89:   VecSetFromOptions(x);
     90:   VecDuplicate(x,&b);
     91:   VecDuplicate(x,&u); /* save the true solution */

     93:   PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);
     94:   switch (isolver) {
     95: #if defined(PETSC_HAVE_MUMPS)
     96:     case 0:
     97:       PetscStrcpy(solver,MATSOLVERMUMPS);
     98:       break;
     99: #endif
    100: #if defined(PETSC_HAVE_MKL_PARDISO)
    101:     case 1:
    102:       PetscStrcpy(solver,MATSOLVERMKL_PARDISO);
    103:       break;
    104: #endif
    105:     default:
    106:       PetscStrcpy(solver,MATSOLVERPETSC);
    107:       break;
    108:   }

    110: #if defined (PETSC_USE_COMPLEX)
    111:   if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
    112:     PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
    113:     PetscScalar val = -1.0;
    114:     val = val + im;
    115:     MatSetValue(A,1,0,val,INSERT_VALUES);
    116:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    117:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    118:   }
    119: #endif

    121:   PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);
    122:   if (sratio < 0. || sratio > 1.) {
    123:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
    124:   }
    125:   size_schur = (PetscInt)(sratio*m);

    127:   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);

    129:   /* Test LU/Cholesky Factorization */
    130:   use_lu = PETSC_FALSE;
    131:   if (!symm) use_lu = PETSC_TRUE;
    132: #if defined (PETSC_USE_COMPLEX)
    133:   if (isolver == 1) use_lu = PETSC_TRUE;
    134: #endif
    135:   if (cuda && symm && !herm) use_lu = PETSC_TRUE;

    137:   if (herm && !use_lu) { /* test also conversion routines inside the solver packages */
    138:     MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
    139:     MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);
    140:   }

    142:   if (use_lu) {
    143:     MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
    144:   } else {
    145:     if (herm) {
    146:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
    147:       MatSetOption(A,MAT_SPD,PETSC_TRUE);
    148:     } else {
    149:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
    150:       MatSetOption(A,MAT_SPD,PETSC_FALSE);
    151:     }
    152:     MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
    153:   }
    154:   ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);
    155:   MatFactorSetSchurIS(F,is_schur);

    157:   ISDestroy(&is_schur);
    158:   if (use_lu) {
    159:     MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
    160:   } else {
    161:     MatCholeskyFactorSymbolic(F,A,NULL,NULL);
    162:   }

    164:   for (nfact = 0; nfact < 3; nfact++) {
    165:     Mat AD;

    167:     if (!nfact) {
    168:       VecSetRandom(x,rand);
    169:       if (symm && herm) {
    170:         VecAbs(x);
    171:       }
    172:       MatDiagonalSet(A,x,ADD_VALUES);
    173:     }
    174:     if (use_lu) {
    175:       MatLUFactorNumeric(F,A,NULL);
    176:     } else {
    177:       MatCholeskyFactorNumeric(F,A,NULL);
    178:     }
    179:     if (cuda) {
    180:       MatFactorGetSchurComplement(F,&S,NULL);
    181:       MatSetType(S,MATSEQDENSECUDA);
    182:       MatCreateVecs(S,&xschur,&bschur);
    183:       MatFactorRestoreSchurComplement(F,&S,MAT_FACTOR_SCHUR_UNFACTORED);
    184:     }
    185:     MatFactorCreateSchurComplement(F,&S,NULL);
    186:     if (!cuda) {
    187:       MatCreateVecs(S,&xschur,&bschur);
    188:     }
    189:     VecDuplicate(xschur,&uschur);
    190:     if (nfact == 1 && (!cuda || (herm && symm))) {
    191:       MatFactorInvertSchurComplement(F);
    192:     }
    193:     for (nsolve = 0; nsolve < 2; nsolve++) {
    194:       VecSetRandom(x,rand);
    195:       VecCopy(x,u);

    197:       if (nsolve) {
    198:         MatMult(A,x,b);
    199:         MatSolve(F,b,x);
    200:       } else {
    201:         MatMultTranspose(A,x,b);
    202:         MatSolveTranspose(F,b,x);
    203:       }
    204:       /* Check the error */
    205:       VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
    206:       VecNorm(u,NORM_2,&norm);
    207:       if (norm > tol) {
    208:         PetscReal resi;
    209:         if (nsolve) {
    210:           MatMult(A,x,u); /* u = A*x */
    211:         } else {
    212:           MatMultTranspose(A,x,u); /* u = A*x */
    213:         }
    214:         VecAXPY(u,-1.0,b);  /* u <- (-1.0)b + u */
    215:         VecNorm(u,NORM_2,&resi);
    216:         if (nsolve) {
    217:           PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
    218:         } else {
    219:           PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
    220:         }
    221:       }
    222:       VecSetRandom(xschur,rand);
    223:       VecCopy(xschur,uschur);
    224:       if (nsolve) {
    225:         MatMult(S,xschur,bschur);
    226:         MatFactorSolveSchurComplement(F,bschur,xschur);
    227:       } else {
    228:         MatMultTranspose(S,xschur,bschur);
    229:         MatFactorSolveSchurComplementTranspose(F,bschur,xschur);
    230:       }
    231:       /* Check the error */
    232:       VecAXPY(uschur,-1.0,xschur);  /* u <- (-1.0)x + u */
    233:       VecNorm(uschur,NORM_2,&norm);
    234:       if (norm > tol) {
    235:         PetscReal resi;
    236:         if (nsolve) {
    237:           MatMult(S,xschur,uschur); /* u = A*x */
    238:         } else {
    239:           MatMultTranspose(S,xschur,uschur); /* u = A*x */
    240:         }
    241:         VecAXPY(uschur,-1.0,bschur);  /* u <- (-1.0)b + u */
    242:         VecNorm(uschur,NORM_2,&resi);
    243:         if (nsolve) {
    244:           PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
    245:         } else {
    246:           PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
    247:         }
    248:       }
    249:     }
    250:     MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);
    251:     if (!nfact) {
    252:       MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);
    253:     } else {
    254:       MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);
    255:     }
    256:     MatDestroy(&AD);
    257:     for (nsolve = 0; nsolve < 2; nsolve++) {
    258:       MatMatSolve(F,RHS,X);

    260:       /* Check the error */
    261:       MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
    262:       MatNorm(X,NORM_FROBENIUS,&norm);
    263:       if (norm > tol) {
    264:         PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
    265:       }
    266:     }
    267:     if (isolver == 0) {
    268:       Mat spRHS,spRHST,RHST;

    270:       MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
    271:       MatConvert(RHST,MATSEQAIJ,MAT_INITIAL_MATRIX,&spRHST);
    272:       MatCreateTranspose(spRHST,&spRHS);
    273:       for (nsolve = 0; nsolve < 2; nsolve++) {
    274:         MatMatSolve(F,spRHS,X);

    276:         /* Check the error */
    277:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
    278:         MatNorm(X,NORM_FROBENIUS,&norm);
    279:         if (norm > tol) {
    280:           PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) sparse MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
    281:         }
    282:       }
    283:       MatDestroy(&spRHST);
    284:       MatDestroy(&spRHS);
    285:       MatDestroy(&RHST);
    286:     }
    287:     MatDestroy(&S);
    288:     VecDestroy(&xschur);
    289:     VecDestroy(&bschur);
    290:     VecDestroy(&uschur);
    291:   }
    292:   /* Free data structures */
    293:   MatDestroy(&A);
    294:   MatDestroy(&C);
    295:   MatDestroy(&F);
    296:   MatDestroy(&X);
    297:   MatDestroy(&RHS);
    298:   PetscRandomDestroy(&rand);
    299:   VecDestroy(&x);
    300:   VecDestroy(&b);
    301:   VecDestroy(&u);
    302:   PetscFinalize();
    303:   return ierr;
    304: }

    307: /*TEST

    309:    testset:
    310:      requires: mkl_pardiso double !complex !define(PETSC_USE_64BIT_INDICES)
    311:      args: -solver 1

    313:      test:
    314:        suffix: mkl_pardiso
    315:      test:
    316:        requires: cuda
    317:        suffix: mkl_pardiso_cuda
    318:        args: -cuda_solve
    319:        output_file: output/ex192_mkl_pardiso.out
    320:      test:
    321:        suffix: mkl_pardiso_1
    322:        args: -symmetric_solve
    323:        output_file: output/ex192_mkl_pardiso_1.out
    324:      test:
    325:        requires: cuda
    326:        suffix: mkl_pardiso_cuda_1
    327:        args: -symmetric_solve -cuda_solve
    328:        output_file: output/ex192_mkl_pardiso_1.out
    329:      test:
    330:        suffix: mkl_pardiso_3
    331:        args: -symmetric_solve -hermitian_solve
    332:        output_file: output/ex192_mkl_pardiso_3.out
    333:      test:
    334:        requires: cuda
    335:        suffix: mkl_pardiso_cuda_3
    336:        args: -symmetric_solve -hermitian_solve -cuda_solve
    337:        output_file: output/ex192_mkl_pardiso_3.out

    339:    testset:
    340:      requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
    341:      args: -solver 0

    343:      test:
    344:        suffix: mumps
    345:      test:
    346:        requires: cuda
    347:        suffix: mumps_cuda
    348:        args: -cuda_solve
    349:        output_file: output/ex192_mumps.out
    350:      test:
    351:        suffix: mumps_2
    352:        args: -symmetric_solve
    353:        output_file: output/ex192_mumps_2.out
    354:      test:
    355:        requires: cuda
    356:        suffix: mumps_cuda_2
    357:        args: -symmetric_solve -cuda_solve
    358:        output_file: output/ex192_mumps_2.out
    359:      test:
    360:        suffix: mumps_3
    361:        args: -symmetric_solve -hermitian_solve
    362:        output_file: output/ex192_mumps_3.out
    363:      test:
    364:        requires: cuda
    365:        suffix: mumps_cuda_3
    366:        args: -symmetric_solve -hermitian_solve -cuda_solve
    367:        output_file: output/ex192_mumps_3.out

    369: TEST*/