2
SUBROUTINE SGBCO (ABD, LDA, N, ML, MU, IPVT, RCOND, Z)
3
C***BEGIN PROLOGUE SGBCO
4
C***PURPOSE Factor a band matrix by Gaussian elimination and
5
C estimate the condition number of the matrix.
6
C***LIBRARY SLATEC (LINPACK)
8
C***TYPE SINGLE PRECISION (SGBCO-S, DGBCO-D, CGBCO-C)
9
C***KEYWORDS BANDED, CONDITION NUMBER, LINEAR ALGEBRA, LINPACK,
10
C MATRIX FACTORIZATION
11
C***AUTHOR Moler, C. B., (U. of New Mexico)
14
C SBGCO factors a real band matrix by Gaussian
15
C elimination and estimates the condition of the matrix.
17
C If RCOND is not needed, SGBFA is slightly faster.
18
C To solve A*X = B , follow SBGCO by SGBSL.
19
C To compute INVERSE(A)*C , follow SBGCO by SGBSL.
20
C To compute DETERMINANT(A) , follow SBGCO by SGBDI.
25
C contains the matrix in band storage. The columns
26
C of the matrix are stored in the columns of ABD and
27
C the diagonals of the matrix are stored in rows
28
C ML+1 through 2*ML+MU+1 of ABD .
29
C See the comments below for details.
32
C the leading dimension of the array ABD .
33
C LDA must be .GE. 2*ML + MU + 1 .
36
C the order of the original matrix.
39
C number of diagonals below the main diagonal.
43
C number of diagonals above the main diagonal.
45
C More efficient if ML .LE. MU .
49
C ABD an upper triangular matrix in band storage and
50
C the multipliers which were used to obtain it.
51
C The factorization can be written A = L*U where
52
C L is a product of permutation and unit lower
53
C triangular matrices and U is upper triangular.
56
C an integer vector of pivot indices.
59
C an estimate of the reciprocal condition of A .
60
C For the system A*X = B , relative perturbations
61
C in A and B of size EPSILON may cause
62
C relative perturbations in X of size EPSILON/RCOND .
63
C If RCOND is so small that the logical expression
64
C 1.0 + RCOND .EQ. 1.0
65
C is true, then A may be singular to working
66
C precision. In particular, RCOND is zero if
67
C exact singularity is detected or the estimate
71
C a work vector whose contents are usually unimportant.
72
C If A is close to a singular matrix, then Z is
73
C an approximate null vector in the sense that
74
C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
78
C If A is a band matrix, the following program segment
79
C will set up the input.
81
C ML = (band width below the diagonal)
82
C MU = (band width above the diagonal)
93
C This uses rows ML+1 through 2*ML+MU+1 of ABD .
94
C In addition, the first ML rows in ABD are used for
95
C elements generated during the triangularization.
96
C The total number of rows needed in ABD is 2*ML+MU+1 .
97
C The ML+MU by ML+MU upper left triangle and the
98
C ML by ML lower right triangle are not referenced.
100
C Example: If the original matrix is
109
C then N = 6, ML = 1, MU = 2, LDA .GE. 5 and ABD should contain
111
C * * * + + + , * = not used
112
C * * 13 24 35 46 , + = used for pivoting
117
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
118
C Stewart, LINPACK Users' Guide, SIAM, 1979.
119
C***ROUTINES CALLED SASUM, SAXPY, SDOT, SGBFA, SSCAL
120
C***REVISION HISTORY (YYMMDD)
121
C 780814 DATE WRITTEN
122
C 890531 Changed all specific intrinsics to generic. (WRB)
123
C 890831 Modified array declarations. (WRB)
124
C 890831 REVISION DATE from Version 3.2
125
C 891214 Prologue converted to Version 4.0 format. (BAB)
126
C 900326 Removed duplicate information from DESCRIPTION section.
128
C 920501 Reformatted the REFERENCES section. (WRB)
129
C***END PROLOGUE SGBCO
130
INTEGER LDA,N,ML,MU,IPVT(*)
134
REAL SDOT,EK,T,WK,WKM
135
REAL ANORM,S,SASUM,SM,YNORM
136
INTEGER IS,INFO,J,JU,K,KB,KP1,L,LA,LM,LZ,M,MM
138
C COMPUTE 1-NORM OF A
140
C***FIRST EXECUTABLE STATEMENT SGBCO
145
ANORM = MAX(ANORM,SASUM(L,ABD(IS,J),1))
146
IF (IS .GT. ML + 1) IS = IS - 1
147
IF (J .LE. MU) L = L + 1
148
IF (J .GE. N - ML) L = L - 1
153
CALL SGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
155
C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
156
C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E .
157
C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE
158
C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE
159
C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID
162
C SOLVE TRANS(U)*W = E
171
IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K))
172
IF (ABS(EK-Z(K)) .LE. ABS(ABD(M,K))) GO TO 30
173
S = ABS(ABD(M,K))/ABS(EK-Z(K))
181
IF (ABD(M,K) .EQ. 0.0E0) GO TO 40
190
JU = MIN(MAX(JU,MU+IPVT(K)),N)
192
IF (KP1 .GT. JU) GO TO 90
195
SM = SM + ABS(Z(J)+WKM*ABD(MM,J))
196
Z(J) = Z(J) + WK*ABD(MM,J)
199
IF (S .GE. SM) GO TO 80
205
Z(J) = Z(J) + T*ABD(MM,J)
211
S = 1.0E0/SASUM(N,Z,1)
214
C SOLVE TRANS(L)*Y = W
219
IF (K .LT. N) Z(K) = Z(K) + SDOT(LM,ABD(M+1,K),1,Z(K+1),1)
220
IF (ABS(Z(K)) .LE. 1.0E0) GO TO 110
229
S = 1.0E0/SASUM(N,Z,1)
242
IF (K .LT. N) CALL SAXPY(LM,T,ABD(M+1,K),1,Z(K+1),1)
243
IF (ABS(Z(K)) .LE. 1.0E0) GO TO 130
249
S = 1.0E0/SASUM(N,Z,1)
257
IF (ABS(Z(K)) .LE. ABS(ABD(M,K))) GO TO 150
258
S = ABS(ABD(M,K))/ABS(Z(K))
262
IF (ABD(M,K) .NE. 0.0E0) Z(K) = Z(K)/ABD(M,K)
263
IF (ABD(M,K) .EQ. 0.0E0) Z(K) = 1.0E0
268
CALL SAXPY(LM,T,ABD(LA,K),1,Z(LZ),1)
271
S = 1.0E0/SASUM(N,Z,1)
275
IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
276
IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
added
removed
deps/openlibm/slatec/sgbco.f
