| [fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,ml,mu,epsfcn,diag,mode,factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,wa3,wa4]=dnsq(fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,ml,mu,epsfcn,diag,mode,factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,wa3,wa4); |
function [fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,ml,mu,epsfcn,diag,mode,factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,wa3,wa4]=dnsq(fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,ml,mu,epsfcn,diag,mode,factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,wa3,wa4);
%***BEGIN PROLOGUE DNSQ
%***PURPOSE Find a zero of a system of a N nonlinear functions in N
% variables by a modification of the Powell hybrid method.
%***LIBRARY SLATEC
%***CATEGORY F2A
%***TYPE doubleprecision (SNSQ-S, DNSQ-D)
%***KEYWORDS NONLINEAR SQUARE SYSTEM, POWELL HYBRID METHOD, ZEROS
%***AUTHOR Hiebert, K. L. (SNLA)
%***DESCRIPTION
%
% 1. Purpose.
%
% The purpose of DNSQ is to find a zero of a system of N nonlinear
% functions in N variables by a modification of the Powell
% hybrid method. The user must provide a subroutine which
% calculates the functions. The user has the option of either to
% provide a subroutine which calculates the Jacobian or to let the
% code calculate it by a forward-difference approximation.
% This code is the combination of the MINPACK codes (Argonne)
% HYBRD and HYBRDJ.
%
% 2. subroutine and Type Statements.
%
% subroutine DNSQ(FCN,JAC,IOPT,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,
% * ML,MU,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,
% * NJEV,R,LR,QTF,WA1,WA2,WA3,WA4)
% INTEGER IOPT,N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,NJEV,LR
% doubleprecision XTOL,EPSFCN,FACTOR
% doubleprecision
% X(N),FVEC(N),DIAG(N),FJAC(LDFJAC,N),R(LR),QTF(N),
% * WA1(N),WA2(N),WA3(N),WA4(N)
% EXTERNAL FCN,JAC
%
% 3. Parameters.
%
% Parameters designated as input parameters must be specified on
% entry to DNSQ and are not changed on exit, while parameters
% designated as output parameters need not be specified on entry
% and are set to appropriate values on exit from DNSQ.
%
% FCN is the name of the user-supplied subroutine which calculates
% the functions. FCN must be declared in an EXTERNAL statement
% in the user calling program, and should be written as follows.
%
% subroutine FCN(N,X,FVEC,IFLAG)
% INTEGER N,IFLAG
% doubleprecision X(N),FVEC(N)
% ----------
% CALCULATE THE FUNCTIONS AT X AND
% RETURN THIS VECTOR IN FVEC.
% ----------
% RETURN
% end
%
% The value of IFLAG should not be changed by FCN unless the
% user wants to terminate execution of DNSQ. In this case set
% IFLAG to a negative integer.
%
% JAC is the name of the user-supplied subroutine which calculates
% the Jacobian. If IOPT=1, then JAC must be declared in an
% EXTERNAL statement in the user calling program, and should be
% written as follows.
%
% subroutine JAC(N,X,FVEC,FJAC,LDFJAC,IFLAG)
% INTEGER N,LDFJAC,IFLAG
% doubleprecision X(N),FVEC(N),FJAC(LDFJAC,N)
% ----------
% Calculate the Jacobian at X and return this
% matrix in FJAC. FVEC contains the function
% values at X and should not be altered.
% ----------
% RETURN
% end
%
% The value of IFLAG should not be changed by JAC unless the
% user wants to terminate execution of DNSQ. In this case set
% IFLAG to a negative integer.
%
% If IOPT=2, JAC can be ignored (treat it as a dummy argument).
%
% IOPT is an input variable which specifies how the Jacobian will
% be calculated. If IOPT=1, then the user must supply the
% Jacobian through the subroutine JAC. If IOPT=2, then the
% code will approximate the Jacobian by forward-differencing.
%
% N is a positive integer input variable set to the number of
% functions and variables.
%
% X is an array of length N. On input X must contain an initial
% estimate of the solution vector. On output X contains the
% final estimate of the solution vector.
%
% FVEC is an output array of length N which contains the functions
% evaluated at the output X.
%
% FJAC is an output N by N array which contains the orthogonal
% matrix Q produced by the QR factorization of the final
% approximate Jacobian.
%
% LDFJAC is a positive integer input variable not less than N
% which specifies the leading dimension of the array FJAC.
%
% XTOL is a nonnegative input variable. Termination occurs when
% the relative error between two consecutive iterates is at most
% XTOL. Therefore, XTOL measures the relative error desired in
% the approximate solution. Section 4 contains more details
% about XTOL.
%
% MAXFEV is a positive integer input variable. Termination occurs
% when the number of calls to FCN is at least MAXFEV by the end
% of an iteration.
%
% ML is a nonnegative integer input variable which specifies the
% number of subdiagonals within the band of the Jacobian matrix.
% If the Jacobian is not banded or IOPT=1, set ML to at
% least N - 1.
%
% MU is a nonnegative integer input variable which specifies the
% number of superdiagonals within the band of the Jacobian
% matrix. If the Jacobian is not banded or IOPT=1, set MU to at
% least N - 1.
%
% EPSFCN is an input variable used in determining a suitable step
% for the forward-difference approximation. This approximation
% assumes that the relative errors in the functions are of the
% order of EPSFCN. If EPSFCN is less than the machine
% precision, it is assumed that the relative errors in the
% functions are of the order of the machine precision. If
% IOPT=1, then EPSFCN can be ignored (treat it as a dummy
% argument).
%
% DIAG is an array of length N. If MODE = 1 (see below), DIAG is
% internally set. If MODE = 2, DIAG must contain positive
% entries that serve as implicit (multiplicative) scale factors
% for the variables.
%
% MODE is an integer input variable. If MODE = 1, the variables
% will be scaled internally. If MODE = 2, the scaling is
% specified by the input DIAG. Other values of MODE are
% equivalent to MODE = 1.
%
% FACTOR is a positive input variable used in determining the
% initial step bound. This bound is set to the product of
% FACTOR and the Euclidean norm of DIAG*X if nonzero, or else to
% FACTOR itself. In most cases FACTOR should lie in the
% interval (.1,100.). 100. is a generally recommended value.
%
% NPRINT is an integer input variable that enables controlled
% printing of iterates if it is positive. In this case, FCN is
% called with IFLAG = 0 at the beginning of the first iteration
% and every NPRINT iterations thereafter and immediately prior
% to return, with X and FVEC available for printing. appropriate
% print statements must be added to FCN(see example). If NPRINT
% is not positive, no special calls of FCN with IFLAG = 0 are
% made.
%
% 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 and JAC. Otherwise, INFO is set as follows.
%
% INFO = 0 Improper input parameters.
%
% INFO = 1 Relative error between two consecutive iterates is
% at most XTOL.
%
% INFO = 2 Number of calls to FCN has reached or exceeded
% MAXFEV.
%
% INFO = 3 XTOL is too small. No further improvement in the
% approximate solution X is possible.
%
% INFO = 4 Iteration is not making good progress, as measured
% by the improvement from the last five Jacobian
% evaluations.
%
% INFO = 5 Iteration is not making good progress, as measured
% by the improvement from the last ten iterations.
%
% Sections 4 and 5 contain more details about INFO.
%
% NFEV is an integer output variable set to the number of calls to
% FCN.
%
% NJEV is an integer output variable set to the number of calls to
% JAC. (If IOPT=2, then NJEV is set to zero.)
%
% R is an output array of length LR which contains the upper
% triangular matrix produced by the QR factorization of the
% final approximate Jacobian, stored rowwise.
%
% LR is a positive integer input variable not less than
% (N*(N+1))/2.
%
% QTF is an output array of length N which contains the vector
% (Q transpose)*FVEC.
%
% WA1, WA2, WA3, and WA4 are work arrays of length N.
%
%
% 4. Successful completion.
%
% The accuracy of DNSQ is controlled by the convergence parameter
% XTOL. This parameter is used in a test which makes a comparison
% between the approximation X and a solution XSOL. DNSQ
% terminates when the test is satisfied. If the convergence
% parameter is less than the machine precision (as defined by the
% function D1MACH(4)), then DNSQ only attempts to satisfy the test
% defined by the machine precision. Further progress is not
% usually possible.
%
% The test assumes that the functions are reasonably well behaved,
% and, if the Jacobian is supplied by the user, that the functions
% and the Jacobian are coded consistently. If these conditions
% are not satisfied, then DNSQ may incorrectly indicate
% convergence. The coding of the Jacobian can be checked by the
% subroutine DCKDER. If the Jacobian is coded correctly or IOPT=2,
% then the validity of the answer can be checked, for example, by
% rerunning DNSQ with a tighter tolerance.
%
% Convergence Test. If DENORM(Z) denotes the Euclidean norm of a
% vector Z and D is the diagonal matrix whose entries are
% defined by the array DIAG, then this test attempts to
% guarantee that
%
% DENORM(D*(X-XSOL)) .LE. XTOL*DENORM(D*XSOL).
%
% If this condition is satisfied with XTOL = 10**(-K), then the
% larger components of D*X have K significant decimal digits and
% INFO is set to 1. There is a danger that the smaller
% components of D*X may have large relative errors, but the fast
% rate of convergence of DNSQ usually avoids this possibility.
% Unless high precision solutions are required, the recommended
% value for XTOL is the square root of the machine precision.
%
%
% 5. Unsuccessful Completion.
%
% Unsuccessful termination of DNSQ can be due to improper input
% parameters, arithmetic interrupts, an excessive number of
% function evaluations, or lack of good progress.
%
% Improper Input Parameters. INFO is set to 0 if IOPT .LT .1,
% or IOPT .GT. 2, or N .LE. 0, or LDFJAC .LT. N, or
% XTOL .LT. 0.0E0, or MAXFEV .LE. 0, or ML .LT. 0, or MU .LT. 0,
% or FACTOR .LE. 0.0E0, or LR .LT. (N*(N+1))/2.
%
% Arithmetic Interrupts. If these interrupts occur in the FCN
% subroutine during an early stage of the computation, they may
% be caused by an unacceptable choice of X by DNSQ. In this
% case, it may be possible to remedy the situation by rerunning
% DNSQ with a smaller value of FACTOR.
%
% Excessive Number of Function Evaluations. A reasonable value
% for MAXFEV is 100*(N+1) for IOPT=1 and 200*(N+1) for IOPT=2.
% If the number of calls to FCN reaches MAXFEV, then this
% indicates that the routine is converging very slowly as
% measured by the progress of FVEC, and INFO is set to 2. This
% situation should be unusual because, as indicated below, lack
% of good progress is usually diagnosed earlier by DNSQ,
% causing termination with info = 4 or INFO = 5.
%
% Lack of Good Progress. DNSQ searches for a zero of the system
% by minimizing the sum of the squares of the functions. In so
% doing, it can become trapped in a region where the minimum
% does not correspond to a zero of the system and, in this
% situation, the iteration eventually fails to make good
% progress. In particular, this will happen if the system does
% not have a zero. If the system has a zero, rerunning DNSQ
% from a different starting point may be helpful.
%
%
% 6. Characteristics of The Algorithm.
%
% DNSQ is a modification of the Powell Hybrid method. Two of its
% main characteristics involve the choice of the correction as a
% convex combination of the Newton and scaled gradient directions,
% and the updating of the Jacobian by the rank-1 method of
% Broyden. The choice of the correction guarantees (under
% reasonable conditions) global convergence for starting points
% far from the solution and a fast rate of convergence. The
% Jacobian is calculated at the starting point by either the
% user-supplied subroutine or a forward-difference approximation,
% but it is not recalculated until the rank-1 method fails to
% produce satisfactory progress.
%
% Timing. The time required by DNSQ to solve a given problem
% depends on N, the behavior of the functions, the accuracy
% requested, and the starting point. The number of arithmetic
% operations needed by DNSQ is about 11.5*(N**2) to process
% each evaluation of the functions (call to FCN) and 1.3*(N**3)
% to process each evaluation of the Jacobian (call to JAC,
% if IOPT = 1). Unless FCN and JAC can be evaluated quickly,
% the timing of DNSQ will be strongly influenced by the time
% spent in FCN and JAC.
%
% Storage. DNSQ requires (3*N**2 + 17*N)/2 single precision
% storage locations, in addition to the storage required by the
% program. There are no internally declared storage arrays.
%
% *Long Description:
%
% 7. Example.
%
% The problem is to determine the values of X(1), X(2), ..., X(9),
% which solve the system of tridiagonal equations
%
% (3-2*X(1))*X(1) -2*X(2) = -1
% -X(I-1) + (3-2*X(I))*X(I) -2*X(I+1) = -1, I=2-8
% -X(8) + (3-2*X(9))*X(9) = -1
% C **********
%
% program TEST
% C
% C Driver for DNSQ example.
% C
% INTEGER J,IOPT,N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR,
% * NWRITE
% doubleprecision XTOL,EPSFCN,FACTOR,FNORM
% doubleprecision X(9),FVEC(9),DIAG(9),FJAC(9,9),R(45),QTF(9),
% * WA1(9),WA2(9),WA3(9),WA4(9)
% doubleprecision DENORM,D1MACH
% EXTERNAL FCN
% DATA NWRITE /6/
% C
% IOPT = 2
% N = 9
% C
% C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH SOLUTION.
% C
% DO 10 J = 1, 9
% X(J) = -1.0E0
% 10 CONTINUE
% C
% LDFJAC = 9
% LR = 45
% C
% C SET XTOL TO THE SQUARE ROOT OF THE MACHINE PRECISION.
% C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED,
% C THIS IS THE RECOMMENDED SETTING.
% C
% XTOL = SQRT(D1MACH(4))
% C
% MAXFEV = 2000
% ML = 1
% MU = 1
% EPSFCN = 0.0E0
% MODE = 2
% DO 20 J = 1, 9
% DIAG(J) = 1.0E0
% 20 CONTINUE
% FACTOR = 1.0E2
% NPRINT = 0
% C
% CALL DNSQ(FCN,JAC,IOPT,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,ML,MU,
% * EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,
% * R,LR,QTF,WA1,WA2,WA3,WA4)
% FNORM = DENORM(N,FVEC)
% WRITE (NWRITE,1000) FNORM,NFEV,INFO,(X(J),J=1,N)
% STOP
% 1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
% * 5X,' NUMBER OF FUNCTION EVALUATIONS',I10 //
% * 5X,' EXIT PARAMETER',16X,I10 //
% * 5X,' FINAL APPROXIMATE SOLUTION' // (5X,3E15.7))
% end
% subroutine FCN(N,X,FVEC,IFLAG)
% INTEGER N,IFLAG
% doubleprecision X(N),FVEC(N)
% INTEGER K
% doubleprecision ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO
% DATA ZERO,ONE,TWO,THREE /0.0E0,1.0E0,2.0E0,3.0E0/
% C
% IF (IFLAG .NE. 0) GO TO 5
% C
% C INSERT PRINT STATEMENTS HERE WHEN NPRINT IS POSITIVE.
% C
% RETURN
% 5 CONTINUE
% DO 10 K = 1, N
% TEMP = (THREE - TWO*X(K))*X(K)
% TEMP1 = ZERO
% IF (K .NE. 1) TEMP1 = X(K-1)
% TEMP2 = ZERO
% IF (K .NE. N) TEMP2 = X(K+1)
% FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE
% 10 CONTINUE
% RETURN
% end
%
% Results obtained with different compilers or machines
% may be slightly different.
%
% Final L2 norm of the residuals 0.1192636E-07
%
% Number of function evaluations 14
%
% Exit parameter 1
%
% Final approximate solution
%
% -0.5706545E+00 -0.6816283E+00 -0.7017325E+00
% -0.7042129E+00 -0.7013690E+00 -0.6918656E+00
% -0.6657920E+00 -0.5960342E+00 -0.4164121E+00
%
%***REFERENCES M. J. D. Powell, A hybrid method for nonlinear equa-
% tions. In Numerical Methods for Nonlinear Algebraic
% Equations, P. Rabinowitz, Editor. Gordon and Breach,
% 1988.
%***ROUTINES CALLED D1MACH, D1MPYQ, D1UPDT, DDOGLG, DENORM, DFDJC1,
% DQFORM, DQRFAC, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800301 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 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)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE DNSQ
persistent actred delta epsmch firstCall fnorm fnorm1 gt i iflag iter iwa j jeval jm1 l ncfail ncsuc nslow1 nslow2 one p0001 p001 p1 p5 pnorm prered ratio sing summlv temp xnorm zero ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(iflag), iflag=0; end;
if isempty(iter), iter=0; end;
if isempty(iwa), iwa=zeros(1,1); end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(l), l=0; end;
if isempty(ncfail), ncfail=0; end;
if isempty(ncsuc), ncsuc=0; end;
if isempty(nslow1), nslow1=0; end;
if isempty(nslow2), nslow2=0; end;
if isempty(gt), gt=0; end;
if isempty(actred), actred=0; end;
if isempty(delta), delta=0; end;
diag_shape=size(diag);diag=reshape(diag,1,[]);
if isempty(epsmch), epsmch=0; end;
fjac_shape=size(fjac);fjac=reshape([fjac(:).',zeros(1,ceil(numel(fjac)./prod([ldfjac])).*prod([ldfjac])-numel(fjac))],ldfjac,[]);
if isempty(fnorm), fnorm=0; end;
if isempty(fnorm1), fnorm1=0; end;
fvec_shape=size(fvec);fvec=reshape(fvec,1,[]);
if isempty(one), one=0; end;
if isempty(p0001), p0001=0; end;
if isempty(p001), p001=0; end;
if isempty(p1), p1=0; end;
if isempty(p5), p5=0; end;
if isempty(pnorm), pnorm=0; end;
if isempty(prered), prered=0; end;
qtf_shape=size(qtf);qtf=reshape(qtf,1,[]);
r_shape=size(r);r=reshape(r,1,[]);
if isempty(ratio), ratio=0; end;
if isempty(summlv), summlv=0; end;
if isempty(temp), temp=0; end;
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,[]);
x_shape=size(x);x=reshape(x,1,[]);
if isempty(xnorm), xnorm=0; end;
if isempty(zero), zero=0; end;
if isempty(jeval), jeval=false; end;
if isempty(sing), sing=false; end;
if firstCall, one =[1.0d0]; end;
if firstCall, p1 =[1.0d-1]; end;
if firstCall, p5 =[5.0d-1]; end;
if firstCall, p001 =[1.0d-3]; end;
if firstCall, p0001 =[1.0d-4]; end;
if firstCall, zero=[0.0d0]; end;
firstCall=0;
%
% BEGIN BLOCK PERMITTING ...EXITS TO 320
%***FIRST EXECUTABLE STATEMENT DNSQ
[epsmch ]=d1mach(4);
%
info = 0;
iflag = 0;
nfev = 0;
njev = 0;
%
% CHECK THE INPUT PARAMETERS FOR ERRORS.
%
% ...EXIT
gt=0;
while (1);
if( iopt>=1 && iopt<=2 && n>0 && xtol>=zero &&maxfev>0 && ml>=0 && mu>=0 && factor>zero &&ldfjac>=n && lr>=fix((n.*(n+1))./2) )
if( mode==2 )
for j = 1 : n;
% .........EXIT
if( diag(j)<=zero )
gt=1;
break;
end;
end;
if(gt==1)
break;
end;
end;
%
% EVALUATE THE FUNCTION AT THE STARTING POINT
% AND CALCULATE ITS NORM.
%
iflag = 1;
[n,x,fvec,iflag]=fcn(n,x,fvec,iflag);
nfev = 1;
% ...EXIT
if( iflag>=0 )
[fnorm ,n,fvec]=denorm(n,fvec);
%
% INITIALIZE ITERATION COUNTER AND MONITORS.
%
iter = 1;
ncsuc = 0;
ncfail = 0;
nslow1 = 0;
nslow2 = 0;
%
% BEGINNING OF THE OUTER LOOP.
%
% BEGIN BLOCK PERMITTING ...EXITS TO 90
while( true );
jeval = true;
%
% CALCULATE THE JACOBIAN MATRIX.
%
if( iopt==2 )
%
% CODE APPROXIMATES THE JACOBIAN
%
iflag = 2;
[fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,wa2]=dfdjc1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,wa2);
nfev = fix(nfev + min(ml+mu+1,n));
else;
%
% USER SUPPLIES JACOBIAN
%
[n,x,fvec,fjac,ldfjac,iflag]=jac(n,x,fvec,fjac,ldfjac,iflag);
njev = fix(njev + 1);
end;
%
% .........EXIT
if( iflag<0 )
break;
end;
%
% COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
%
n_orig=n; [n,dumvar2,fjac,ldfjac,dumvar5,iwa,dumvar7,wa1,wa2,wa3]=dqrfac(n,n,fjac,ldfjac,false,iwa,1,wa1,wa2,wa3); n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);
%
% ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
% TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
%
% ...EXIT
if( iter==1 )
if( mode~=2 )
for j = 1 : n;
diag(j) = wa2(j);
if( wa2(j)==zero )
diag(j) = one;
end;
end; j = fix(n+1);
end;
%
% ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED
% X AND INITIALIZE THE STEP BOUND DELTA.
%
for j = 1 : n;
wa3(j) = diag(j).*x(j);
end; j = fix(n+1);
[xnorm ,n,wa3]=denorm(n,wa3);
delta = factor.*xnorm;
if( delta==zero )
delta = factor;
end;
end;
%
% FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.
%
for i = 1 : n;
qtf(i) = fvec(i);
end; i = fix(n+1);
for j = 1 : n;
if( fjac(j,j)~=zero )
summlv = zero;
for i = j : n;
summlv = summlv + fjac(i,j).*qtf(i);
end; i = fix(n+1);
temp = -summlv./fjac(j,j);
for i = j : n;
qtf(i) = qtf(i) + fjac(i,j).*temp;
end; i = fix(n+1);
end;
end; j = fix(n+1);
%
% COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.
%
sing = false;
for j = 1 : n;
l = fix(j);
jm1 = fix(j - 1);
if( jm1>=1 )
for i = 1 : jm1;
r(l) = fjac(i,j);
l = fix(l + n - i);
end; i = fix(jm1+1);
end;
r(l) = wa1(j);
if( wa1(j)==zero )
sing = true;
end;
end; j = fix(n+1);
%
% ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.
%
n_orig=n; [n,dumvar2,fjac,ldfjac,wa1]=dqform(n,n,fjac,ldfjac,wa1); n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);
%
% RESCALE IF NECESSARY.
%
if( mode~=2 )
for j = 1 : n;
diag(j) = max(diag(j),wa2(j));
end; j = fix(n+1);
end;
%
% BEGINNING OF THE INNER LOOP.
%
%
% IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
%
while( true );
if( nprint>0 )
iflag = 0;
if( rem(iter-1,nprint)==0 )
[n,x,fvec,iflag]=fcn(n,x,fvec,iflag);
end;
% ............EXIT
if( iflag<0 )
gt=1;
break;
end;
end;
%
% DETERMINE THE DIRECTION P.
%
[n,r,lr,diag,qtf,delta,wa1,wa2,wa3]=ddoglg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3);
%
% STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
%
for j = 1 : n;
wa1(j) = -wa1(j);
wa2(j) = x(j) + wa1(j);
wa3(j) = diag(j).*wa1(j);
end; j = fix(n+1);
[pnorm ,n,wa3]=denorm(n,wa3);
%
% ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
%
if( iter==1 )
delta = min(delta,pnorm);
end;
%
% EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
%
iflag = 1;
[n,wa2,wa4,iflag]=fcn(n,wa2,wa4,iflag);
nfev = fix(nfev + 1);
% .........EXIT
if( iflag<0 )
gt=1;
break;
end;
[fnorm1 ,n,wa4]=denorm(n,wa4);
%
% COMPUTE THE SCALED ACTUAL REDUCTION.
%
actred = -one;
if( fnorm1<fnorm )
actred = one -(fnorm1./fnorm).^2;
end;
%
% COMPUTE THE SCALED PREDICTED REDUCTION.
%
l = 1;
for i = 1 : n;
summlv = zero;
for j = i : n;
summlv = summlv + r(l).*wa1(j);
l = fix(l + 1);
end; j = fix(n+1);
wa3(i) = qtf(i) + summlv;
end; i = fix(n+1);
[temp ,n,wa3]=denorm(n,wa3);
prered = zero;
if( temp<fnorm )
prered = one -(temp./fnorm).^2;
end;
%
% COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
% REDUCTION.
%
ratio = zero;
if( prered>zero )
ratio = actred./prered;
end;
%
% UPDATE THE STEP BOUND.
%
if( ratio>=p1 )
ncfail = 0;
ncsuc = fix(ncsuc + 1);
if( ratio>=p5 || ncsuc>1 )
delta = max(delta,pnorm./p5);
end;
if( abs(ratio-one)<=p1 )
delta = pnorm./p5;
end;
else;
ncsuc = 0;
ncfail = fix(ncfail + 1);
delta = p5.*delta;
end;
%
% TEST FOR SUCCESSFUL ITERATION.
%
if( ratio>=p0001 )
%
% SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
%
for j = 1 : n;
x(j) = wa2(j);
wa2(j) = diag(j).*x(j);
fvec(j) = wa4(j);
end; j = fix(n+1);
[xnorm ,n,wa2]=denorm(n,wa2);
fnorm = fnorm1;
iter = fix(iter + 1);
end;
%
% DETERMINE THE PROGRESS OF THE ITERATION.
%
nslow1 = fix(nslow1 + 1);
if( actred>=p001 )
nslow1 = 0;
end;
if( jeval )
nslow2 = fix(nslow2 + 1);
end;
if( actred>=p1 )
nslow2 = 0;
end;
%
% TEST FOR CONVERGENCE.
%
if( delta<=xtol.*xnorm || fnorm==zero )
info = 1;
end;
% .........EXIT
if( info~=0 )
gt=1;
break;
end;
%
% TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
%
if( nfev>=maxfev )
info = 2;
end;
if( p1.*max(p1.*delta,pnorm)<=epsmch.*xnorm )
info = 3;
end;
if( nslow2==5 )
info = 4;
end;
if( nslow1==10 )
info = 5;
end;
% .........EXIT
if( info~=0 )
gt=1;
break;
end;
%
% CRITERION FOR RECALCULATING JACOBIAN
%
% ...EXIT
if( ncfail==2 )
break;
end;
%
% CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN
% AND UPDATE QTF IF NECESSARY.
%
for j = 1 : n;
summlv = zero;
for i = 1 : n;
summlv = summlv + fjac(i,j).*wa4(i);
end; i = fix(n+1);
wa2(j) =(summlv-wa3(j))./pnorm;
wa1(j) = diag(j).*((diag(j).*wa1(j))./pnorm);
if( ratio>=p0001 )
qtf(j) = summlv;
end;
end; j = fix(n+1);
%
% COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.
%
n_orig=n; [n,dumvar2,r,lr,wa1,wa2,wa3,sing]=d1updt(n,n,r,lr,wa1,wa2,wa3,sing); n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);
n_orig=n; [n,dumvar2,fjac,ldfjac,wa2,wa3]=d1mpyq(n,n,fjac,ldfjac,wa2,wa3); n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);
[dumvar1,n,qtf,dumvar4,wa2,wa3]=d1mpyq(1,n,qtf,1,wa2,wa3);
%
% end OF THE INNER LOOP.
%
jeval = false;
%
% end OF THE OUTER LOOP.
%
end;
if(gt==1)
break;
end;
end;
if(gt==1)
break;
end;
end;
end;
break;
end;
%
% TERMINATION, EITHER NORMAL OR USER IMPOSED.
%
if( iflag<0 )
info = fix(iflag);
end;
iflag = 0;
if( nprint>0 )
[n,x,fvec,iflag]=fcn(n,x,fvec,iflag);
end;
if( info<0 )
xermsg('SLATEC','DNSQ','EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.',1,1);
end;
if( info==0 )
xermsg('SLATEC','DNSQ','INVALID INPUT PARAMETER.',2,1);
end;
if( info==2 )
xermsg('SLATEC','DNSQ','TOO MANY FUNCTION EVALUATIONS.',9,1);
end;
if( info==3 )
xermsg('SLATEC','DNSQ','XTOL TOO SMALL. NO FURTHER IMPROVEMENT POSSIBLE.',3,1);
end;
if( info>4 )
xermsg('SLATEC','DNSQ','ITERATION NOT MAKING GOOD PROGRESS.',1,1);
end;
%
% LAST CARD OF SUBROUTINE DNSQ.
%
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
fjac_shape=zeros(fjac_shape);fjac_shape(:)=fjac(1:numel(fjac_shape));fjac=fjac_shape;
fvec_shape=zeros(fvec_shape);fvec_shape(:)=fvec(1:numel(fvec_shape));fvec=fvec_shape;
qtf_shape=zeros(qtf_shape);qtf_shape(:)=qtf(1:numel(qtf_shape));qtf=qtf_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_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;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end %subroutine dnsq
%DECK DOGLEG
|
|
