Actual source code: ex27.c

  1: /*$Id: ex10.c,v 1.58 2001/09/11 16:33:29 bsmith Exp $*/

  3: static char help[] = "Reads a PETSc matrix and vector from a file and solves the normal equations.\n\n";
  4: /*T
  5:    Concepts: KSP^solving a linear system
  6:    Concepts: Normal equations
  7:    Processors: n
  8: T*/

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


 23: int main(int argc,char **args)
 24: {
 25:   KSP           ksp;             /* linear solver context */
 26:   Mat            A,N;                /* matrix */
 27:   Vec            x,b,u,Ab;          /* approx solution, RHS, exact solution */
 28:   PetscViewer    fd;               /* viewer */
 29:   char           file[128];     /* input file name */
 30:   int            ierr,its,ierrp,n,m;
 31:   PetscReal      norm;
 32:   PetscScalar    zero = 0.0,none = -1.0;
 33:   KSP            ksp;

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


 38:   /* 
 39:      Determine files from which we read the linear system
 40:      (matrix and right-hand-side vector).
 41:   */
 42:   PetscOptionsGetString(PETSC_NULL,"-f",file,127,PETSC_NULL);

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

 58:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 59:                            Load system
 60:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

133:     KSPSetOperators(ksp,N,N,SAME_NONZERO_PATTERN);
134:     KSPSetFromOptions(ksp);

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

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

151:     /*
152:        Begin profiling next stage
153:     */
154:     PreLoadStage("KSPSolve");

156:     /*
157:        Solve linear system
158:     */
159:     KSPSolve(ksp);

161:    /* 
162:        Conclude profiling this stage
163:     */
164:     PreLoadStage("Cleanup");

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

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

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

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