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