Actual source code: ex125.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  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,size;
 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, &size);

 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;
 89:       PetscReal   cntl;
 90: 
 91:       icntl = 2;        /* sequential matrix ordering */
 92:       MatMumpsSetIcntl(F,7,icntl);

 94:       cntl = 1.e-6; /* threshhold for row pivot detection */
 95:       MatMumpsSetIcntl(F,24,1);
 96:       MatMumpsSetCntl(F,3,cntl);
 97:     }
 98:     break;
 99: #endif
100:   default:
101:     if (!rank) printf(" PETSC LU:\n");
102:     MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
103:   }

105:   MatFactorInfoInitialize(&info);
106:   info.fill      = 5.0;
107:   info.shifttype = (PetscReal) MAT_SHIFT_NONE;
108:   MatLUFactorSymbolic(F,A,perm,iperm,&info);

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

114:     /* Test MatMatSolve() */
115:     /*
116:     if ((ipack == 0 || ipack == 2) && testMatMatSolve) {
117:       printf("   MatMatSolve() is not implemented for this package. Skip the testing.\n");
118:       testMatMatSolve = PETSC_FALSE;
119:     }
120:      */
121:     if (testMatMatSolve) {
122:       if (!nfact) {
123:         MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
124:       } else {
125:         MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
126:       }
127:       for (nsolve = 0; nsolve < 2; nsolve++) {
128:         if (!rank) printf("   %d-the MatMatSolve \n",nsolve);
129:         MatMatSolve(F,RHS,X);

131:         /* Check the error */
132:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
133:         MatNorm(X,NORM_FROBENIUS,&norm);
134:         if (norm > tol) {
135:           if (!rank) {
136:             PetscPrintf(PETSC_COMM_SELF,"1st MatMatSolve: Norm of error %g, nsolve %d\n",norm,nsolve);
137:           }
138:         }
139:       }
140:     }

142:     /* Test MatSolve() */
143:     if (testMatSolve) {
144:       for (nsolve = 0; nsolve < 2; nsolve++) {
145:         VecGetArray(x,&array);
146:         for (i=0; i<m; i++) {
147:           PetscRandomGetValue(rand,&rval);
148:           array[i] = rval;
149:         }
150:         VecRestoreArray(x,&array);
151:         VecCopy(x,u);
152:         MatMult(A,x,b);

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

157:         /* Check the error */
158:         VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
159:         VecNorm(u,NORM_2,&norm);
160:         if (norm > tol) {
161:           PetscReal resi;
162:           MatMult(A,x,u); /* u = A*x */
163:           VecAXPY(u,-1.0,b);  /* u <- (-1.0)b + u */
164:           VecNorm(u,NORM_2,&resi);
165:           if (!rank) {
166:             PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g, resi %g, LU numfact %d\n",norm,resi,nfact);
167:           }
168:         }
169:       }
170:     }
171:   }

173:   /* Free data structures */
174:   MatDestroy(&A);
175:   MatDestroy(&C);
176:   MatDestroy(&F);
177:   MatDestroy(&X);
178:   if (testMatMatSolve) {
179:     MatDestroy(&RHS);
180:   }

182:   PetscRandomDestroy(&rand);
183:   ISDestroy(&perm);
184:   ISDestroy(&iperm);
185:   VecDestroy(&x);
186:   VecDestroy(&b);
187:   VecDestroy(&u);
188:   PetscFinalize();
189:   return 0;
190: }