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>


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

 31:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 32:   /*
 33:      Determine files from which we read the linear system
 34:      (matrix and right-hand-side vector).
 35:   */
 36:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);

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

 52:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 53:                          Load system
 54:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 56:   /*
 57:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 58:      reading from this file.
 59:   */
 60:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

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


 83:   VecDuplicate(b,&u);
 84:   VecCreateMPI(PETSC_COMM_WORLD,n,PETSC_DECIDE,&x);
 85:   VecDuplicate(x,&Ab);
 86:   VecSet(x,0.0);

 88:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 89:                     Setup solve for system
 90:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 92:   /*
 93:      Conclude profiling last stage; begin profiling next stage.
 94:   */
 95:   PetscPreLoadStage("KSPSetUp");

 97:   MatCreateNormal(A,&N);
 98:   MatMultTranspose(A,b,Ab);

100:   /*
101:      Create linear solver; set operators; set runtime options.
102:   */
103:   KSPCreate(PETSC_COMM_WORLD,&ksp);

105:   KSPSetOperators(ksp,N,N);
106:   KSPSetFromOptions(ksp);

108:   /*
109:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
110:      enable more precise profiling of setting up the preconditioner.
111:      These calls are optional, since both will be called within
112:      KSPSolve() if they haven't been called already.
113:   */
114:   KSPSetUp(ksp);
115:   KSPSetUpOnBlocks(ksp);

117:   /*
118:                          Solve system
119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   /*
122:      Begin profiling next stage
123:   */
124:   PetscPreLoadStage("KSPSolve");

126:   /*
127:      Solve linear system
128:   */
129:   KSPSolve(ksp,Ab,x);

131:   /*
132:       Conclude profiling this stage
133:    */
134:   PetscPreLoadStage("Cleanup");

136:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
137:           Check error, print output, free data structures.
138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

140:   /*
141:      Check error
142:   */
143:   MatMult(A,x,u);
144:   VecAXPY(u,-1.0,b);
145:   VecNorm(u,NORM_2,&norm);
146:   KSPGetIterationNumber(ksp,&its);
147:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
148:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

150:   /*
151:      Free work space.  All PETSc objects should be destroyed when they
152:      are no longer needed.
153:   */
154:   MatDestroy(&A); VecDestroy(&b);
155:   MatDestroy(&N); VecDestroy(&Ab);
156:   VecDestroy(&u); VecDestroy(&x);
157:   KSPDestroy(&ksp);
158:   PetscPreLoadEnd();
159:   /* -----------------------------------------------------------
160:                       End of linear solver loop
161:      ----------------------------------------------------------- */

163:   PetscFinalize();
164:   return ierr;
165: }