Actual source code: ex16f.F90
petsc-3.13.4 2020-08-01
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*(128) 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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
38: if (ierr .ne. 0) then
39: print*,'Unable to initialize PETSc'
40: stop
41: endif
43: ! Read in matrix and RHS
44: call PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-f',f,flg,ierr);CHKERRA(ierr)
45: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,f,FILE_MODE_READ,fd,ierr);CHKERRA(ierr)
47: call MatCreate(PETSC_COMM_WORLD,A,ierr)
48: call MatSetType(A, MATSEQAIJ,ierr)
49: call MatLoad(A,fd,ierr)
51: call VecCreate(PETSC_COMM_WORLD,b,ierr)
52: call VecLoad(b,fd,ierr)
53: call PetscViewerDestroy(fd,ierr)
55: ! Set up solution
56: call VecDuplicate(b,x,ierr)
57: call VecDuplicate(b,u,ierr)
59: ! Solve system-1
60: call KSPCreate(PETSC_COMM_WORLD,ksp1,ierr)
61: call KSPSetOptionsPrefix(ksp1,'a',ierr)
62: call KSPAppendOptionsPrefix(ksp1,'_',ierr)
63: call KSPSetOperators(ksp1,A,A,ierr)
64: call KSPSetFromOptions(ksp1,ierr)
65: call KSPSolve(ksp1,b,x,ierr)
67: ! Show result
68: call MatMult(A,x,u,ierr)
69: call VecAXPY(u,none,b,ierr)
70: call VecNorm(u,NORM_2,norm,ierr)
71: call KSPGetIterationNumber(ksp1,its,ierr)
74: write(6,100) norm,its
75: 100 format('Residual norm ',e11.4,' iterations ',i5)
77: ! Create system 2 by striping off some rows of the matrix
78: call ISCreateStride(PETSC_COMM_SELF,ifive,izero,ione,isrow,ierr)
79: call MatZeroRowsIS(A,isrow,five,PETSC_NULL_VEC, &
80: & PETSC_NULL_VEC,ierr)
82: ! Solve system-2
83: call KSPCreate(PETSC_COMM_WORLD,ksp2,ierr)
84: call KSPSetOptionsPrefix(ksp2,'b',ierr)
85: call KSPAppendOptionsPrefix(ksp2,'_',ierr)
86: call KSPSetOperators(ksp2,A,A,ierr)
87: call KSPSetFromOptions(ksp2,ierr)
88: call KSPSolve(ksp2,b,x,ierr)
90: ! Show result
91: call MatMult(A,x,u,ierr)
92: call VecAXPY(u,none,b,ierr)
93: call VecNorm(u,NORM_2,norm,ierr)
94: call KSPGetIterationNumber(ksp2,its,ierr)
95: write(6,100) norm,its
97: ! Cleanup
98: call KSPDestroy(ksp1,ierr)
99: call KSPDestroy(ksp2,ierr)
100: call VecDestroy(b,ierr)
101: call VecDestroy(x,ierr)
102: call VecDestroy(u,ierr)
103: call MatDestroy(A,ierr)
104: call ISDestroy(isrow,ierr)
106: call PetscFinalize(ierr)
107: end
109: !/*TEST
110: !
111: ! test:
112: ! args: -f ${DATAFILESPATH}/matrices/arco1 -options_left no
113: ! requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
114: !
115: !TEST*/