Actual source code: ex125.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist).\n\
      3: Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n";

      5: #include <petscmat.h>

      9: int main(int argc,char **args)
     10: {
     11:   Mat            A,RHS,C,F,X;
     12:   Vec            u,x,b;
     14:   PetscMPIInt    rank,nproc;
     15:   PetscInt       i,m,n,nfact,nsolve,nrhs,ipack=0;
     16:   PetscScalar    *array,rval;
     17:   PetscReal      norm,tol=1.e-12;
     18:   IS             perm,iperm;
     19:   MatFactorInfo  info;
     20:   PetscRandom    rand;
     21:   PetscBool      flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE;
     22:   PetscViewer    fd;              /* viewer */
     23:   char           file[PETSC_MAX_PATH_LEN]; /* input file name */

     25:   PetscInitialize(&argc,&args,(char*)0,help);
     26:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
     27:   MPI_Comm_size(PETSC_COMM_WORLD, &nproc);

     29:   /* Determine file from which we read the matrix A */
     30:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
     31:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

     33:   /* Load matrix A */
     34:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
     35:   MatCreate(PETSC_COMM_WORLD,&A);
     36:   MatLoad(A,fd);
     37:   PetscViewerDestroy(&fd);
     38:   MatGetLocalSize(A,&m,&n);
     39:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);

     41:   /* Create dense matrix C and X; C holds true solution with identical colums */
     42:   nrhs = 2;
     43:   PetscOptionsGetInt(NULL,"-nrhs",&nrhs,NULL);
     44:   if (!rank) printf("ex125: nrhs %d\n",nrhs);
     45:   MatCreate(PETSC_COMM_WORLD,&C);
     46:   MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
     47:   MatSetType(C,MATDENSE);
     48:   MatSetFromOptions(C);
     49:   MatSetUp(C);

     51:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
     52:   PetscRandomSetFromOptions(rand);
     53:   MatSetRandom(C,rand);
     54:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);

     56:   /* Create vectors */
     57:   VecCreate(PETSC_COMM_WORLD,&x);
     58:   VecSetSizes(x,n,PETSC_DECIDE);
     59:   VecSetFromOptions(x);
     60:   VecDuplicate(x,&b);
     61:   VecDuplicate(x,&u); /* save the true solution */

     63:   /* Test LU Factorization */
     64:   MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
     65:   /*ISView(perm,PETSC_VIEWER_STDOUT_WORLD);*/
     66:   /*ISView(perm,PETSC_VIEWER_STDOUT_SELF);*/

     68:   PetscOptionsGetInt(NULL,"-mat_solver_package",&ipack,NULL);
     69:   switch (ipack) {
     70: #if defined(PETSC_HAVE_SUPERLU)
     71:   case 0:
     72:     if (!rank) printf(" SUPERLU LU:\n");
     73:     MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
     74:     break;
     75: #endif
     76: #if defined(PETSC_HAVE_SUPERLU_DIST)
     77:   case 1:
     78:     if (!rank) printf(" SUPERLU_DIST LU:\n");
     79:     MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
     80:     break;
     81: #endif
     82: #if defined(PETSC_HAVE_MUMPS)
     83:   case 2:
     84:     if (!rank) printf(" MUMPS LU:\n");
     85:     MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
     86:     {
     87:       /* test mumps options */
     88:       PetscInt icntl_7 = 5;
     89:       MatMumpsSetIcntl(F,7,icntl_7);
     90:     }
     91:     break;
     92: #endif
     93:   default:
     94:     if (!rank) printf(" PETSC LU:\n");
     95:     MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
     96:   }

     98:   info.fill = 5.0;
     99:   info.shifttype = (PetscReal) MAT_SHIFT_NONE;
    100:   MatLUFactorSymbolic(F,A,perm,iperm,&info);

    102:   for (nfact = 0; nfact < 2; nfact++) {
    103:     if (!rank) printf(" %d-the LU numfactorization \n",nfact);
    104:     MatLUFactorNumeric(F,A,&info);

    106:     /* Test MatMatSolve() */
    107:     /*
    108:     if ((ipack == 0 || ipack == 2) && testMatMatSolve) {
    109:       printf("   MatMatSolve() is not implemented for this package. Skip the testing.\n");
    110:       testMatMatSolve = PETSC_FALSE;
    111:     }
    112:      */
    113:     if (testMatMatSolve) {
    114:       if (!nfact) {
    115:         MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
    116:       } else {
    117:         MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
    118:       }
    119:       for (nsolve = 0; nsolve < 2; nsolve++) {
    120:         if (!rank) printf("   %d-the MatMatSolve \n",nsolve);
    121:         MatMatSolve(F,RHS,X);

    123:         /* Check the error */
    124:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
    125:         MatNorm(X,NORM_FROBENIUS,&norm);
    126:         if (norm > tol) {
    127:           if (!rank) {
    128:             PetscPrintf(PETSC_COMM_SELF,"1st MatMatSolve: Norm of error %g, nsolve %d\n",norm,nsolve);
    129:           }
    130:         }
    131:       }
    132:     }

    134:     /* Test MatSolve() */
    135:     if (testMatSolve) {
    136:       for (nsolve = 0; nsolve < 2; nsolve++) {
    137:         VecGetArray(x,&array);
    138:         for (i=0; i<m; i++) {
    139:           PetscRandomGetValue(rand,&rval);
    140:           array[i] = rval;
    141:         }
    142:         VecRestoreArray(x,&array);
    143:         VecCopy(x,u);
    144:         MatMult(A,x,b);

    146:         if (!rank) printf("   %d-the MatSolve \n",nsolve);
    147:         MatSolve(F,b,x);

    149:         /* Check the error */
    150:         VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
    151:         VecNorm(u,NORM_2,&norm);
    152:         if (norm > tol) {
    153:           MatMult(A,x,u); /* u = A*x */
    154:           PetscReal resi;
    155:           VecAXPY(u,-1.0,b);  /* u <- (-1.0)b + u */
    156:           VecNorm(u,NORM_2,&resi);
    157:           if (!rank) {
    158:             PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g, resi %g, LU numfact %d\n",norm,resi,nfact);
    159:           }
    160:         }
    161:       }
    162:     }
    163:   }

    165:   /* Free data structures */
    166:   MatDestroy(&A);
    167:   MatDestroy(&C);
    168:   MatDestroy(&F);
    169:   MatDestroy(&X);
    170:   if (testMatMatSolve) {
    171:     MatDestroy(&RHS);
    172:   }

    174:   PetscRandomDestroy(&rand);
    175:   ISDestroy(&perm);
    176:   ISDestroy(&iperm);
    177:   VecDestroy(&x);
    178:   VecDestroy(&b);
    179:   VecDestroy(&u);
    180:   PetscFinalize();
    181:   return 0;
    182: }