Actual source code: ex7.c
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: PetscTruth 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(PETSC_NULL,"-f0",file[0],PETSC_MAX_PATH_LEN-1,&flg);
43: if (!flg) SETERRQ(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: MatLoad(fd,MATSEQBAIJ,&A);
56: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&B);
57: VecLoad(fd,PETSC_NULL,&b);
58: PetscViewerDestroy(fd);
60: /*
61: If the loaded matrix is larger than the vector (due to being padded
62: to match the block size of the system), then create a new padded vector.
63: */
64: {
65: PetscInt m,n,j,mvec,start,end,idx;
66: Vec tmp;
67: PetscScalar *bold;
69: /* Create a new vector b by padding the old one */
70: MatGetLocalSize(A,&m,&n);
71: VecCreate(PETSC_COMM_WORLD,&tmp);
72: VecSetSizes(tmp,m,PETSC_DECIDE);
73: VecSetFromOptions(tmp);
74: VecGetOwnershipRange(b,&start,&end);
75: VecGetLocalSize(b,&mvec);
76: VecGetArray(b,&bold);
77: for (j=0; j<mvec; j++) {
78: idx = start+j;
79: VecSetValues(tmp,1,&idx,bold+j,INSERT_VALUES);
80: }
81: VecRestoreArray(b,&bold);
82: VecDestroy(b);
83: VecAssemblyBegin(tmp);
84: VecAssemblyEnd(tmp);
85: b = tmp;
86: }
87: VecDuplicate(b,&x);
88: VecDuplicate(b,&u);
89: VecSet(x,0.0);
91: /*
92: Create linear solver; set operators; set runtime options.
93: */
94: KSPCreate(PETSC_COMM_WORLD,&ksp);
95: KSPSetOperators(ksp,A,B,DIFFERENT_NONZERO_PATTERN);
96: KSPSetFromOptions(ksp);
98: /*
99: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
100: enable more precise profiling of setting up the preconditioner.
101: These calls are optional, since both will be called within
102: KSPSolve() if they haven't been called already.
103: */
104: KSPSetUp(ksp);
105: KSPSetUpOnBlocks(ksp);
106: KSPSolve(ksp,b,x);
108: /*
109: Check error, print output, free data structures.
110: This stage is not profiled separately.
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: /*
114: Check error
115: */
116: MatMult(A,x,u);
117: VecAXPY(u,-1.0,b);
118: VecNorm(u,NORM_2,&norm);
119: KSPGetIterationNumber(ksp,&its);
120: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
121: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %A\n",norm);
123: /*
124: Free work space. All PETSc objects should be destroyed when they
125: are no longer needed.
126: */
127: MatDestroy(A);
128: MatDestroy(B);
129: VecDestroy(b);
130: VecDestroy(u); VecDestroy(x);
131: KSPDestroy(ksp);
134: PetscFinalize();
135: return 0;
136: }