Parent Directory | Revision Log

Revision **1** -
(**show annotations**)
(**download**)

*Fri Oct 29 20:54:12 2004 UTC*
(17 years, 11 months ago)
by *aw0a*

File size: 5418 byte(s)

File size: 5418 byte(s)

Setting up web subdirectory in repository

1 | C dgbfa.f |

2 | C is freely available from netlib. It is not subject to any GNU License |

3 | C set by the authors of the ASCEND math programming system. |

4 | C $Date: 1996/04/30 18:17:11 $ $Revision: 1.1.1.1 $ |

5 | C |

6 | subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info) |

7 | integer lda,n,ml,mu,ipvt(1),info |

8 | double precision abd(lda,1) |

9 | c |

10 | c dgbfa factors a double precision band matrix by elimination. |

11 | c |

12 | c dgbfa is usually called by dgbco, but it can be called |

13 | c directly with a saving in time if rcond is not needed. |

14 | c |

15 | c on entry |

16 | c |

17 | c abd double precision(lda, n) |

18 | c contains the matrix in band storage. the columns |

19 | c of the matrix are stored in the columns of abd and |

20 | c the diagonals of the matrix are stored in rows |

21 | c ml+1 through 2*ml+mu+1 of abd . |

22 | c see the comments below for details. |

23 | c |

24 | c lda integer |

25 | c the leading dimension of the array abd . |

26 | c lda must be .ge. 2*ml + mu + 1 . |

27 | c |

28 | c n integer |

29 | c the order of the original matrix. |

30 | c |

31 | c ml integer |

32 | c number of diagonals below the main diagonal. |

33 | c 0 .le. ml .lt. n . |

34 | c |

35 | c mu integer |

36 | c number of diagonals above the main diagonal. |

37 | c 0 .le. mu .lt. n . |

38 | c more efficient if ml .le. mu . |

39 | c on return |

40 | c |

41 | c abd an upper triangular matrix in band storage and |

42 | c the multipliers which were used to obtain it. |

43 | c the factorization can be written a = l*u where |

44 | c l is a product of permutation and unit lower |

45 | c triangular matrices and u is upper triangular. |

46 | c |

47 | c ipvt integer(n) |

48 | c an integer vector of pivot indices. |

49 | c |

50 | c info integer |

51 | c = 0 normal value. |

52 | c = k if u(k,k) .eq. 0.0 . this is not an error |

53 | c condition for this subroutine, but it does |

54 | c indicate that dgbsl will divide by zero if |

55 | c called. use rcond in dgbco for a reliable |

56 | c indication of singularity. |

57 | c |

58 | c band storage |

59 | c |

60 | c if a is a band matrix, the following program segment |

61 | c will set up the input. |

62 | c |

63 | c ml = (band width below the diagonal) |

64 | c mu = (band width above the diagonal) |

65 | c m = ml + mu + 1 |

66 | c do 20 j = 1, n |

67 | c i1 = max0(1, j-mu) |

68 | c i2 = min0(n, j+ml) |

69 | c do 10 i = i1, i2 |

70 | c k = i - j + m |

71 | c abd(k,j) = a(i,j) |

72 | c 10 continue |

73 | c 20 continue |

74 | c |

75 | c this uses rows ml+1 through 2*ml+mu+1 of abd . |

76 | c in addition, the first ml rows in abd are used for |

77 | c elements generated during the triangularization. |

78 | c the total number of rows needed in abd is 2*ml+mu+1 . |

79 | c the ml+mu by ml+mu upper left triangle and the |

80 | c ml by ml lower right triangle are not referenced. |

81 | c |

82 | c linpack. this version dated 08/14/78 . |

83 | c cleve moler, university of new mexico, argonne national lab. |

84 | c |

85 | c subroutines and functions |

86 | c |

87 | c blas daxpy,dscal,idamax |

88 | c fortran max0,min0 |

89 | c |

90 | c internal variables |

91 | c |

92 | double precision t |

93 | integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1 |

94 | c |

95 | c |

96 | m = ml + mu + 1 |

97 | info = 0 |

98 | c |

99 | c zero initial fill-in columns |

100 | c |

101 | j0 = mu + 2 |

102 | j1 = min0(n,m) - 1 |

103 | if (j1 .lt. j0) go to 30 |

104 | do 20 jz = j0, j1 |

105 | i0 = m + 1 - jz |

106 | do 10 i = i0, ml |

107 | abd(i,jz) = 0.0d0 |

108 | 10 continue |

109 | 20 continue |

110 | 30 continue |

111 | jz = j1 |

112 | ju = 0 |

113 | c |

114 | c gaussian elimination with partial pivoting |

115 | c |

116 | nm1 = n - 1 |

117 | if (nm1 .lt. 1) go to 130 |

118 | do 120 k = 1, nm1 |

119 | kp1 = k + 1 |

120 | c |

121 | c zero next fill-in column |

122 | c |

123 | jz = jz + 1 |

124 | if (jz .gt. n) go to 50 |

125 | if (ml .lt. 1) go to 50 |

126 | do 40 i = 1, ml |

127 | abd(i,jz) = 0.0d0 |

128 | 40 continue |

129 | 50 continue |

130 | c |

131 | c find l = pivot index |

132 | c |

133 | lm = min0(ml,n-k) |

134 | l = idamax(lm+1,abd(m,k),1) + m - 1 |

135 | ipvt(k) = l + k - m |

136 | c |

137 | c zero pivot implies this column already triangularized |

138 | c |

139 | if (abd(l,k) .eq. 0.0d0) go to 100 |

140 | c |

141 | c interchange if necessary |

142 | c |

143 | if (l .eq. m) go to 60 |

144 | t = abd(l,k) |

145 | abd(l,k) = abd(m,k) |

146 | abd(m,k) = t |

147 | 60 continue |

148 | c |

149 | c compute multipliers |

150 | c |

151 | t = -1.0d0/abd(m,k) |

152 | call dscal(lm,t,abd(m+1,k),1) |

153 | c |

154 | c row elimination with column indexing |

155 | c |

156 | ju = min0(max0(ju,mu+ipvt(k)),n) |

157 | mm = m |

158 | if (ju .lt. kp1) go to 90 |

159 | do 80 j = kp1, ju |

160 | l = l - 1 |

161 | mm = mm - 1 |

162 | t = abd(l,j) |

163 | if (l .eq. mm) go to 70 |

164 | abd(l,j) = abd(mm,j) |

165 | abd(mm,j) = t |

166 | 70 continue |

167 | call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1) |

168 | 80 continue |

169 | 90 continue |

170 | go to 110 |

171 | 100 continue |

172 | info = k |

173 | 110 continue |

174 | 120 continue |

175 | 130 continue |

176 | ipvt(n) = n |

177 | if (abd(m,n) .eq. 0.0d0) info = n |

178 | return |

179 | end |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |