| [fcn,iopt,m,n,x,fvec,r,ldr,info,wa1,wa2,wa3,wa4]=dcov(fcn,iopt,m,n,x,fvec,r,ldr,info,wa1,wa2,wa3,wa4); |
function [fcn,iopt,m,n,x,fvec,r,ldr,info,wa1,wa2,wa3,wa4]=dcov(fcn,iopt,m,n,x,fvec,r,ldr,info,wa1,wa2,wa3,wa4);
%***BEGIN PROLOGUE DCOV
%***PURPOSE Calculate the covariance matrix for a nonlinear data
% fitting problem. It is intended to be used after a
% successful return from either DNLS1 or DNLS1E.
%***LIBRARY SLATEC
%***CATEGORY K1B1
%***TYPE doubleprecision (SCOV-S, DCOV-D)
%***KEYWORDS COVARIANCE MATRIX, NONLINEAR DATA FITTING,
% NONLINEAR LEAST SQUARES
%***AUTHOR Hiebert, K. L., (SNLA)
%***DESCRIPTION
%
% 1. Purpose.
%
% DCOV calculates the covariance matrix for a nonlinear data
% fitting problem. It is intended to be used after a
% successful return from either DNLS1 or DNLS1E. DCOV
% and DNLS1 (and DNLS1E) have compatible parameters. The
% required external subroutine, FCN, is the same
% for all three codes, DCOV, DNLS1, and DNLS1E.
%
% 2. subroutine and Type Statements.
%
% subroutine DCOV(FCN,IOPT,M,N,X,FVEC,R,LDR,INFO,
% WA1,WA2,WA3,WA4)
% INTEGER IOPT,M,N,LDR,INFO
% doubleprecision X(N),FVEC(M),R(LDR,N),WA1(N),WA2(N),WA3(N),WA4(M)
% EXTERNAL FCN
%
% 3. Parameters. All TYPE REAL parameters are doubleprecision
%
% FCN is the name of the user-supplied subroutine which calculates
% the functions. If the user wants to supply the Jacobian
% (IOPT=2 or 3), then FCN must be written to calculate the
% Jacobian, as well as the functions. See the explanation
% of the IOPT argument below.
% If the user wants the iterates printed in DNLS1 or DNLS1E,
% then FCN must do the printing. See the explanation of NPRINT
% in DNLS1 or DNLS1E. FCN must be declared in an EXTERNAL
% statement in the calling program and should be written as
% follows.
%
% subroutine FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
% INTEGER IFLAG,LDFJAC,M,N
% doubleprecision X(N),FVEC(M)
% ----------
% FJAC and LDFJAC may be ignored , if IOPT=1.
% doubleprecision FJAC(LDFJAC,N) , if IOPT=2.
% doubleprecision FJAC(N) , if IOPT=3.
% ----------
% If IFLAG=0, the values in X and FVEC are available
% for printing in DNLS1 or DNLS1E.
% IFLAG will never be zero when FCN is called by DCOV.
% The values of X and FVEC must not be changed.
% RETURN
% ----------
% If IFLAG=1, calculate the functions at X and return
% this vector in FVEC.
% RETURN
% ----------
% If IFLAG=2, calculate the full Jacobian at X and return
% this matrix in FJAC. Note that IFLAG will never be 2 unless
% IOPT=2. FVEC contains the function values at X and must
% not be altered. FJAC(I,J) must be set to the derivative
% of FVEC(I) with respect to X(J).
% RETURN
% ----------
% If IFLAG=3, calculate the LDFJAC-th row of the Jacobian
% and return this vector in FJAC. Note that IFLAG will
% never be 3 unless IOPT=3. FJAC(J) must be set to
% the derivative of FVEC(LDFJAC) with respect to X(J).
% RETURN
% ----------
% end
%
%
% The value of IFLAG should not be changed by FCN unless the
% user wants to terminate execution of DCOV. In this case, set
% IFLAG to a negative integer.
%
%
% IOPT is an input variable which specifies how the Jacobian will
% be calculated. If IOPT=2 or 3, then the user must supply the
% Jacobian, as well as the function values, through the
% subroutine FCN. If IOPT=2, the user supplies the full
% Jacobian with one call to FCN. If IOPT=3, the user supplies
% one row of the Jacobian with each call. (In this manner,
% storage can be saved because the full Jacobian is not stored.)
% If IOPT=1, the code will approximate the Jacobian by forward
% differencing.
%
% M is a positive integer input variable set to the number of
% functions.
%
% N is a positive integer input variable set to the number of
% variables. N must not exceed M.
%
% X is an array of length N. On input X must contain the value
% at which the covariance matrix is to be evaluated. This is
% usually the value for X returned from a successful run of
% DNLS1 (or DNLS1E). The value of X will not be changed.
%
% FVEC is an output array of length M which contains the functions
% evaluated at X.
%
% R is an output array. For IOPT=1 and 2, R is an M by N array.
% For IOPT=3, R is an N by N array. On output, if INFO=1,
% the upper N by N submatrix of R contains the covariance
% matrix evaluated at X.
%
% LDR is a positive integer input variable which specifies
% the leading dimension of the array R. For IOPT=1 and 2,
% LDR must not be less than M. For IOPT=3, LDR must not
% be less than N.
%
% INFO is an integer output variable. If the user has terminated
% execution, INFO is set to the (negative) value of IFLAG. See
% description of FCN. Otherwise, INFO is set as follows.
%
% INFO = 0 Improper input parameters (M.LE.0 or N.LE.0).
%
% INFO = 1 Successful return. The covariance matrix has been
% calculated and stored in the upper N by N
% submatrix of R.
%
% INFO = 2 The Jacobian matrix is singular for the input value
% of X. The covariance matrix cannot be calculated.
% The upper N by N submatrix of R contains the QR
% factorization of the Jacobian (probably not of
% interest to the user).
%
% WA1,WA2 are work arrays of length N.
% and WA3
%
% WA4 is a work array of length M.
%
%***REFERENCES (NONE)
%***ROUTINES CALLED DENORM, DFDJC3, DQRFAC, DWUPDT, XERMSG
%***REVISION HISTORY (YYMMDD)
% 810522 DATE WRITTEN
% 890831 Modified array declarations. (WRB)
% 891006 Cosmetic changes to prologue. (WRB)
% 891006 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 900510 Fixed an error message. (RWC)
%***end PROLOGUE DCOV
%
% REVISED 850601-1100
% REVISED YYMMDD HHMM
%
persistent firstCall gt i idum iflag j k kp1 nm1 nrow one sigma sing temp zero ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(idum), idum=0; end;
if isempty(iflag), iflag=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kp1), kp1=0; end;
if isempty(nm1), nm1=0; end;
if isempty(nrow), nrow=0; end;
if isempty(gt), gt=0; end;
x_shape=size(x);x=reshape(x,1,[]);
r_shape=size(r);r=reshape([r(:).',zeros(1,ceil(numel(r)./prod([ldr])).*prod([ldr])-numel(r))],ldr,[]);
fvec_shape=size(fvec);fvec=reshape(fvec,1,[]);
wa1_shape=size(wa1);wa1=reshape(wa1,1,[]);
wa2_shape=size(wa2);wa2=reshape(wa2,1,[]);
wa3_shape=size(wa3);wa3=reshape(wa3,1,[]);
wa4_shape=size(wa4);wa4=reshape(wa4,1,[]);
if isempty(one), one=0; end;
if isempty(sigma), sigma=0; end;
if isempty(temp), temp=0; end;
if isempty(zero), zero=0; end;
if isempty(sing), sing=false; end;
if firstCall, zero=[0.0d0]; end;
if firstCall, one=[1.0d0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DCOV
gt=0;
sing = false;
iflag = 0;
while (1);
if( m>0 && n>0 )
%
% CALCULATE SIGMA = (SUM OF THE SQUARED RESIDUALS) / (M-N)
iflag = 1;
[iflag,m,n,x,fvec,r,ldr]=fcn(iflag,m,n,x,fvec,r,ldr);
if( iflag>=0 )
[temp ,m,fvec]=denorm(m,fvec);
sigma = one;
if( m~=n )
sigma = temp.*temp./(m-n);
end;
%
% CALCULATE THE JACOBIAN
if( iopt==3 )
%
% COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX CALCULATED ONE
% ROW AT A TIME AND STORED IN THE UPPER TRIANGLE OF R.
% ( (Q TRANSPOSE)*FVEC IS ALSO CALCULATED BUT NOT USED.)
for j = 1 : n;
wa2(j) = zero;
for i = 1 : n;
r(i,j) = zero;
end; i = fix(n+1);
end; j = fix(n+1);
iflag = 3;
for i = 1 : m;
nrow = fix(i);
[iflag,m,n,x,fvec,wa1,nrow]=fcn(iflag,m,n,x,fvec,wa1,nrow);
if( iflag<0 )
gt=1;
break;
end;
temp = fvec(i);
[n,r,ldr,wa1,wa2,temp,wa3,wa4]=dwupdt(n,r,ldr,wa1,wa2,temp,wa3,wa4);
end;
if(gt==1)
break;
end;
else;
%
% STORE THE FULL JACOBIAN USING M*N STORAGE
if( iopt==1 )
%
% CODE APPROXIMATES THE JACOBIAN
[fcn,m,n,x,fvec,r,ldr,iflag,zero,wa4]=dfdjc3(fcn,m,n,x,fvec,r,ldr,iflag,zero,wa4);
else;
%
% USER SUPPLIES THE JACOBIAN
iflag = 2;
[iflag,m,n,x,fvec,r,ldr]=fcn(iflag,m,n,x,fvec,r,ldr);
end;
if( iflag<0 )
gt=1;
break;
end;
%
% COMPUTE THE QR DECOMPOSITION
wa1_orig=wa1; wa1_orig=wa1; [m,n,r,ldr,dumvar5,idum,dumvar7,wa1,dumvar9,dumvar10]=dqrfac(m,n,r,ldr,false,idum,1,wa1,wa1,wa1); wa1(dumvar10~=wa1_orig)=dumvar10(dumvar10~=wa1_orig); wa1(dumvar9~=wa1_orig)=dumvar9(dumvar9~=wa1_orig);
for i = 1 : n;
r(i,i) = wa1(i);
end; i = fix(n+1);
end;
%
% CHECK IF R IS SINGULAR.
for i = 1 : n;
if( r(i,i)==zero )
sing = true;
end;
end; i = fix(n+1);
if( ~(sing) )
%
% R IS UPPER TRIANGULAR. CALCULATE (R TRANSPOSE) INVERSE AND STORE
% IN THE UPPER TRIANGLE OF R.
if( n~=1 )
nm1 = fix(n - 1);
for k = 1 : nm1;
%
% INITIALIZE THE RIGHT-HAND SIDE (WA1(*)) AS THE K-TH COLUMN OF THE
% IDENTITY MATRIX.
for i = 1 : n;
wa1(i) = zero;
end; i = fix(n+1);
wa1(k) = one;
%
r(k,k) = wa1(k)./r(k,k);
kp1 = fix(k + 1);
for i = kp1 : n;
%
% SUBTRACT R(K,I-1)*R(I-1,*) FROM THE RIGHT-HAND SIDE, WA1(*).
for j = i : n;
wa1(j) = wa1(j) - r(k,i-1).*r(i-1,j);
end; j = fix(n+1);
r(k,i) = wa1(i)./r(i,i);
end; i = fix(n+1);
end; k = fix(nm1+1);
end;
r(n,n) = one./r(n,n);
%
% CALCULATE R-INVERSE * (R TRANSPOSE) INVERSE AND STORE IN THE UPPER
% TRIANGLE OF R.
for i = 1 : n;
for j = i : n;
temp = zero;
for k = j : n;
temp = temp + r(i,k).*r(j,k);
end; k = fix(n+1);
r(i,j) = temp.*sigma;
end; j = fix(n+1);
end; i = fix(n+1);
info = 1;
end;
end;
end;
break;
end;
%
if( m<=0 || n<=0 )
info = 0;
end;
if( iflag<0 )
info = fix(iflag);
end;
if( sing )
info = 2;
end;
if( info<0 )
xermsg('SLATEC','DCOV','EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.',1,1);
end;
if( info==0 )
xermsg('SLATEC','DCOV','INVALID INPUT PARAMETER.',2,1);
end;
if( info==2 )
xermsg('SLATEC','DCOV',['SINGULAR JACOBIAN MATRIX, COVARIANCE MATRIX CANNOT BE ','CALCULATED.'],1,1);
end;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
fvec_shape=zeros(fvec_shape);fvec_shape(:)=fvec(1:numel(fvec_shape));fvec=fvec_shape;
wa1_shape=zeros(wa1_shape);wa1_shape(:)=wa1(1:numel(wa1_shape));wa1=wa1_shape;
wa2_shape=zeros(wa2_shape);wa2_shape(:)=wa2(1:numel(wa2_shape));wa2=wa2_shape;
wa3_shape=zeros(wa3_shape);wa3_shape(:)=wa3(1:numel(wa3_shape));wa3=wa3_shape;
wa4_shape=zeros(wa4_shape);wa4_shape(:)=wa4(1:numel(wa4_shape));wa4=wa4_shape;
end %subroutine dcov
%DECK DCPPLT
|
|
