Actual source code: ex125.c

  1: 
  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,k,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(PETSC_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) {
 40:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
 41:   }
 42: 
 43:   /* Create dense matrix C and X; C holds true solution with identical colums */
 44:   nrhs = 2;
 45:   PetscOptionsGetInt(PETSC_NULL,"-nrhs",&nrhs,PETSC_NULL);
 46:   if (!rank) printf("ex125: nrhs %d\n",nrhs);
 47:   MatCreate(PETSC_COMM_WORLD,&C);
 48:   MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
 49:   MatSetType(C,MATDENSE);
 50:   MatSetFromOptions(C);
 51: 
 52:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 53:   PetscRandomSetFromOptions(rand);
 54:   MatGetArray(C,&array);
 55:   for (i=0; i<m; i++){
 56:     PetscRandomGetValue(rand,&rval);
 57:     array[i] = rval;
 58:   }
 59:   if (nrhs > 1){
 60:     for (k=1; k<nrhs; k++){
 61:       for (i=0; i<m; i++){
 62:         array[m*k+i] = array[i];
 63:       }
 64:     }
 65:   }
 66:   MatRestoreArray(C,&array);
 67:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 69:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
 70: 
 71:   /* Create vectors */
 72:   VecCreate(PETSC_COMM_WORLD,&x);
 73:   VecSetSizes(x,n,PETSC_DECIDE);
 74:   VecSetFromOptions(x);
 75:   VecDuplicate(x,&b);
 76:   VecDuplicate(x,&u); /* save the true solution */

 78:   /* Test LU Factorization */
 79:   MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
 80:   //ISView(perm,PETSC_VIEWER_STDOUT_WORLD);
 81:   //ISView(perm,PETSC_VIEWER_STDOUT_SELF);
 82: 
 83:   PetscOptionsGetInt(PETSC_NULL,"-mat_solver_package",&ipack,PETSC_NULL);
 84:   switch (ipack){
 85:   case 0:
 86: #ifdef PETSC_HAVE_SUPERLU
 87:     if (!rank) printf(" SUPERLU LU:\n");
 88:     MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
 89:     break;
 90: #endif
 91:   case 1:
 92: #ifdef PETSC_HAVE_SUPERLU_DIST
 93:     if (!rank) printf(" SUPERLU_DIST LU:\n");
 94:     MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
 95:     break;
 96: #endif
 97:   case 2:
 98: #ifdef PETSC_HAVE_MUMPS 
 99:     if (!rank) printf(" MUMPS LU:\n");
100:     MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
101:     {
102:       /* test mumps options */
103:       PetscInt icntl_7 = 5;
104:       MatMumpsSetIcntl(F,7,icntl_7);
105:     }
106:     break;
107: #endif
108:   default:
109:     if (!rank) printf(" PETSC LU:\n");
110:     MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
111:   }

113:   info.fill = 5.0;
114:   MatLUFactorSymbolic(F,A,perm,iperm,&info);

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

120:     /* Test MatMatSolve() */
121:     if ((ipack == 0 || ipack == 2) && testMatMatSolve){
122:       printf("   MaMattSolve() is not implemented for this package. Skip the testing.\n");
123:       testMatMatSolve = PETSC_FALSE;
124:     }
125:     if (testMatMatSolve){
126:       if (!nfact){
127:         MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
128:       } else {
129:         MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
130:       }
131:       for (nsolve = 0; nsolve < 2; nsolve++){
132:         if (!rank) printf("   %d-the MatMatSolve \n",nsolve);
133:         MatMatSolve(F,RHS,X);
134: 
135:         /* Check the error */
136:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
137:         MatNorm(X,NORM_FROBENIUS,&norm);
138:         if (norm > tol){
139:           if (!rank){
140:             PetscPrintf(PETSC_COMM_SELF,"1st MatMatSolve: Norm of error %g, nsolve %d\n",norm,nsolve);
141:           }
142:         }
143:       }
144:     }

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

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

176:     /* Test MatMatSolve() */
177:     if (testMatMatSolve){
178:       for (nsolve = 0; nsolve < 2; nsolve++){
179:         if (!rank) printf("   %d-the MatMatSolve \n",nsolve);
180:         MatMatSolve(F,RHS,X);
181: 
182:         /* Check the error */
183:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
184:         MatNorm(X,NORM_FROBENIUS,&norm);
185:         if (norm > tol){
186:           if (!rank){
187:             PetscPrintf(PETSC_COMM_SELF,"2nd MatMatSolve: Norm of error %g, nsolve %d\n",norm,nsolve);
188:           }
189:         }
190:       }
191:     }
192:   }
193: 
194:   /* Free data structures */
195:   MatDestroy(&A);
196:   MatDestroy(&C);
197:   MatDestroy(&F);
198:   MatDestroy(&X);
199:   if (testMatMatSolve){
200:     MatDestroy(&RHS);
201:   }
202: 
203:   PetscRandomDestroy(&rand);
204:   ISDestroy(&perm);
205:   ISDestroy(&iperm);
206:   VecDestroy(&x);
207:   VecDestroy(&b);
208:   VecDestroy(&u);
209:   PetscFinalize();
210:   return 0;
211: }