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: petscsys.h - base PETSc routines petscvec.h - vectors
13: 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: }