Actual source code: ex27.c

  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:      petsc.h       - base PETSc routines   petscvec.h - vectors
 13:      petscsys.h    - system routines       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(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,PETSC_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:   PreLoadBegin(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:     MatLoad(fd,MATMPIAIJ,&A);
 70:     PetscPushErrorHandler(PetscIgnoreErrorHandler,PETSC_NULL);
 71:     ierrp = VecLoad(fd,PETSC_NULL,&b);
 72:     PetscPopErrorHandler();
 73:     MatGetLocalSize(A,&m,&n);
 74:     if (ierrp) { /* if file contains no RHS, then use a vector of all ones */
 75:       PetscScalar one = 1.0;
 76:       VecCreate(PETSC_COMM_WORLD,&b);
 77:       VecSetSizes(b,m,PETSC_DECIDE);
 78:       VecSetFromOptions(b);
 79:       VecSet(b,one);
 80:     }
 81:     PetscViewerDestroy(fd);

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

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

114:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
115:                       Setup solve for system
116:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

118:     /*
119:        Conclude profiling last stage; begin profiling next stage.
120:     */
121:     PreLoadStage("KSPSetUp");

123:     MatCreateNormal(A,&N);
124:     MatMultTranspose(A,b,Ab);

126:     /*
127:        Create linear solver; set operators; set runtime options.
128:     */
129:     KSPCreate(PETSC_COMM_WORLD,&ksp);

131:     KSPSetOperators(ksp,N,N,SAME_NONZERO_PATTERN);
132:     KSPSetFromOptions(ksp);

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

143:     /*
144:                            Solve system
145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

147:     /*
148:        Begin profiling next stage
149:     */
150:     PreLoadStage("KSPSolve");

152:     /*
153:        Solve linear system
154:     */
155:     KSPSolve(ksp,Ab,x);

157:    /* 
158:        Conclude profiling this stage
159:     */
160:     PreLoadStage("Cleanup");

162:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
163:             Check error, print output, free data structures.
164:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

189:   PetscFinalize();
190:   return 0;
191: }