Actual source code: ex16f.F90
1: !
2: program main
3: #include <petsc/finclude/petscksp.h>
4: use petscksp
5: implicit none
7: !
8: ! This example is a modified Fortran version of ex6.c. It tests the use of
9: ! options prefixes in PETSc. Two linear problems are solved in this program.
10: ! The first problem is read from a file. The second problem is constructed
11: ! from the first, by eliminating some of the entries of the linear matrix 'A'.
13: ! Each solve is distinguished by a unique prefix - 'a' for the first, 'b'
14: ! for the second. With the prefix the user can distinguish between the various
15: ! options (command line, from .petscrc file, etc.) for each of the solvers.
16: ! Input arguments are:
17: ! -f <input_file> : file to load. For a 5X5 example of the 5-pt. stencil
18: ! use the file petsc/src/mat/examples/mat.ex.binary
20: PetscErrorCode ierr
21: PetscInt its,ione,ifive,izero
22: PetscBool flg
23: PetscScalar none,five
24: PetscReal norm
25: Vec x,b,u
26: Mat A
27: KSP ksp1,ksp2
28: character*(PETSC_MAX_PATH_LEN) f
29: PetscViewer fd
30: IS isrow
31: none = -1.0
32: five = 5.0
33: ifive = 5
34: ione = 1
35: izero = 0
37: PetscCallA(PetscInitialize(ierr))
39: ! Read in matrix and RHS
40: PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-f',f,flg,ierr))
41: PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD,f,FILE_MODE_READ,fd,ierr))
43: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
44: PetscCallA(MatSetType(A, MATSEQAIJ,ierr))
45: PetscCallA(MatLoad(A,fd,ierr))
47: PetscCallA(VecCreate(PETSC_COMM_WORLD,b,ierr))
48: PetscCallA(VecLoad(b,fd,ierr))
49: PetscCallA(PetscViewerDestroy(fd,ierr))
51: ! Set up solution
52: PetscCallA(VecDuplicate(b,x,ierr))
53: PetscCallA(VecDuplicate(b,u,ierr))
55: ! Solve system-1
56: PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp1,ierr))
57: PetscCallA(KSPSetOptionsPrefix(ksp1,'a',ierr))
58: PetscCallA(KSPAppendOptionsPrefix(ksp1,'_',ierr))
59: PetscCallA(KSPSetOperators(ksp1,A,A,ierr))
60: PetscCallA(KSPSetFromOptions(ksp1,ierr))
61: PetscCallA(KSPSolve(ksp1,b,x,ierr))
63: ! Show result
64: PetscCallA(MatMult(A,x,u,ierr))
65: PetscCallA(VecAXPY(u,none,b,ierr))
66: PetscCallA(VecNorm(u,NORM_2,norm,ierr))
67: PetscCallA(KSPGetIterationNumber(ksp1,its,ierr))
69: write(6,100) norm,its
70: 100 format('Residual norm ',e11.4,' iterations ',i5)
72: ! Create system 2 by striping off some rows of the matrix
73: PetscCallA(ISCreateStride(PETSC_COMM_SELF,ifive,izero,ione,isrow,ierr))
74: PetscCallA(MatZeroRowsIS(A,isrow,five,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr))
76: ! Solve system-2
77: PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp2,ierr))
78: PetscCallA(KSPSetOptionsPrefix(ksp2,'b',ierr))
79: PetscCallA(KSPAppendOptionsPrefix(ksp2,'_',ierr))
80: PetscCallA(KSPSetOperators(ksp2,A,A,ierr))
81: PetscCallA(KSPSetFromOptions(ksp2,ierr))
82: PetscCallA(KSPSolve(ksp2,b,x,ierr))
84: ! Show result
85: PetscCallA(MatMult(A,x,u,ierr))
86: PetscCallA(VecAXPY(u,none,b,ierr))
87: PetscCallA(VecNorm(u,NORM_2,norm,ierr))
88: PetscCallA(KSPGetIterationNumber(ksp2,its,ierr))
89: write(6,100) norm,its
91: ! Cleanup
92: PetscCallA(KSPDestroy(ksp1,ierr))
93: PetscCallA(KSPDestroy(ksp2,ierr))
94: PetscCallA(VecDestroy(b,ierr))
95: PetscCallA(VecDestroy(x,ierr))
96: PetscCallA(VecDestroy(u,ierr))
97: PetscCallA(MatDestroy(A,ierr))
98: PetscCallA(ISDestroy(isrow,ierr))
100: PetscCallA(PetscFinalize(ierr))
101: end
103: !/*TEST
104: !
105: ! test:
106: ! args: -f ${DATAFILESPATH}/matrices/arco1 -options_left no
107: ! requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
108: !
109: !TEST*/