| [fcn,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,epsfcn,diag,mode,factor,nprint,info,nfev,njev,ipvt,qtf,wa1,wa2,wa3,wa4]=dnls1(fcn,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,epsfcn,diag,mode,factor,nprint,info,nfev,njev,ipvt,qtf,wa1,wa2,wa3,w |
function [fcn,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,epsfcn,diag,mode,factor,nprint,info,nfev,njev,ipvt,qtf,wa1,wa2,wa3,wa4]=dnls1(fcn,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,epsfcn,diag,mode,factor,nprint,info,nfev,njev,ipvt,qtf,wa1,wa2,wa3,wa4);
% C
% X(1) = 1.0E0
% X(2) = 1.0E0
% X(3) = 1.0E0
% C
% LDFJAC = 15
% C
% C Set FTOL and XTOL to the square root of the machine precision
% C and GTOL to zero. Unless high precision solutions are
% C required, these are the recommended settings.
% C
% FTOL = SQRT(R1MACH(4))
% XTOL = SQRT(R1MACH(4))
% GTOL = 0.0E0
% C
% MAXFEV = 400
% EPSFCN = 0.0
% MODE = 1
% FACTOR = 1.0E2
% NPRINT = 0
% C
% CALL DNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,
% * GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,
% * INFO,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4)
% FNORM = ENORM(M,FVEC)
% WRITE (NWRITE,1000) FNORM,NFEV,NJEV,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,' NUMBER OF JACOBIAN EVALUATIONS',I10 //
% * 5X,' EXIT PARAMETER',16X,I10 //
% * 5X,' FINAL APPROXIMATE SOLUTION' // 5X,3E15.7)
% end
% subroutine FCN(IFLAG,M,N,X,FVEC,DUM,IDUM)
% C This is the form of the FCN routine if IOPT=1,
% C that is, if the user does not calculate the Jacobian.
% INTEGER I,M,N,IFLAG
% doubleprecision X(N),FVEC(M),Y(15)
% doubleprecision TMP1,TMP2,TMP3,TMP4
% DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
% * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
% * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
% * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
% C
% IF (IFLAG .NE. 0) GO TO 5
% C
% C Insert print statements here when NPRINT is positive.
% C
% RETURN
% 5 CONTINUE
% DO 10 I = 1, M
% TMP1 = I
% TMP2 = 16 - I
% TMP3 = TMP1
% IF (I .GT. 8) TMP3 = TMP2
% FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
% 10 CONTINUE
% RETURN
% end
%
%
% Results obtained with different compilers or machines
% may be slightly different.
%
% FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01
%
% NUMBER OF FUNCTION EVALUATIONS 25
%
% NUMBER OF JACOBIAN EVALUATIONS 0
%
% EXIT PARAMETER 1
%
% FINAL APPROXIMATE SOLUTION
%
% 0.8241058E-01 0.1133037E+01 0.2343695E+01
%
%
% For IOPT=2, FCN would be modified as follows to also
% calculate the full Jacobian when IFLAG=2.
%
% subroutine FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
% C
% C This is the form of the FCN routine if IOPT=2,
% C that is, if the user calculates the full Jacobian.
% C
% INTEGER I,LDFJAC,M,N,IFLAG
% doubleprecision X(N),FVEC(M),FJAC(LDFJAC,N),Y(15)
% doubleprecision TMP1,TMP2,TMP3,TMP4
% DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
% * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
% * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
% * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
% C
% IF (IFLAG .NE. 0) GO TO 5
% C
% C Insert print statements here when NPRINT is positive.
% C
% RETURN
% 5 CONTINUE
% IF(IFLAG.NE.1) GO TO 20
% DO 10 I = 1, M
% TMP1 = I
% TMP2 = 16 - I
% TMP3 = TMP1
% IF (I .GT. 8) TMP3 = TMP2
% FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
% 10 CONTINUE
% RETURN
% C
% C Below, calculate the full Jacobian.
% C
% 20 CONTINUE
% C
% DO 30 I = 1, M
% TMP1 = I
% TMP2 = 16 - I
% TMP3 = TMP1
% IF (I .GT. 8) TMP3 = TMP2
% TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
% FJAC(I,1) = -1.0E0
% FJAC(I,2) = TMP1*TMP2/TMP4
% FJAC(I,3) = TMP1*TMP3/TMP4
% 30 CONTINUE
% RETURN
% end
%
%
% For IOPT = 3, FJAC would be dimensioned as FJAC(3,3),
% LDFJAC would be set to 3, and FCN would be written as
% follows to calculate a row of the Jacobian when IFLAG=3.
%
% subroutine FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
% C This is the form of the FCN routine if IOPT=3,
% C that is, if the user calculates the Jacobian row by row.
% INTEGER I,M,N,IFLAG
% doubleprecision X(N),FVEC(M),FJAC(N),Y(15)
% doubleprecision TMP1,TMP2,TMP3,TMP4
% DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
% * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
% * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
% * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
% C
% IF (IFLAG .NE. 0) GO TO 5
% C
% C Insert print statements here when NPRINT is positive.
% C
% RETURN
% 5 CONTINUE
% IF( IFLAG.NE.1) GO TO 20
% DO 10 I = 1, M
% TMP1 = I
% TMP2 = 16 - I
% TMP3 = TMP1
% IF (I .GT. 8) TMP3 = TMP2
% FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
% 10 CONTINUE
% RETURN
% C
% C Below, calculate the LDFJAC-th row of the Jacobian.
% C
% 20 CONTINUE
%
% I = LDFJAC
% TMP1 = I
% TMP2 = 16 - I
% TMP3 = TMP1
% IF (I .GT. 8) TMP3 = TMP2
% TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
% FJAC(1) = -1.0E0
% FJAC(2) = TMP1*TMP2/TMP4
% FJAC(3) = TMP1*TMP3/TMP4
% RETURN
% end
%
%***REFERENCES Jorge J. More, The Levenberg-Marquardt algorithm:
% implementation and theory. In Numerical Analysis
% Proceedings (Dundee, June 28 - July 1, 1977, G. A.
% Watson, Editor), Lecture Notes in Mathematics 630,
% Springer-Verlag, 1978.
%***ROUTINES CALLED D1MACH, DCKDER, DENORM, DFDJC3, DMPAR, DQRFAC,
% DWUPDT, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800301 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 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 Convert XERRWV calls to XERMSG calls. (RWC)
% 920205 Corrected XERN1 declaration. (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE DNLS1
persistent actred chklim delta dirder epsmch err firstCall fnorm fnorm1 gnorm gt i iflag ijunk iter j l modech nrow one p0001 p1 p25 p5 p75 par pnorm prered ratio sing summlv temp temp1 temp2 xern1 xern3 xnorm zero ; if isempty(firstCall),firstCall=1;end;
if isempty(gt), gt=0; end;
if isempty(ijunk), ijunk=0; end;
if isempty(nrow), nrow=0; end;
ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
fvec_shape=size(fvec);fvec=reshape(fvec,1,[]);
fjac_shape=size(fjac);fjac=reshape([fjac(:).',zeros(1,ceil(numel(fjac)./prod([ldfjac])).*prod([ldfjac])-numel(fjac))],ldfjac,[]);
diag_shape=size(diag);diag=reshape(diag,1,[]);
qtf_shape=size(qtf);qtf=reshape(qtf,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(sing), sing=false; end;
if isempty(i), i=0; end;
if isempty(iflag), iflag=0; end;
if isempty(iter), iter=0; end;
if isempty(j), j=0; end;
if isempty(l), l=0; end;
if isempty(modech), modech=0; end;
if isempty(actred), actred=0; end;
if isempty(delta), delta=0; end;
if isempty(dirder), dirder=0; end;
if isempty(epsmch), epsmch=0; end;
if isempty(fnorm), fnorm=0; end;
if isempty(fnorm1), fnorm1=0; end;
if isempty(gnorm), gnorm=0; end;
if isempty(one), one=0; end;
if isempty(par), par=0; end;
if isempty(pnorm), pnorm=0; end;
if isempty(prered), prered=0; end;
if isempty(p1), p1=0; end;
if isempty(p5), p5=0; end;
if isempty(p25), p25=0; end;
if isempty(p75), p75=0; end;
if isempty(p0001), p0001=0; end;
if isempty(ratio), ratio=0; end;
if isempty(summlv), summlv=0; end;
if isempty(temp), temp=0; end;
if isempty(temp1), temp1=0; end;
if isempty(temp2), temp2=0; end;
if isempty(xnorm), xnorm=0; end;
if isempty(zero), zero=0; end;
if isempty(err), err=0; end;
if isempty(chklim), chklim=0; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,16); end;
%
if firstCall, chklim=[.1d0]; end;
if firstCall, one =[1.0d0]; end;
if firstCall, p1 =[1.0d-1]; end;
if firstCall, p5 =[5.0d-1]; end;
if firstCall, p25 =[2.5d-1]; end;
if firstCall, p75 =[7.5d-1]; end;
if firstCall, p0001 =[1.0d-4]; end;
if firstCall, zero=[0.0d0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DNLS1
[epsmch ]=d1mach(4);
%
info = 0;
iflag = 0;
nfev = 0;
njev = 0;
%
% CHECK THE INPUT PARAMETERS FOR ERRORS.
%
while (1);
if( iopt>=1 && iopt<=3 && n>0 && m>=n &&ldfjac>=n && ftol>=zero && xtol>=zero &>ol>=zero && maxfev>0 && factor>zero )
if( iopt>=3 || ldfjac>=m )
if( mode==2 )
for j = 1 : n;
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;
ijunk = 1;
[iflag,m,n,x,fvec,fjac,ijunk]=fcn(iflag,m,n,x,fvec,fjac,ijunk);
nfev = 1;
if( iflag>=0 )
[fnorm ,m,fvec]=denorm(m,fvec);
%
% INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER.
%
par = zero;
iter = 1;
%
% BEGINNING OF THE OUTER LOOP.
%
%
% IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
%
while( true );
if( nprint>0 )
iflag = 0;
if( rem(iter-1,nprint)==0 )
[iflag,m,n,x,fvec,fjac,ijunk]=fcn(iflag,m,n,x,fvec,fjac,ijunk);
end;
if( iflag<0 )
break;
end;
end;
%
% CALCULATE THE JACOBIAN MATRIX.
%
if( iopt==3 )
%
% ACCUMULATE THE JACOBIAN BY ROWS IN ORDER TO SAVE STORAGE.
% COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX
% CALCULATED ONE ROW AT A TIME, WHILE SIMULTANEOUSLY
% FORMING (Q TRANSPOSE)*FVEC AND STORING THE FIRST
% N COMPONENTS IN QTF.
%
for j = 1 : n;
qtf(j) = zero;
for i = 1 : n;
fjac(i,j) = zero;
end; i = fix(n+1);
end; j = fix(n+1);
for i = 1 : m;
nrow = fix(i);
iflag = 3;
[iflag,m,n,x,fvec,wa3,nrow]=fcn(iflag,m,n,x,fvec,wa3,nrow);
if( iflag<0 )
gt=1;
break;
end;
%
% ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN.
%
if( iter<=1 )
%
% GET THE INCREMENTED X-VALUES INTO WA1(*).
%
modech = 1;
[m,n,x,fvec,fjac,ldfjac,wa1,wa4,modech,err]=dckder(m,n,x,fvec,fjac,ldfjac,wa1,wa4,modech,err);
%
% EVALUATE AT INCREMENTED VALUES, IF NOT ALREADY EVALUATED.
%
if( i==1 )
%
% EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT INTO WA4(*).
%
iflag = 1;
[iflag,m,n,wa1,wa4,fjac,nrow]=fcn(iflag,m,n,wa1,wa4,fjac,nrow);
nfev = fix(nfev + 1);
if( iflag<0 )
gt=1;
break;
end;
end;
modech = 2;
[dumvar1,n,x,fvec(sub2ind(size(fvec),max(i,1)):end),wa3,dumvar6,wa1,wa4(sub2ind(size(wa4),max(i,1)):end),modech,err]=dckder(1,n,x,fvec(sub2ind(size(fvec),max(i,1)):end),wa3,1,wa1,wa4(sub2ind(size(wa4),max(i,1)):end),modech,err);
if( err<chklim )
xern1=sprintf(['%8i'], i);
xern3=sprintf([repmat('%15.6f',1,1)], err);
xermsg('SLATEC','DNLS1',['DERIVATIVE OF FUNCTION ',[xern1,[' MAY BE WRONG, ERR = ',[xern3,' TOO CLOSE TO 0.']]]],7,0);
end;
end;
%
temp = fvec(i);
[n,fjac,ldfjac,wa3,qtf,temp,wa1,wa2]=dwupdt(n,fjac,ldfjac,wa3,qtf,temp,wa1,wa2);
end;
if(gt==1)
break;
end;
njev = fix(njev + 1);
%
% IF THE JACOBIAN IS RANK DEFICIENT, CALL DQRFAC TO
% REORDER ITS COLUMNS AND UPDATE THE COMPONENTS OF QTF.
%
sing = false;
for j = 1 : n;
if( fjac(j,j)==zero )
sing = true;
end;
ipvt(j) = fix(j);
[wa2(j) ,j,fjac(sub2ind(size(fjac),1,j):end)]=denorm(j,fjac(sub2ind(size(fjac),1,j):end));
end; j = fix(n+1);
if( sing )
n_orig=n; n_orig=n; [n,dumvar2,fjac,ldfjac,dumvar5,ipvt,dumvar7,wa1,wa2,wa3]=dqrfac(n,n,fjac,ldfjac,true,ipvt,n,wa1,wa2,wa3); n(dumvar7~=n_orig)=dumvar7(dumvar7~=n_orig); n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);
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;
fjac(j,j) = wa1(j);
end; j = fix(n+1);
end;
else;
%
% STORE THE FULL JACOBIAN USING M*N STORAGE
%
if( iopt==1 )
%
% THE CODE APPROXIMATES THE JACOBIAN
%
iflag = 1;
[fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4]=dfdjc3(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4);
nfev = fix(nfev + n);
else;
%
% THE USER SUPPLIES THE JACOBIAN
%
iflag = 2;
[iflag,m,n,x,fvec,fjac,ldfjac]=fcn(iflag,m,n,x,fvec,fjac,ldfjac);
njev = fix(njev + 1);
%
% ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN
%
if( iter<=1 )
if( iflag<0 )
break;
end;
%
% GET THE INCREMENTED X-VALUES INTO WA1(*).
%
modech = 1;
[m,n,x,fvec,fjac,ldfjac,wa1,wa4,modech,err]=dckder(m,n,x,fvec,fjac,ldfjac,wa1,wa4,modech,err);
%
% EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT IN WA4(*).
%
iflag = 1;
[iflag,m,n,wa1,wa4,fjac,ldfjac]=fcn(iflag,m,n,wa1,wa4,fjac,ldfjac);
nfev = fix(nfev + 1);
if( iflag<0 )
break;
end;
for i = 1 : m;
%
modech = 2;
[dumvar1,n,x,fvec(sub2ind(size(fvec),max(i,1)):end),fjac(sub2ind(size(fjac),i,1):end),ldfjac,wa1,wa4(sub2ind(size(wa4),max(i,1)):end),modech,err]=dckder(1,n,x,fvec(sub2ind(size(fvec),max(i,1)):end),fjac(sub2ind(size(fjac),i,1):end),ldfjac,wa1,wa4(sub2ind(size(wa4),max(i,1)):end),modech,err);
if( err<chklim )
xern1=sprintf(['%8i'], i);
xern3=sprintf([repmat('%15.6f',1,1)], err);
xermsg('SLATEC','DNLS1',['DERIVATIVE OF ',['FUNCTION ',[xern1,[' MAY BE WRONG, ERR = ',[xern3,' TOO CLOSE TO 0.']]]]],7,0);
end;
end; i = fix(m+1);
end;
end;
if( iflag<0 )
break;
end;
%
% COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
%
n_orig=n; [m,n,fjac,ldfjac,dumvar5,ipvt,dumvar7,wa1,wa2,wa3]=dqrfac(m,n,fjac,ldfjac,true,ipvt,n,wa1,wa2,wa3); n(dumvar7~=n_orig)=dumvar7(dumvar7~=n_orig);
%
% FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN
% QTF.
%
for i = 1 : m;
wa4(i) = fvec(i);
end; i = fix(m+1);
for j = 1 : n;
if( fjac(j,j)~=zero )
summlv = zero;
for i = j : m;
summlv = summlv + fjac(i,j).*wa4(i);
end; i = fix(m+1);
temp = -summlv./fjac(j,j);
for i = j : m;
wa4(i) = wa4(i) + fjac(i,j).*temp;
end; i = fix(m+1);
end;
fjac(j,j) = wa1(j);
qtf(j) = wa4(j);
end; j = fix(n+1);
end;
%
% ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
% TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
%
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;
%
% COMPUTE THE NORM OF THE SCALED GRADIENT.
%
gnorm = zero;
if( fnorm~=zero )
for j = 1 : n;
l = fix(ipvt(j));
if( wa2(l)~=zero )
summlv = zero;
for i = 1 : j;
summlv = summlv + fjac(i,j).*(qtf(i)./fnorm);
end; i = fix(j+1);
gnorm = max(gnorm,abs(summlv./wa2(l)));
end;
end; j = fix(n+1);
end;
%
% TEST FOR CONVERGENCE OF THE GRADIENT NORM.
%
if( gnorm<=gtol )
info = 4;
end;
if( info~=0 )
break;
end;
%
% 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.
%
%
% DETERMINE THE LEVENBERG-MARQUARDT PARAMETER.
%
while( true );
[n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,wa3,wa4]=dmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,wa3,wa4);
%
% 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;
[iflag,m,n,wa2,wa4,fjac,ijunk]=fcn(iflag,m,n,wa2,wa4,fjac,ijunk);
nfev = fix(nfev + 1);
if( iflag<0 )
gt=1;
break;
end;
[fnorm1 ,m,wa4]=denorm(m,wa4);
%
% COMPUTE THE SCALED ACTUAL REDUCTION.
%
actred = -one;
if( p1.*fnorm1<fnorm )
actred = one -(fnorm1./fnorm).^2;
end;
%
% COMPUTE THE SCALED PREDICTED REDUCTION AND
% THE SCALED DIRECTIONAL DERIVATIVE.
%
for j = 1 : n;
wa3(j) = zero;
l = fix(ipvt(j));
temp = wa1(l);
for i = 1 : j;
wa3(i) = wa3(i) + fjac(i,j).*temp;
end; i = fix(j+1);
end; j = fix(n+1);
temp1 = denorm(n,wa3)./fnorm;
temp2 =(sqrt(par).*pnorm)./fnorm;
prered = temp1.^2 + temp2.^2./p5;
dirder = -(temp1.^2+temp2.^2);
%
% 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<=p25 )
if( actred>=zero )
temp = p5;
end;
if( actred<zero )
temp = p5.*dirder./(dirder+p5.*actred);
end;
if( p1.*fnorm1>=fnorm || temp<p1 )
temp = p1;
end;
delta = temp.*min(delta,pnorm./p1);
par = par./temp;
elseif( par==zero || ratio>=p75 ) ;
delta = pnorm./p5;
par = p5.*par;
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);
end; j = fix(n+1);
for i = 1 : m;
fvec(i) = wa4(i);
end; i = fix(m+1);
[xnorm ,n,wa2]=denorm(n,wa2);
fnorm = fnorm1;
iter = fix(iter + 1);
end;
%
% TESTS FOR CONVERGENCE.
%
if( abs(actred)<=ftol && prered<=ftol &&p5.*ratio<=one )
info = 1;
end;
if( delta<=xtol.*xnorm )
info = 2;
end;
if( abs(actred)<=ftol && prered<=ftol &&p5.*ratio<=one && info==2 )
info = 3;
end;
if( info~=0 )
gt=1;
break;
end;
%
% TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
%
if( nfev>=maxfev )
info = 5;
end;
if( abs(actred)<=epsmch && prered<=epsmch &&p5.*ratio<=one )
info = 6;
end;
if( delta<=epsmch.*xnorm )
info = 7;
end;
if( gnorm<=epsmch )
info = 8;
end;
if( info~=0 )
gt=1;
break;
end;
%
% end OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL.
%
%
% end OF THE OUTER LOOP.
%
if( ratio>=p0001 )
break;
end;
end;
if(gt==1)
break;
end;
end;
if(gt==1)
break;
end;
end;
end;
end;
break;
end;
%
% TERMINATION, EITHER NORMAL OR USER IMPOSED.
%
if( iflag<0 )
info = fix(iflag);
end;
iflag = 0;
if( nprint>0 )
[iflag,m,n,x,fvec,fjac,ijunk]=fcn(iflag,m,n,x,fvec,fjac,ijunk);
end;
if( info<0 )
xermsg('SLATEC','DNLS1','EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.',1,1);
end;
if( info==0 )
xermsg('SLATEC','DNLS1','INVALID INPUT PARAMETER.',2,1);
end;
if( info==4 )
xermsg('SLATEC','DNLS1','THIRD CONVERGENCE CONDITION, CHECK RESULTS BEFORE ACCEPTING.',1,1);
end;
if( info==5 )
xermsg('SLATEC','DNLS1','TOO MANY FUNCTION EVALUATIONS.',9,1);
end;
if( info>=6 )
xermsg('SLATEC','DNLS1','TOLERANCES TOO SMALL, NO FURTHER IMPROVEMENT POSSIBLE.',3,1);
end;
%
% LAST CARD OF SUBROUTINE DNLS1.
%
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
fvec_shape=zeros(fvec_shape);fvec_shape(:)=fvec(1:numel(fvec_shape));fvec=fvec_shape;
fjac_shape=zeros(fjac_shape);fjac_shape(:)=fjac(1:numel(fjac_shape));fjac=fjac_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
qtf_shape=zeros(qtf_shape);qtf_shape(:)=qtf(1:numel(qtf_shape));qtf=qtf_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 dnls1
%DECK DNRM2
|
|
