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.nn";
  4: /*T
  5:    Concepts: SLES^solving a linear system
  6:    Concepts: Normal equations
  7:    Processors: n
  8: T*/

 10: /* 
 11:   Include "petscsles.h" so that we can use SLES 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 petscsles.h


 21: #undef __FUNCT__
 23: int main(int argc,char **args)
 24: {
 25:   SLES           sles;             /* 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;

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


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

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

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

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

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

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

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

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

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

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

127:     /*
128:        Create linear solver; set operators; set runtime options.
129:     */
130:     SLESCreate(PETSC_COMM_WORLD,&sles);

132:     SLESSetOperators(sles,N,N,SAME_NONZERO_PATTERN);
133:     SLESSetFromOptions(sles);

135:     /* 
136:        Here we explicitly call SLESSetUp() and SLESSetUpOnBlocks() to
137:        enable more precise profiling of setting up the preconditioner.
138:        These calls are optional, since both will be called within
139:        SLESSolve() if they haven't been called already.
140:     */
141:     SLESSetUp(sles,Ab,x);
142:     SLESSetUpOnBlocks(sles);

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

148:     /*
149:        Begin profiling next stage
150:     */
151:     PreLoadStage("SLESSolve");

153:     /*
154:        Solve linear system
155:     */
156:     SLESSolve(sles,Ab,x,&its);

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

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

167:     /* 
168:        Check error
169:     */
170:     MatMult(A,x,u);
171:     VecAXPY(&none,b,u);
172:     VecNorm(u,NORM_2,&norm);

174:     PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3dn",its);
175:     PetscPrintf(PETSC_COMM_WORLD,"Residual norm %An",norm);

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

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