Actual source code: ex7.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: Tests inplace factorization for SeqBAIJ. Input parameters include\n\
4: -f0 <input_file> : first file to load (small system)\n\n";
6: /*T
7: Concepts: KSP^solving a linear system
8: Concepts: PetscLog^profiling multiple stages of code;
9: Processors: n
10: T*/
12: /*
13: Include "petscksp.h" so that we can use KSP solvers. Note that this file
14: automatically includes:
15: petscsys.h - base PETSc routines petscvec.h - vectors
16: petscmat.h - matrices
17: petscis.h - index sets petscksp.h - Krylov subspace methods
18: petscviewer.h - viewers petscpc.h - preconditioners
19: */
20: #include <petscksp.h>
24: int main(int argc,char **args)
25: {
26: KSP ksp; /* linear solver context */
27: Mat A,B; /* matrix */
28: Vec x,b,u; /* approx solution, RHS, exact solution */
29: PetscViewer fd; /* viewer */
30: char file[2][PETSC_MAX_PATH_LEN]; /* input file name */
32: PetscInt its;
33: PetscBool flg;
34: PetscReal norm;
36: PetscInitialize(&argc,&args,(char*)0,help);
38: /*
39: Determine files from which we read the two linear systems
40: (matrix and right-hand-side vector).
41: */
42: PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
43: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 option");
46: /*
47: Open binary file. Note that we use FILE_MODE_READ to indicate
48: reading from this file.
49: */
50: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
52: /*
53: Load the matrix and vector; then destroy the viewer.
54: */
55: MatCreate(PETSC_COMM_WORLD,&A);
56: MatSetType(A,MATSEQBAIJ);
57: MatLoad(A,fd);
58: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&B);
59: VecCreate(PETSC_COMM_WORLD,&b);
60: VecLoad(b,fd);
61: PetscViewerDestroy(&fd);
63: /*
64: If the loaded matrix is larger than the vector (due to being padded
65: to match the block size of the system), then create a new padded vector.
66: */
67: {
68: PetscInt m,n,j,mvec,start,end,idx;
69: Vec tmp;
70: PetscScalar *bold;
72: /* Create a new vector b by padding the old one */
73: MatGetLocalSize(A,&m,&n);
74: VecCreate(PETSC_COMM_WORLD,&tmp);
75: VecSetSizes(tmp,m,PETSC_DECIDE);
76: VecSetFromOptions(tmp);
77: VecGetOwnershipRange(b,&start,&end);
78: VecGetLocalSize(b,&mvec);
79: VecGetArray(b,&bold);
80: for (j=0; j<mvec; j++) {
81: idx = start+j;
82: VecSetValues(tmp,1,&idx,bold+j,INSERT_VALUES);
83: }
84: VecRestoreArray(b,&bold);
85: VecDestroy(&b);
86: VecAssemblyBegin(tmp);
87: VecAssemblyEnd(tmp);
88: b = tmp;
89: }
90: VecDuplicate(b,&x);
91: VecDuplicate(b,&u);
92: VecSet(x,0.0);
94: /*
95: Create linear solver; set operators; set runtime options.
96: */
97: KSPCreate(PETSC_COMM_WORLD,&ksp);
98: KSPSetOperators(ksp,A,B,DIFFERENT_NONZERO_PATTERN);
99: KSPSetFromOptions(ksp);
101: /*
102: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
103: enable more precise profiling of setting up the preconditioner.
104: These calls are optional, since both will be called within
105: KSPSolve() if they haven't been called already.
106: */
107: KSPSetUp(ksp);
108: KSPSetUpOnBlocks(ksp);
109: KSPSolve(ksp,b,x);
111: /*
112: Check error, print output, free data structures.
113: This stage is not profiled separately.
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: /*
117: Check error
118: */
119: MatMult(A,x,u);
120: VecAXPY(u,-1.0,b);
121: VecNorm(u,NORM_2,&norm);
122: KSPGetIterationNumber(ksp,&its);
123: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
124: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %G\n",norm);
126: /*
127: Free work space. All PETSc objects should be destroyed when they
128: are no longer needed.
129: */
130: MatDestroy(&A);
131: MatDestroy(&B);
132: VecDestroy(&b);
133: VecDestroy(&u); VecDestroy(&x);
134: KSPDestroy(&ksp);
137: PetscFinalize();
138: return 0;
139: }