1 ! 2 subroutine dgtsl(n, c, d, e, b, info) 3 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 15:38:13 1/21/99 4 !...Switches: 5 implicit none 6 !----------------------------------------------- 7 ! D u m m y A r g u m e n t s 8 !----------------------------------------------- 9 integer n, info 10 double precision, dimension(n) :: c, d, e, b 11 !----------------------------------------------- 12 ! L o c a l V a r i a b l e s 13 !----------------------------------------------- 14 integer :: k, kb, kp1, nm1, nm2 15 double precision :: t,tc,td,te,tb 16 !----------------------------------------------- 17 ! 18 ! dgtsl given a general tridiagonal matrix and a right hand 19 ! side will find the solution. 20 ! 21 ! on entry 22 ! 23 ! n integer 24 ! is the order of the tridiagonal matrix. 25 ! 26 ! c double precision(n) 27 ! is the subdiagonal of the tridiagonal matrix. 28 ! c(2) through c(n) should contain the subdiagonal. 29 ! on output c is destroyed. 30 ! 31 ! d double precision(n) 32 ! is the diagonal of the tridiagonal matrix. 33 ! on output d is destroyed. 34 ! 35 ! e double precision(n) 36 ! is the superdiagonal of the tridiagonal matrix. 37 ! e(1) through e(n-1) should contain the superdiagonal. 38 ! on output e is destroyed. 39 ! 40 ! b double precision(n) 41 ! is the right hand side vector. 42 ! 43 ! on return 44 ! 45 ! b is the solution vector. 46 ! 47 ! info integer 48 ! = 0 normal value. 49 ! = k if the k-th element of the diagonal becomes 50 ! exactly zero. the subroutine returns when 51 ! this is detected. 52 ! 53 ! linpack. this version dated 08/14/78 . 54 ! jack dongarra, argonne national laboratory. 55 ! 56 ! no externals 57 ! fortran dabs 58 ! 59 ! internal variables 60 ! 61 ! begin block permitting ...exits to 100 62 ! 63 info = 0 64 c(1) = d(1) 65 nm1 = n - 1 66 if (nm1 >= 1) then 67 d(1) = e(1) 68 e(1) = 0.0D0 69 e(n) = 0.0D0 70 ! 71 do k = 1, nm1 72 kp1 = k + 1 73 ! 74 ! find the largest of the two rows 75 ! 76 if (dabs(c(kp1)) >= dabs(c(k))) then 77 ! 78 ! interchange row 79 ! 80 tc = c(kp1) 81 c(kp1) = c(k) 82 c(k) = tc 83 84 td = d(kp1) 85 d(kp1) = d(k) 86 d(k) = td 87 88 te = e(kp1) 89 e(kp1) = e(k) 90 e(k) = te 91 92 tb = b(kp1) 93 b(kp1) = b(k) 94 b(k) = tb 95 endif 96 if (c(k) == 0.0D0) then 97 info = k 98 ! ............exit 99 go to 100 100 endif 101 t = -c(kp1)/c(k) 102 c(kp1) = d(kp1) + t*d(k) 103 d(kp1) = e(kp1) + t*e(k) 104 e(kp1) = 0.0D0 105 b(kp1) = b(kp1) + t*b(k) 106 end do 107 endif 108 if (c(n) == 0.0D0) then 109 info = n 110 else 111 nm2 = n - 2 112 b(n) = b(n)/c(n) 113 if (n /= 1) then 114 b(nm1) = (b(nm1)-d(nm1)*b(n))/c(nm1) 115 if (nm2 >= 1) then 116 do kb = 1, nm2 117 k = nm2 - kb + 1 118 b(k) = (b(k)-d(k)*b(k+1)-e(k)*b(k+2))/c(k) 119 end do 120 endif 121 endif 122 endif 123 100 continue 124 return 125 end subroutine dgtsl 126