1 | !> \file tridiag_mod-0.3.f90 |
---|
2 | !! Module pour la resolution d'un systeme tri-diagonal |
---|
3 | !< |
---|
4 | |
---|
5 | !> \namespace TRIDIAGMOD |
---|
6 | !! Module pour la resolution d'un systeme tri-diagonal |
---|
7 | !! \author ... |
---|
8 | !! \date ... |
---|
9 | !< |
---|
10 | MODULE TRIDIAGMOD |
---|
11 | |
---|
12 | CONTAINS |
---|
13 | |
---|
14 | !*************************************************************** |
---|
15 | ! |
---|
16 | ! tridiag : routine de resolution d'un systeme tri-diagonal |
---|
17 | ! M * U = R ou M est tridiag. (ssdiag :A, diag:B, surdiag : C) |
---|
18 | ! U est la solution recherchee |
---|
19 | ! |
---|
20 | !****************************************************************** |
---|
21 | |
---|
22 | !> SUBROUTINE: TRIDIAG() |
---|
23 | !! \author ... |
---|
24 | !! \date ... |
---|
25 | !! @note Cette routine permet de resoudre d'un systeme tri-diagonal |
---|
26 | !! @note M * U = R ou M est tridiag. (ssdiag :A, diag:B, surdiag : C) |
---|
27 | !! @note U est la solution recherchee |
---|
28 | !! \param A sous-diagonale l'indice 1 ne sert pas |
---|
29 | !! \param B diagonale |
---|
30 | !! \param C sur-diagonale, l'indice n ne sert pas |
---|
31 | !! \param R vecteur membre de droite |
---|
32 | !! \param U solution |
---|
33 | !! \param n |
---|
34 | !! \param iFAIL permet de detecter une erreur |
---|
35 | !< |
---|
36 | |
---|
37 | |
---|
38 | subroutine TRIDIAG (A,B,C,R,U,n,iFAIL) |
---|
39 | |
---|
40 | |
---|
41 | IMPLICIT NONE |
---|
42 | |
---|
43 | INTEGER, intent(in) :: n |
---|
44 | REAL,dimension(n), intent(in) :: a !< sous-diagonale l'indice 1 ne sert pas |
---|
45 | REAL,dimension(n), intent(inout) :: b !< diagonale |
---|
46 | REAL,dimension(n), intent(in) :: c !< sur-diagonale, l'indice n ne sert pas |
---|
47 | REAL,dimension(n), intent(in) :: r !< vecteur membre de droite |
---|
48 | REAL,dimension(n), intent(out) :: u !< solution |
---|
49 | REAL,dimension(n) :: gam !< tableau de travail. |
---|
50 | REAL :: BET |
---|
51 | INTEGER,intent(out) :: iFAIL !< permet de detecter une erreur |
---|
52 | INTEGER :: JJ |
---|
53 | |
---|
54 | |
---|
55 | |
---|
56 | |
---|
57 | iFAIL=0 |
---|
58 | |
---|
59 | if (abs(B(1)).lt.1.e-20) then |
---|
60 | B(1)=1.0 |
---|
61 | iFAIL=1 |
---|
62 | endif |
---|
63 | |
---|
64 | |
---|
65 | BET=1./B(1) |
---|
66 | U(1)=R(1)*BET |
---|
67 | |
---|
68 | |
---|
69 | !BET=B(1) |
---|
70 | !U(1)=R(1)/BET |
---|
71 | |
---|
72 | !GAM(:)=EOSHIFT(C,shift=-1) |
---|
73 | |
---|
74 | do jj=2,n |
---|
75 | GAM(jj)=C(jj-1)*BET |
---|
76 | BET=1./(B(jj)-A(jj)*GAM(jj)) |
---|
77 | |
---|
78 | ! GAM(JJ)=C(JJ-1)/BET |
---|
79 | ! BET=B(JJ)-A(JJ)*GAM(JJ) |
---|
80 | |
---|
81 | |
---|
82 | if (abs(BET)>1.e20) then |
---|
83 | ! if (abs(BET)<1.e-20) then |
---|
84 | BET=1.0 |
---|
85 | iFAIL=1 |
---|
86 | endif |
---|
87 | |
---|
88 | U(jj) = (R(JJ)-A(JJ)*U(JJ-1))*BET |
---|
89 | ! U(jj) = (R(JJ)-A(JJ)*U(JJ-1))/BET |
---|
90 | end do |
---|
91 | |
---|
92 | do JJ=N-1,1,-1 |
---|
93 | U(JJ)=U(JJ)-GAM(JJ+1)*U(JJ+1) |
---|
94 | end do |
---|
95 | |
---|
96 | if (iFAIL.eq.1) then |
---|
97 | write(6,*) 'a B element is nul. TRIDIAG cannot work.' |
---|
98 | endif |
---|
99 | |
---|
100 | end subroutine TRIDIAG |
---|
101 | |
---|
102 | end module tridiagmod |
---|