Actual source code: ex27.c

petsc-dev 2014-08-28
Report Typos and Errors
  2: static char help[] = "Reads a PETSc matrix and vector from a file and solves the normal equations.\n\n";
  3: /*T
  4:    Concepts: KSP^solving a linear system
  5:    Concepts: Normal equations
  6:    Processors: n
  7: T*/

  9: /*
 10:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 11:   automatically includes:
 12:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 13:      petscmat.h - matrices
 14:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 15:      petscviewer.h - viewers               petscpc.h  - preconditioners
 16: */
 17: #include <petscksp.h>


 22: int main(int argc,char **args)
 23: {
 24:   KSP            ksp;             /* linear solver context */
 25:   Mat            A,N;                /* matrix */
 26:   Vec            x,b,u,Ab;          /* approx solution, RHS, exact solution */
 27:   PetscViewer    fd;               /* viewer */
 28:   char           file[PETSC_MAX_PATH_LEN];     /* input file name */
 29:   PetscErrorCode ierr,ierrp;
 30:   PetscInt       its,n,m;
 31:   PetscReal      norm;

 33:   PetscInitialize(&argc,&args,(char*)0,help);


 36:   /*
 37:      Determine files from which we read the linear system
 38:      (matrix and right-hand-side vector).
 39:   */
 40:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);

 42:   /* -----------------------------------------------------------
 43:                   Beginning of linear solver loop
 44:      ----------------------------------------------------------- */
 45:   /*
 46:      Loop through the linear solve 2 times.
 47:       - The intention here is to preload and solve a small system;
 48:         then load another (larger) system and solve it as well.
 49:         This process preloads the instructions with the smaller
 50:         system so that more accurate performance monitoring (via
 51:         -log_summary) can be done with the larger one (that actually
 52:         is the system of interest).
 53:   */
 54:   PetscPreLoadBegin(PETSC_FALSE,"Load system");

 56:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 57:                          Load system
 58:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 60:   /*
 61:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 62:      reading from this file.
 63:   */
 64:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 66:   /*
 67:      Load the matrix and vector; then destroy the viewer.
 68:   */
 69:   MatCreate(PETSC_COMM_WORLD,&A);
 70:   MatSetType(A,MATMPIAIJ);
 71:   MatLoad(A,fd);
 72:   VecCreate(PETSC_COMM_WORLD,&b);
 73:   PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
 74:   ierrp = VecLoad(b,fd);
 75:   PetscPopErrorHandler();
 76:   MatGetLocalSize(A,&m,&n);
 77:   if (ierrp) {   /* if file contains no RHS, then use a vector of all ones */
 78:     PetscScalar one = 1.0;
 79:     VecCreate(PETSC_COMM_WORLD,&b);
 80:     VecSetSizes(b,m,PETSC_DECIDE);
 81:     VecSetFromOptions(b);
 82:     VecSet(b,one);
 83:   }
 84:   PetscViewerDestroy(&fd);

 86:   /*
 87:      If the loaded matrix is larger than the vector (due to being padded
 88:      to match the block size of the system), then create a new padded vector.
 89:   */
 90:   {
 91:     PetscInt    j,mvec,start,end,idex;
 92:     Vec         tmp;
 93:     PetscScalar *bold;

 95:     /* Create a new vector b by padding the old one */
 96:     VecCreate(PETSC_COMM_WORLD,&tmp);
 97:     VecSetSizes(tmp,m,PETSC_DECIDE);
 98:     VecSetFromOptions(tmp);
 99:     VecGetOwnershipRange(b,&start,&end);
100:     VecGetLocalSize(b,&mvec);
101:     VecGetArray(b,&bold);
102:     for (j=0; j<mvec; j++) {
103:       idex = start+j;
104:       VecSetValues(tmp,1,&idex,bold+j,INSERT_VALUES);
105:     }
106:     VecRestoreArray(b,&bold);
107:     VecDestroy(&b);
108:     VecAssemblyBegin(tmp);
109:     VecAssemblyEnd(tmp);
110:     b    = tmp;
111:   }
112:   VecDuplicate(b,&u);
113:   VecCreateMPI(PETSC_COMM_WORLD,n,PETSC_DECIDE,&x);
114:   VecDuplicate(x,&Ab);
115:   VecSet(x,0.0);

117:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
118:                     Setup solve for system
119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   /*
122:      Conclude profiling last stage; begin profiling next stage.
123:   */
124:   PetscPreLoadStage("KSPSetUp");

126:   MatCreateNormal(A,&N);
127:   MatMultTranspose(A,b,Ab);

129:   /*
130:      Create linear solver; set operators; set runtime options.
131:   */
132:   KSPCreate(PETSC_COMM_WORLD,&ksp);

134:   KSPSetOperators(ksp,N,N);
135:   KSPSetFromOptions(ksp);

137:   /*
138:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
139:      enable more precise profiling of setting up the preconditioner.
140:      These calls are optional, since both will be called within
141:      KSPSolve() if they haven't been called already.
142:   */
143:   KSPSetUp(ksp);
144:   KSPSetUpOnBlocks(ksp);

146:   /*
147:                          Solve system
148:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

150:   /*
151:      Begin profiling next stage
152:   */
153:   PetscPreLoadStage("KSPSolve");

155:   /*
156:      Solve linear system
157:   */
158:   KSPSolve(ksp,Ab,x);

160:   /*
161:       Conclude profiling this stage
162:    */
163:   PetscPreLoadStage("Cleanup");

165:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
166:           Check error, print output, free data structures.
167:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

169:   /*
170:      Check error
171:   */
172:   MatMult(A,x,u);
173:   VecAXPY(u,-1.0,b);
174:   VecNorm(u,NORM_2,&norm);
175:   KSPGetIterationNumber(ksp,&its);
176:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
177:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

179:   /*
180:      Free work space.  All PETSc objects should be destroyed when they
181:      are no longer needed.
182:   */
183:   MatDestroy(&A); VecDestroy(&b);
184:   MatDestroy(&N); VecDestroy(&Ab);
185:   VecDestroy(&u); VecDestroy(&x);
186:   KSPDestroy(&ksp);
187:   PetscPreLoadEnd();
188:   /* -----------------------------------------------------------
189:                       End of linear solver loop
190:      ----------------------------------------------------------- */

192:   PetscFinalize();
193:   return 0;
194: }