Actual source code: ex16f.F
1: !
2: ! "$Id: ex16.F,v 1.29 2001/08/07 03:03:57 balay Exp $";
3: !
4: program main
5: implicit none
7: #include include/finclude/petsc.h
8: #include include/finclude/petscvec.h
9: #include include/finclude/petscmat.h
10: #include include/finclude/petscpc.h
11: #include include/finclude/petscksp.h
12: #include include/finclude/petscviewer.h
13: #include include/finclude/petscis.h
14: !
15: ! This example is a modified Fortran version of ex6.c. It tests the use of
16: ! options prefixes in PETSc. Two linear problems are solved in this program.
17: ! The first problem is read from a file. The second problem is constructed
18: ! from the first, by eliminating some of the entries of the linear matrix 'A'.
20: ! Each solve is distinguished by a unique prefix - 'a' for the first, 'b'
21: ! for the second. With the prefix the user can distinguish between the various
22: ! options (command line, from .petscrc file, etc.) for each of the solvers.
23: ! Input arguments are:
24: ! -f <input_file> : file to load. For a 5X5 example of the 5-pt. stencil
25: ! use the file petsc/src/mat/examples/mat.ex.binary
27: integer ierr,its,flg
28: PetscScalar norm,none,five
29: Vec x,b,u
30: Mat A
31: KSP ksp1,ksp2
32: character*(128) f
33: PetscViewer fd
34: IS isrow
35: PetscLogDouble time1,time2
36: KSP ksp
37: none = -1.0
38: five = 5.0
39: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
41: ! Read in matrix and RHS
42: call PetscOptionsGetString(PETSC_NULL_CHARACTER,'-f',f,flg,ierr)
43: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,f,PETSC_FILE_RDONLY, &
44: & fd,ierr)
45: if (ierr .ne. 0) then
46: print*, 'Unable to open file ',f
47: SETERRQ(1,' ',ierr)
48: endif
50: call MatLoad(fd,MATSEQAIJ,A,ierr)
51: if (ierr .ne. 0) then
52: print*, 'Unable to load matrix '
53: SETERRQ(1,' ',ierr)
54: endif
56: call VecLoad(fd,PETSC_NULL_CHARACTER,b,ierr)
57: call PetscViewerDestroy(fd,ierr)
59: ! Set up solution
60: call VecDuplicate(b,x,ierr)
61: call VecDuplicate(b,u,ierr)
63: ! Solve system-1
64: call KSPCreate(PETSC_COMM_WORLD,ksp1,ierr)
65: call KSPSetOptionsPrefix(ksp1,'a',ierr)
66: call KSPAppendOptionsPrefix(ksp1,'_',ierr)
67: call KSPSetOperators(ksp1,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
68: call KSPSetFromOptions(ksp1,ierr)
69: call PetscGetTime(time1,ierr)
70: call KSPSetRhs(ksp1,b,ierr)
71: call KSPSetSolution(ksp1,x,ierr)
72: call KSPSolve(ksp1,ierr)
73: call PetscGetTime(time2,ierr)
75: ! Show result
76: call MatMult(A,x,u,ierr)
77: call VecAXPY(none,b,u,ierr)
78: call VecNorm(u,NORM_2,norm,ierr)
79: call KSPGetIterationNumber(ksp,its,ierr)
81: ! print*, 'Time for solve = ',time2-time1
83: write(6,100) norm,its
84: 100 format('Residual norm ',e10.4,' iterations ',i5)
86: ! Create system 2 by striping off some rows of the matrix
87: call ISCreateStride(PETSC_COMM_SELF,5,0,1,isrow,ierr)
88: call MatZeroRows(A,isrow,five,ierr)
90: ! Solve system-2
91: call KSPCreate(PETSC_COMM_WORLD,ksp2,ierr)
92: call KSPSetOptionsPrefix(ksp2,'b',ierr)
93: call KSPAppendOptionsPrefix(ksp2,'_',ierr)
94: call KSPSetOperators(ksp2,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
95: call KSPSetFromOptions(ksp2,ierr)
96: call PetscGetTime(time1,ierr)
97: call KSPSetRhs(ksp1,b,ierr)
98: call KSPSetSolution(ksp1,x,ierr)
99: call KSPSolve(ksp2,ierr)
100: call PetscGetTime(time2,ierr)
102: ! Show result
103: call MatMult(A,x,u,ierr)
104: call VecAXPY(none,b,u,ierr)
105: call VecNorm(u,NORM_2,norm,ierr)
106: call KSPGetIterationNumber(ksp,its,ierr)
107: ! print*, 'Time for solve = ',time2-time1
108: write(6,100) norm,its
110: ! Cleanup
111: call KSPDestroy(ksp1,ierr)
112: call KSPDestroy(ksp2,ierr)
113: call VecDestroy(b,ierr)
114: call VecDestroy(x,ierr)
115: call VecDestroy(u,ierr)
116: call MatDestroy(A,ierr)
117: call ISDestroy(isrow,ierr)
119: call PetscFinalize(ierr)
120: end