2
SUBROUTINE DBNDAC (G, MDG, NB, IP, IR, MT, JT)
3
C***BEGIN PROLOGUE DBNDAC
4
C***PURPOSE Compute the LU factorization of a banded matrices using
5
C sequential accumulation of rows of the data matrix.
6
C Exactly one right-hand side vector is permitted.
9
C***TYPE DOUBLE PRECISION (BNDACC-S, DBNDAC-D)
10
C***KEYWORDS BANDED MATRIX, CURVE FITTING, LEAST SQUARES
11
C***AUTHOR Lawson, C. L., (JPL)
12
C Hanson, R. J., (SNLA)
15
C These subroutines solve the least squares problem Ax = b for
16
C banded matrices A using sequential accumulation of rows of the
17
C data matrix. Exactly one right-hand side vector is permitted.
19
C These subroutines are intended for the type of least squares
20
C systems that arise in applications such as curve or surface
21
C fitting of data. The least squares equations are accumulated and
22
C processed using only part of the data. This requires a certain
23
C user interaction during the solution of Ax = b.
25
C Specifically, suppose the data matrix (A B) is row partitioned
26
C into Q submatrices. Let (E F) be the T-th one of these
27
C submatrices where E = (0 C 0). Here the dimension of E is MT by N
28
C and the dimension of C is MT by NB. The value of NB is the
29
C bandwidth of A. The dimensions of the leading block of zeros in E
32
C The user of the subroutine DBNDAC provides MT,JT,C and F for
33
C T=1,...,Q. Not all of this data must be supplied at once.
35
C Following the processing of the various blocks (E F), the matrix
36
C (A B) has been transformed to the form (R D) where R is upper
37
C triangular and banded with bandwidth NB. The least squares
38
C system Rx = d is then easily solved using back substitution by
39
C executing the statement CALL DBNDSL(1,...). The sequence of
40
C values for JT must be nondecreasing. This may require some
41
C preliminary interchanges of rows and columns of the matrix A.
43
C The primary reason for these subroutines is that the total
44
C processing can take place in a working array of dimension MU by
45
C NB+1. An acceptable value for MU is
47
C MU = MAX(MT + N + 1),
49
C where N is the number of unknowns.
51
C Here the maximum is taken over all values of MT for T=1,...,Q.
52
C Notice that MT can be taken to be a small as one, showing that
53
C MU can be as small as N+2. The subprogram DBNDAC processes the
54
C rows more efficiently if MU is large enough so that each new
55
C block (C F) has a distinct value of JT.
57
C The four principle parts of these algorithms are obtained by the
58
C following call statements
60
C CALL DBNDAC(...) Introduce new blocks of data.
62
C CALL DBNDSL(1,...)Compute solution vector and length of
65
C CALL DBNDSL(2,...)Given any row vector H solve YR = H for the
68
C CALL DBNDSL(3,...)Given any column vector W solve RZ = W for
69
C the column vector Z.
71
C The dots in the above call statements indicate additional
72
C arguments that will be specified in the following paragraphs.
74
C The user must dimension the array appearing in the call list..
77
C Description of calling sequence for DBNDAC..
79
C The entire set of parameters for DBNDAC are
81
C Input.. All Type REAL variables are DOUBLE PRECISION
83
C G(*,*) The working array into which the user will
84
C place the MT by NB+1 block (C F) in rows IR
85
C through IR+MT-1, columns 1 through NB+1.
86
C See descriptions of IR and MT below.
88
C MDG The number of rows in the working array
89
C G(*,*). The value of MDG should be .GE. MU.
90
C The value of MU is defined in the abstract
91
C of these subprograms.
93
C NB The bandwidth of the data matrix A.
95
C IP Set by the user to the value 1 before the
96
C first call to DBNDAC. Its subsequent value
97
C is controlled by DBNDAC to set up for the
98
C next call to DBNDAC.
100
C IR Index of the row of G(*,*) where the user is
101
C to place the new block of data (C F). Set by
102
C the user to the value 1 before the first call
103
C to DBNDAC. Its subsequent value is controlled
104
C by DBNDAC. A value of IR .GT. MDG is considered
107
C MT,JT Set by the user to indicate respectively the
108
C number of new rows of data in the block and
109
C the index of the first nonzero column in that
110
C set of rows (E F) = (0 C 0 F) being processed.
112
C Output.. All Type REAL variables are DOUBLE PRECISION
114
C G(*,*) The working array which will contain the
115
C processed rows of that part of the data
116
C matrix which has been passed to DBNDAC.
118
C IP,IR The values of these arguments are advanced by
119
C DBNDAC to be ready for storing and processing
120
C a new block of data in G(*,*).
122
C Description of calling sequence for DBNDSL..
124
C The user must dimension the arrays appearing in the call list..
128
C The entire set of parameters for DBNDSL are
130
C Input.. All Type REAL variables are DOUBLE PRECISION
132
C MODE Set by the user to one of the values 1, 2, or
133
C 3. These values respectively indicate that
134
C the solution of AX = B, YR = H or RZ = W is
137
C G(*,*),MDG, These arguments all have the same meaning and
138
C NB,IP,IR contents as following the last call to DBNDAC.
140
C X(*) With mode=2 or 3 this array contains,
141
C respectively, the right-side vectors H or W of
142
C the systems YR = H or RZ = W.
144
C N The number of variables in the solution
145
C vector. If any of the N diagonal terms are
146
C zero the subroutine DBNDSL prints an
147
C appropriate message. This condition is
148
C considered an error.
150
C Output.. All Type REAL variables are DOUBLE PRECISION
152
C X(*) This array contains the solution vectors X,
153
C Y or Z of the systems AX = B, YR = H or
154
C RZ = W depending on the value of MODE=1,
157
C RNORM If MODE=1 RNORM is the Euclidean length of the
158
C residual vector AX-B. When MODE=2 or 3 RNORM
163
C To obtain the upper triangular matrix and transformed right-hand
164
C side vector D so that the super diagonals of R form the columns
165
C of G(*,*), execute the following Fortran statements.
177
C CALL DBNDAC(G,MDG,NB,IP,IR,MT,JT)
179
C***REFERENCES C. L. Lawson and R. J. Hanson, Solving Least Squares
180
C Problems, Prentice-Hall, Inc., 1974, Chapter 27.
181
C***ROUTINES CALLED DH12, XERMSG
182
C***REVISION HISTORY (YYMMDD)
183
C 790101 DATE WRITTEN
184
C 890531 Changed all specific intrinsics to generic. (WRB)
185
C 891006 Cosmetic changes to prologue. (WRB)
186
C 891006 REVISION DATE from Version 3.2
187
C 891214 Prologue converted to Version 4.0 format. (BAB)
188
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
189
C 920501 Reformatted the REFERENCES section. (WRB)
190
C***END PROLOGUE DBNDAC
191
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
193
C***FIRST EXECUTABLE STATEMENT DBNDAC
196
C ALG. STEPS 1-4 ARE PERFORMED EXTERNAL TO THIS SUBROUTINE.
199
IF (MT.LE.0.OR.NB.LE.0) RETURN
201
IF(.NOT.MDG.LT.IR) GO TO 5
204
CALL XERMSG ('SLATEC', 'DBNDAC', 'MDG.LT.IR, PROBABLE ERROR.',
210
IF (JT.EQ.IP) GO TO 70
212
IF (JT.LE.IR) GO TO 30
230
30 MU=MIN(NB-1,IR-IP-1)
231
IF (MU.EQ.0) GO TO 60
255
CALL DH12 (1,I,MAX(I+1,IR-IP+1),MH,G(IP,I),1,RHO,
256
1 G(IP,I+1),1,MDG,NBP1-I)
261
IF (KH.LT.NBP1) GO TO 100
added
removed
deps/openlibm/slatec/dbndac.f
