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: }