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: }