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