Actual source code: ex27.c
petsc-3.5.4 2015-05-23
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(NULL,"-f",file,PETSC_MAX_PATH_LEN,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: PetscPreLoadBegin(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: MatCreate(PETSC_COMM_WORLD,&A);
70: MatSetType(A,MATMPIAIJ);
71: MatLoad(A,fd);
72: VecCreate(PETSC_COMM_WORLD,&b);
73: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
74: ierrp = VecLoad(b,fd);
75: PetscPopErrorHandler();
76: MatGetLocalSize(A,&m,&n);
77: if (ierrp) { /* if file contains no RHS, then use a vector of all ones */
78: PetscScalar one = 1.0;
79: VecCreate(PETSC_COMM_WORLD,&b);
80: VecSetSizes(b,m,PETSC_DECIDE);
81: VecSetFromOptions(b);
82: VecSet(b,one);
83: }
84: PetscViewerDestroy(&fd);
86: /*
87: If the loaded matrix is larger than the vector (due to being padded
88: to match the block size of the system), then create a new padded vector.
89: */
90: {
91: PetscInt j,mvec,start,end,idex;
92: Vec tmp;
93: PetscScalar *bold;
95: /* Create a new vector b by padding the old one */
96: VecCreate(PETSC_COMM_WORLD,&tmp);
97: VecSetSizes(tmp,m,PETSC_DECIDE);
98: VecSetFromOptions(tmp);
99: VecGetOwnershipRange(b,&start,&end);
100: VecGetLocalSize(b,&mvec);
101: VecGetArray(b,&bold);
102: for (j=0; j<mvec; j++) {
103: idex = start+j;
104: VecSetValues(tmp,1,&idex,bold+j,INSERT_VALUES);
105: }
106: VecRestoreArray(b,&bold);
107: VecDestroy(&b);
108: VecAssemblyBegin(tmp);
109: VecAssemblyEnd(tmp);
110: b = tmp;
111: }
112: VecDuplicate(b,&u);
113: VecCreateMPI(PETSC_COMM_WORLD,n,PETSC_DECIDE,&x);
114: VecDuplicate(x,&Ab);
115: VecSet(x,0.0);
117: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
118: Setup solve for system
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: /*
122: Conclude profiling last stage; begin profiling next stage.
123: */
124: PetscPreLoadStage("KSPSetUp");
126: MatCreateNormal(A,&N);
127: MatMultTranspose(A,b,Ab);
129: /*
130: Create linear solver; set operators; set runtime options.
131: */
132: KSPCreate(PETSC_COMM_WORLD,&ksp);
134: KSPSetOperators(ksp,N,N);
135: KSPSetFromOptions(ksp);
137: /*
138: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
139: enable more precise profiling of setting up the preconditioner.
140: These calls are optional, since both will be called within
141: KSPSolve() if they haven't been called already.
142: */
143: KSPSetUp(ksp);
144: KSPSetUpOnBlocks(ksp);
146: /*
147: Solve system
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: /*
151: Begin profiling next stage
152: */
153: PetscPreLoadStage("KSPSolve");
155: /*
156: Solve linear system
157: */
158: KSPSolve(ksp,Ab,x);
160: /*
161: Conclude profiling this stage
162: */
163: PetscPreLoadStage("Cleanup");
165: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
166: Check error, print output, free data structures.
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: /*
170: Check error
171: */
172: MatMult(A,x,u);
173: VecAXPY(u,-1.0,b);
174: VecNorm(u,NORM_2,&norm);
175: KSPGetIterationNumber(ksp,&its);
176: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
177: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
179: /*
180: Free work space. All PETSc objects should be destroyed when they
181: are no longer needed.
182: */
183: MatDestroy(&A); VecDestroy(&b);
184: MatDestroy(&N); VecDestroy(&Ab);
185: VecDestroy(&u); VecDestroy(&x);
186: KSPDestroy(&ksp);
187: PetscPreLoadEnd();
188: /* -----------------------------------------------------------
189: End of linear solver loop
190: ----------------------------------------------------------- */
192: PetscFinalize();
193: return 0;
194: }