2
SUBROUTINE DCOV (FCN, IOPT, M, N, X, FVEC, R, LDR, INFO, WA1, WA2,
4
C***BEGIN PROLOGUE DCOV
5
C***PURPOSE Calculate the covariance matrix for a nonlinear data
6
C fitting problem. It is intended to be used after a
7
C successful return from either DNLS1 or DNLS1E.
10
C***TYPE DOUBLE PRECISION (SCOV-S, DCOV-D)
11
C***KEYWORDS COVARIANCE MATRIX, NONLINEAR DATA FITTING,
12
C NONLINEAR LEAST SQUARES
13
C***AUTHOR Hiebert, K. L., (SNLA)
18
C DCOV calculates the covariance matrix for a nonlinear data
19
C fitting problem. It is intended to be used after a
20
C successful return from either DNLS1 or DNLS1E. DCOV
21
C and DNLS1 (and DNLS1E) have compatible parameters. The
22
C required external subroutine, FCN, is the same
23
C for all three codes, DCOV, DNLS1, and DNLS1E.
25
C 2. Subroutine and Type Statements.
27
C SUBROUTINE DCOV(FCN,IOPT,M,N,X,FVEC,R,LDR,INFO,
29
C INTEGER IOPT,M,N,LDR,INFO
30
C DOUBLE PRECISION X(N),FVEC(M),R(LDR,N),WA1(N),WA2(N),WA3(N),WA4(M)
33
C 3. Parameters. All TYPE REAL parameters are DOUBLE PRECISION
35
C FCN is the name of the user-supplied subroutine which calculates
36
C the functions. If the user wants to supply the Jacobian
37
C (IOPT=2 or 3), then FCN must be written to calculate the
38
C Jacobian, as well as the functions. See the explanation
39
C of the IOPT argument below.
40
C If the user wants the iterates printed in DNLS1 or DNLS1E,
41
C then FCN must do the printing. See the explanation of NPRINT
42
C in DNLS1 or DNLS1E. FCN must be declared in an EXTERNAL
43
C statement in the calling program and should be written as
46
C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
47
C INTEGER IFLAG,LDFJAC,M,N
48
C DOUBLE PRECISION X(N),FVEC(M)
50
C FJAC and LDFJAC may be ignored , if IOPT=1.
51
C DOUBLE PRECISION FJAC(LDFJAC,N) , if IOPT=2.
52
C DOUBLE PRECISION FJAC(N) , if IOPT=3.
54
C If IFLAG=0, the values in X and FVEC are available
55
C for printing in DNLS1 or DNLS1E.
56
C IFLAG will never be zero when FCN is called by DCOV.
57
C The values of X and FVEC must not be changed.
60
C If IFLAG=1, calculate the functions at X and return
61
C this vector in FVEC.
64
C If IFLAG=2, calculate the full Jacobian at X and return
65
C this matrix in FJAC. Note that IFLAG will never be 2 unless
66
C IOPT=2. FVEC contains the function values at X and must
67
C not be altered. FJAC(I,J) must be set to the derivative
68
C of FVEC(I) with respect to X(J).
71
C If IFLAG=3, calculate the LDFJAC-th row of the Jacobian
72
C and return this vector in FJAC. Note that IFLAG will
73
C never be 3 unless IOPT=3. FJAC(J) must be set to
74
C the derivative of FVEC(LDFJAC) with respect to X(J).
80
C The value of IFLAG should not be changed by FCN unless the
81
C user wants to terminate execution of DCOV. In this case, set
82
C IFLAG to a negative integer.
85
C IOPT is an input variable which specifies how the Jacobian will
86
C be calculated. If IOPT=2 or 3, then the user must supply the
87
C Jacobian, as well as the function values, through the
88
C subroutine FCN. If IOPT=2, the user supplies the full
89
C Jacobian with one call to FCN. If IOPT=3, the user supplies
90
C one row of the Jacobian with each call. (In this manner,
91
C storage can be saved because the full Jacobian is not stored.)
92
C If IOPT=1, the code will approximate the Jacobian by forward
95
C M is a positive integer input variable set to the number of
98
C N is a positive integer input variable set to the number of
99
C variables. N must not exceed M.
101
C X is an array of length N. On input X must contain the value
102
C at which the covariance matrix is to be evaluated. This is
103
C usually the value for X returned from a successful run of
104
C DNLS1 (or DNLS1E). The value of X will not be changed.
106
C FVEC is an output array of length M which contains the functions
109
C R is an output array. For IOPT=1 and 2, R is an M by N array.
110
C For IOPT=3, R is an N by N array. On output, if INFO=1,
111
C the upper N by N submatrix of R contains the covariance
112
C matrix evaluated at X.
114
C LDR is a positive integer input variable which specifies
115
C the leading dimension of the array R. For IOPT=1 and 2,
116
C LDR must not be less than M. For IOPT=3, LDR must not
119
C INFO is an integer output variable. If the user has terminated
120
C execution, INFO is set to the (negative) value of IFLAG. See
121
C description of FCN. Otherwise, INFO is set as follows.
123
C INFO = 0 Improper input parameters (M.LE.0 or N.LE.0).
125
C INFO = 1 Successful return. The covariance matrix has been
126
C calculated and stored in the upper N by N
129
C INFO = 2 The Jacobian matrix is singular for the input value
130
C of X. The covariance matrix cannot be calculated.
131
C The upper N by N submatrix of R contains the QR
132
C factorization of the Jacobian (probably not of
133
C interest to the user).
135
C WA1,WA2 are work arrays of length N.
138
C WA4 is a work array of length M.
140
C***REFERENCES (NONE)
141
C***ROUTINES CALLED DENORM, DFDJC3, DQRFAC, DWUPDT, XERMSG
142
C***REVISION HISTORY (YYMMDD)
143
C 810522 DATE WRITTEN
144
C 890831 Modified array declarations. (WRB)
145
C 891006 Cosmetic changes to prologue. (WRB)
146
C 891006 REVISION DATE from Version 3.2
147
C 891214 Prologue converted to Version 4.0 format. (BAB)
148
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
149
C 900510 Fixed an error message. (RWC)
150
C***END PROLOGUE DCOV
152
C REVISED 850601-1100
153
C REVISED YYMMDD HHMM
155
INTEGER I,IDUM,IFLAG,INFO,IOPT,J,K,KP1,LDR,M,N,NM1,NROW
156
DOUBLE PRECISION X(*),R(LDR,*),FVEC(*),WA1(*),WA2(*),WA3(*),
159
DOUBLE PRECISION ONE,SIGMA,TEMP,ZERO,DENORM
162
DATA ZERO/0.D0/,ONE/1.D0/
163
C***FIRST EXECUTABLE STATEMENT DCOV
166
IF (M.LE.0 .OR. N.LE.0) GO TO 300
168
C CALCULATE SIGMA = (SUM OF THE SQUARED RESIDUALS) / (M-N)
170
CALL FCN(IFLAG,M,N,X,FVEC,R,LDR)
171
IF (IFLAG.LT.0) GO TO 300
174
IF (M.NE.N) SIGMA=TEMP*TEMP/(M-N)
176
C CALCULATE THE JACOBIAN
177
IF (IOPT.EQ.3) GO TO 200
179
C STORE THE FULL JACOBIAN USING M*N STORAGE
180
IF (IOPT.EQ.1) GO TO 100
182
C USER SUPPLIES THE JACOBIAN
184
CALL FCN(IFLAG,M,N,X,FVEC,R,LDR)
187
C CODE APPROXIMATES THE JACOBIAN
188
100 CALL DFDJC3(FCN,M,N,X,FVEC,R,LDR,IFLAG,ZERO,WA4)
189
110 IF (IFLAG.LT.0) GO TO 300
191
C COMPUTE THE QR DECOMPOSITION
192
CALL DQRFAC(M,N,R,LDR,.FALSE.,IDUM,1,WA1,WA1,WA1)
197
C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX CALCULATED ONE
198
C ROW AT A TIME AND STORED IN THE UPPER TRIANGLE OF R.
199
C ( (Q TRANSPOSE)*FVEC IS ALSO CALCULATED BUT NOT USED.)
210
CALL FCN(IFLAG,M,N,X,FVEC,WA1,NROW)
211
IF (IFLAG.LT.0) GO TO 300
213
CALL DWUPDT(N,R,LDR,WA1,WA2,TEMP,WA3,WA4)
216
C CHECK IF R IS SINGULAR.
219
IF (R(I,I).EQ.ZERO) SING=.TRUE.
223
C R IS UPPER TRIANGULAR. CALCULATE (R TRANSPOSE) INVERSE AND STORE
224
C IN THE UPPER TRIANGLE OF R.
225
IF (N.EQ.1) GO TO 275
229
C INITIALIZE THE RIGHT-HAND SIDE (WA1(*)) AS THE K-TH COLUMN OF THE
240
C SUBTRACT R(K,I-1)*R(I-1,*) FROM THE RIGHT-HAND SIDE, WA1(*).
242
WA1(J)=WA1(J)-R(K,I-1)*R(I-1,J)
247
275 R(N,N)=ONE/R(N,N)
249
C CALCULATE R-INVERSE * (R TRANSPOSE) INVERSE AND STORE IN THE UPPER
255
TEMP=TEMP+R(I,K)*R(J,K)
262
IF (M.LE.0 .OR. N.LE.0) INFO=0
263
IF (IFLAG.LT.0) INFO=IFLAG
265
IF (INFO .LT. 0) CALL XERMSG ('SLATEC', 'DCOV',
266
+ 'EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.', 1, 1)
267
IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'DCOV',
268
+ 'INVALID INPUT PARAMETER.', 2, 1)
269
IF (INFO .EQ. 2) CALL XERMSG ('SLATEC', 'DCOV',
270
+ 'SINGULAR JACOBIAN MATRIX, COVARIANCE MATRIX CANNOT BE ' //
271
+ 'CALCULATED.', 1, 1)
added
removed
deps/openlibm/slatec/dcov.f
