2
SUBROUTINE SNBIR (ABE, LDA, N, ML, MU, V, ITASK, IND, WORK, IWORK)
3
C***BEGIN PROLOGUE SNBIR
4
C***PURPOSE Solve a general nonsymmetric banded system of linear
5
C equations. Iterative refinement is used to obtain an error
9
C***TYPE SINGLE PRECISION (SNBIR-S, CNBIR-C)
10
C***KEYWORDS BANDED, LINEAR EQUATIONS, NONSYMMETRIC
11
C***AUTHOR Voorhees, E. A., (LANL)
14
C Subroutine SNBIR solves a general nonsymmetric banded NxN
15
C system of single precision real linear equations using
16
C SLATEC subroutines SNBFA and SNBSL. These are adaptations
17
C of the LINPACK subroutines SGBFA and SGBSL, which require
18
C a different format for storing the matrix elements.
19
C One pass of iterative refinement is used only to obtain an
20
C estimate of the accuracy. If A is an NxN real banded
21
C matrix and if X and B are real N-vectors, then SNBIR
26
C A band matrix is a matrix whose nonzero elements are all
27
C fairly near the main diagonal, specifically A(I,J) = 0
28
C if I-J is greater than ML or J-I is greater than
29
C MU . The integers ML and MU are called the lower and upper
30
C band widths and M = ML+MU+1 is the total band width.
31
C SNBIR uses less time and storage than the corresponding
32
C program for general matrices (SGEIR) if 2*ML+MU .LT. N .
34
C The matrix A is first factored into upper and lower tri-
35
C angular matrices U and L using partial pivoting. These
36
C factors and the pivoting information are used to find the
37
C solution vector X . Then the residual vector is found and used
38
C to calculate an estimate of the relative error, IND . IND esti-
39
C mates the accuracy of the solution only when the input matrix
40
C and the right hand side are represented exactly in the computer
41
C and does not take into account any errors in the input data.
43
C If the equation A*X=B is to be solved for more than one vector
44
C B, the factoring of A does not need to be performed again and
45
C the option to only solve (ITASK .GT. 1) will be faster for
46
C the succeeding solutions. In this case, the contents of A, LDA,
47
C N, work and IWORK must not have been altered by the user follow-
48
C ing factorization (ITASK=1). IND will not be changed by SNBIR
54
C If A is a band matrix, the following program segment
55
C will set up the input.
57
C ML = (band width below the diagonal)
58
C MU = (band width above the diagonal)
68
C This uses columns 1 Through ML+MU+1 of ABE .
70
C Example: If the original matrix is
79
C then N = 6, ML = 1, MU = 2, LDA .GE. 5 and ABE should contain
81
C * 11 12 13 , * = not used
89
C Argument Description ***
92
C on entry, contains the matrix in band storage as
93
C described above. MM must not be less than M =
94
C ML+MU+1 . The user is cautioned to dimension ABE
95
C with care since MM is not an argument and cannot
96
C be checked by SNBIR. The rows of the original
97
C matrix are stored in the rows of ABE and the
98
C diagonals of the original matrix are stored in
99
C columns 1 through ML+MU+1 of ABE . ABE is
100
C not altered by the program.
102
C the leading dimension of array ABE. LDA must be great-
103
C er than or equal to N. (terminal error message IND=-1)
105
C the order of the matrix A. N must be greater
106
C than or equal to 1 . (terminal error message IND=-2)
108
C the number of diagonals below the main diagonal.
109
C ML must not be less than zero nor greater than or
110
C equal to N . (terminal error message IND=-5)
112
C the number of diagonals above the main diagonal.
113
C MU must not be less than zero nor greater than or
114
C equal to N . (terminal error message IND=-6)
116
C on entry, the singly subscripted array(vector) of di-
117
C mension N which contains the right hand side B of a
118
C system of simultaneous linear equations A*X=B.
119
C on return, V contains the solution vector, X .
121
C If ITASK=1, the matrix A is factored and then the
122
C linear equation is solved.
123
C If ITASK .GT. 1, the equation is solved using the existing
124
C factored matrix A and IWORK.
125
C If ITASK .LT. 1, then terminal error message IND=-3 is
128
C GT. 0 IND is a rough estimate of the number of digits
129
C of accuracy in the solution, X . IND=75 means
130
C that the solution vector X is zero.
131
C LT. 0 See error message corresponding to IND below.
132
C WORK REAL(N*(NC+1))
133
C a singly subscripted array of dimension at least
134
C N*(NC+1) where NC = 2*ML+MU+1 .
136
C a singly subscripted array of dimension at least N.
138
C Error Messages Printed ***
140
C IND=-1 terminal N is greater than LDA.
141
C IND=-2 terminal N is less than 1.
142
C IND=-3 terminal ITASK is less than 1.
143
C IND=-4 terminal the matrix A is computationally singular.
144
C A solution has not been computed.
145
C IND=-5 terminal ML is less than zero or is greater than
147
C IND=-6 terminal MU is less than zero or is greater than
149
C IND=-10 warning the solution has no apparent significance.
150
C The solution may be inaccurate or the matrix
151
C A may be poorly scaled.
153
C Note- The above terminal(*fatal*) error messages are
154
C designed to be handled by XERMSG in which
155
C LEVEL=1 (recoverable) and IFLAG=2 . LEVEL=0
156
C for warning error messages from XERMSG. Unless
157
C the user provides otherwise, an error message
158
C will be printed followed by an abort.
160
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
161
C Stewart, LINPACK Users' Guide, SIAM, 1979.
162
C***ROUTINES CALLED R1MACH, SASUM, SCOPY, SDSDOT, SNBFA, SNBSL, XERMSG
163
C***REVISION HISTORY (YYMMDD)
164
C 800815 DATE WRITTEN
165
C 890531 Changed all specific intrinsics to generic. (WRB)
166
C 890831 Modified array declarations. (WRB)
167
C 890831 REVISION DATE from Version 3.2
168
C 891214 Prologue converted to Version 4.0 format. (BAB)
169
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
170
C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
171
C 920501 Reformatted the REFERENCES section. (WRB)
172
C***END PROLOGUE SNBIR
174
INTEGER LDA,N,ITASK,IND,IWORK(*),INFO,J,K,KK,L,M,ML,MU,NC
175
REAL ABE(LDA,*),V(*),WORK(N,*),XNORM,DNORM,SDSDOT,SASUM,R1MACH
176
CHARACTER*8 XERN1, XERN2
177
C***FIRST EXECUTABLE STATEMENT SNBIR
180
WRITE (XERN1, '(I8)') LDA
181
WRITE (XERN2, '(I8)') N
182
CALL XERMSG ('SLATEC', 'SNBIR', 'LDA = ' // XERN1 //
183
* ' IS LESS THAN N = ' // XERN2, -1, 1)
189
WRITE (XERN1, '(I8)') N
190
CALL XERMSG ('SLATEC', 'SNBIR', 'N = ' // XERN1 //
191
* ' IS LESS THAN 1', -2, 1)
197
WRITE (XERN1, '(I8)') ITASK
198
CALL XERMSG ('SLATEC', 'SNBIR', 'ITASK = ' // XERN1 //
199
* ' IS LESS THAN 1', -3, 1)
203
IF (ML.LT.0 .OR. ML.GE.N) THEN
205
WRITE (XERN1, '(I8)') ML
206
CALL XERMSG ('SLATEC', 'SNBIR',
207
* 'ML = ' // XERN1 // ' IS OUT OF RANGE', -5, 1)
211
IF (MU.LT.0 .OR. MU.GE.N) THEN
213
WRITE (XERN1, '(I8)') MU
214
CALL XERMSG ('SLATEC', 'SNBIR',
215
* 'MU = ' // XERN1 // ' IS OUT OF RANGE', -6, 1)
222
C MOVE MATRIX ABE TO WORK
226
CALL SCOPY(N,ABE(1,J),1,WORK(1,J),1)
229
C FACTOR MATRIX A INTO LU
231
CALL SNBFA(WORK,N,N,ML,MU,IWORK,INFO)
233
C CHECK FOR COMPUTATIONALLY SINGULAR MATRIX
237
CALL XERMSG ('SLATEC', 'SNBIR',
238
* 'SINGULAR MATRIX A - NO SOLUTION', -4, 1)
243
C SOLVE WHEN FACTORING COMPLETE
244
C MOVE VECTOR B TO WORK
246
CALL SCOPY(N,V(1),1,WORK(1,NC+1),1)
247
CALL SNBSL(WORK,N,N,ML,MU,IWORK,V,0)
251
XNORM = SASUM(N,V(1),1)
252
IF (XNORM.EQ.0.0) THEN
262
L = MIN(J-1,ML)+MIN(N-J,MU)+1
263
WORK(J,NC+1) = SDSDOT(L,-WORK(J,NC+1),ABE(J,K),LDA,V(KK),1)
268
CALL SNBSL(WORK,N,N,ML,MU,IWORK,WORK(1,NC+1),0)
272
DNORM = SASUM(N,WORK(1,NC+1),1)
274
C COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS)
275
C AND CHECK FOR IND GREATER THAN ZERO
277
IND = -LOG10(MAX(R1MACH(4),DNORM/XNORM))
280
CALL XERMSG ('SLATEC', 'SNBIR',
281
* 'SOLUTION MAY HAVE NO SIGNIFICANCE', -10, 0)
added
removed
deps/openlibm/slatec/snbir.f
