Actual source code: ex125.c

petsc-3.3-p6 2013-02-11
  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:   MatSetUp(C);
 52: 
 53:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 54:   PetscRandomSetFromOptions(rand);
 55:   MatGetArray(C,&array);
 56:   for (i=0; i<m; i++){
 57:     PetscRandomGetValue(rand,&rval);
 58:     array[i] = rval;
 59:   }
 60:   if (nrhs > 1){
 61:     for (k=1; k<nrhs; k++){
 62:       for (i=0; i<m; i++){
 63:         array[m*k+i] = array[i];
 64:       }
 65:     }
 66:   }
 67:   MatRestoreArray(C,&array);
 68:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 70:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
 71: 
 72:   /* Create vectors */
 73:   VecCreate(PETSC_COMM_WORLD,&x);
 74:   VecSetSizes(x,n,PETSC_DECIDE);
 75:   VecSetFromOptions(x);
 76:   VecDuplicate(x,&b);
 77:   VecDuplicate(x,&u); /* save the true solution */

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

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

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

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

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

161:         if (!rank) printf("   %d-the MatSolve \n",nsolve);
162:         MatSolve(F,b,x);
163: 
164:         /* Check the error */
165:         VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
166:         VecNorm(u,NORM_2,&norm);
167:         if (norm > tol){
168:           MatMult(A,x,u); /* u = A*x */
169:           PetscReal resi;
170:           VecAXPY(u,-1.0,b);  /* u <- (-1.0)b + u */
171:           VecNorm(u,NORM_2,&resi);
172:           if (!rank){
173:             PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g, resi %g, LU numfact %d\n",norm,resi,nfact);
174:           }
175:         }
176:       }
177:     }
178:   }
179: 
180:   /* Free data structures */
181:   MatDestroy(&A);
182:   MatDestroy(&C);
183:   MatDestroy(&F);
184:   MatDestroy(&X);
185:   if (testMatMatSolve){
186:     MatDestroy(&RHS);
187:   }
188: 
189:   PetscRandomDestroy(&rand);
190:   ISDestroy(&perm);
191:   ISDestroy(&iperm);
192:   VecDestroy(&x);
193:   VecDestroy(&b);
194:   VecDestroy(&u);
195:   PetscFinalize();
196:   return 0;
197: }