Actual source code: ex27.c
petsc-3.7.5 2017-01-01
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,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);
87: VecDuplicate(b,&u);
88: VecCreateMPI(PETSC_COMM_WORLD,n,PETSC_DECIDE,&x);
89: VecDuplicate(x,&Ab);
90: VecSet(x,0.0);
92: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
93: Setup solve for system
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: /*
97: Conclude profiling last stage; begin profiling next stage.
98: */
99: PetscPreLoadStage("KSPSetUp");
101: MatCreateNormal(A,&N);
102: MatMultTranspose(A,b,Ab);
104: /*
105: Create linear solver; set operators; set runtime options.
106: */
107: KSPCreate(PETSC_COMM_WORLD,&ksp);
109: KSPSetOperators(ksp,N,N);
110: KSPSetFromOptions(ksp);
112: /*
113: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
114: enable more precise profiling of setting up the preconditioner.
115: These calls are optional, since both will be called within
116: KSPSolve() if they haven't been called already.
117: */
118: KSPSetUp(ksp);
119: KSPSetUpOnBlocks(ksp);
121: /*
122: Solve system
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: /*
126: Begin profiling next stage
127: */
128: PetscPreLoadStage("KSPSolve");
130: /*
131: Solve linear system
132: */
133: KSPSolve(ksp,Ab,x);
135: /*
136: Conclude profiling this stage
137: */
138: PetscPreLoadStage("Cleanup");
140: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
141: Check error, print output, free data structures.
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: /*
145: Check error
146: */
147: MatMult(A,x,u);
148: VecAXPY(u,-1.0,b);
149: VecNorm(u,NORM_2,&norm);
150: KSPGetIterationNumber(ksp,&its);
151: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
152: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
154: /*
155: Free work space. All PETSc objects should be destroyed when they
156: are no longer needed.
157: */
158: MatDestroy(&A); VecDestroy(&b);
159: MatDestroy(&N); VecDestroy(&Ab);
160: VecDestroy(&u); VecDestroy(&x);
161: KSPDestroy(&ksp);
162: PetscPreLoadEnd();
163: /* -----------------------------------------------------------
164: End of linear solver loop
165: ----------------------------------------------------------- */
167: PetscFinalize();
168: return 0;
169: }