C ALGORITHM 768, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 232,NO. 21, June, 1997, P. 174--195. C #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Src # This archive created: Mon Sep 1 10:20:24 1997 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'makefile' then echo shar: will not over-write existing file "'makefile'" else cat << \SHAR_EOF > 'makefile' # # TENSOLVE directory # FFLAGS = -O -u # Files for TENSOLVE # The program runs nonlinear equations and nonlinear # least squares problems EXAMPLE1 = driver1.o EXAMPLE2 = driver2.o TENSOLVE = tensolve.o UNCMIN = uncmin.o BLAS = blas.o ex1 : $(EXAMPLE1) $(TENSOLVE) $(UNCMIN) $(BLAS) f77 $(FFLAGS) $(EXAMPLE1) $(TENSOLVE) $(UNCMIN) \ $(BLAS) -o tensolve ex2 : $(EXAMPLE2) $(TENSOLVE) $(UNCMIN) $(BLAS) f77 $(FFLAGS) $(EXAMPLE2) $(TENSOLVE) $(UNCMIN) \ $(BLAS) -o tensolve SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'RES1' then echo shar: will not over-write existing file "'RES1'" else cat << \SHAR_EOF > 'RES1' TSNESV TYPICAL X TSNESV 0.1000000000000D+01 0.1000000000000D+01 TSNESV DIAGONAL SCALING MATRIX FOR X TSNESV 0.1000000000000D+01 0.1000000000000D+01 TSNESV TYPICAL F TSNESV 0.1000000000000D+01 0.1000000000000D+01 TSNESV DIAGONAL SCALING MATRIX FOR F TSNESV 0.1000000000000D+01 0.1000000000000D+01 TSNESV JACOBIAN FLAG = 0 TSNESV METHOD USED = 1 TSNESV GLOBAL STRATEGY = 0 TSNESV ITERATION LIMIT = 150 TSNESV MACHINE EPSILON = 0.2220446049250D-15 TSNESV STEP TOLERANCE = 0.3666852862501D-10 TSNESV GRADIENT TOLERANCE = 0.6055454452393D-05 TSNESV FUNCTION TOLERANCE = 0.3666852862501D-10 TSNESV MAXIMUM STEP SIZE = 0.1000000000000D+04 TSNESV TRUST REG RADIUS =-0.1000000000000D+01 TSRSLT ITERATION K = 0 TSRSLT X(K) TSRSLT -0.1200000000000D+01 0.1000000000000D+01 TSRSLT FUNCTION AT X(K) TSRSLT 0.1210000000000D+02 TSRSLT GRADIENT AT X(K) TSRSLT -0.1077999998579D+03 -0.4400000000000D+02 TSNSTP RELATIVE GRADIENT CLOSE TO ZERO TSRSLT ITERATION K = 7 TSRSLT X(K) TSRSLT 0.1000000000353D+01 0.1000000000710D+01 TSRSLT FUNCTION AT X(K) TSRSLT 0.6294894734552D-19 TSRSLT GRADIENT AT X(K) TSRSLT -0.3079716530014D-09 0.3306244164987D-09 SHAR_EOF fi # end of overwriting check if test -f 'RES2' then echo shar: will not over-write existing file "'RES2'" else cat << \SHAR_EOF > 'RES2' TSNESV TYPICAL X TSNESV 0.1000000000000D+01 0.1000000000000D+01 0.1000000000000D+01 TSNESV 0.1000000000000D+01 TSNESV DIAGONAL SCALING MATRIX FOR X TSNESV 0.1000000000000D+01 0.1000000000000D+01 0.1000000000000D+01 TSNESV 0.1000000000000D+01 TSNESV TYPICAL F TSNESV 0.1000000000000D+01 0.1000000000000D+01 0.1000000000000D+01 TSNESV 0.1000000000000D+01 0.1000000000000D+01 0.1000000000000D+01 TSNESV TSNESV DIAGONAL SCALING MATRIX FOR F TSNESV 0.1000000000000D+01 0.1000000000000D+01 0.1000000000000D+01 TSNESV 0.1000000000000D+01 0.1000000000000D+01 0.1000000000000D+01 TSNESV TSNESV JACOBIAN FLAG = 0 TSNESV METHOD USED = 1 TSNESV GLOBAL STRATEGY = 1 TSNESV ITERATION LIMIT = 150 TSNESV MACHINE EPSILON = 0.2220446049250D-15 TSNESV STEP TOLERANCE = 0.1000000000000D-08 TSNESV GRADIENT TOLERANCE = 0.1000000000000D-04 TSNESV FUNCTION TOLERANCE = 0.1000000000000D-08 TSNESV MAXIMUM STEP SIZE = 0.1000000000000D+04 TSNESV TRUST REG RADIUS =-0.1000000000000D+01 TSRSLT ITERATION K = 0 TSRSLT X(K) TSRSLT -0.3000000000000D+03 -0.1000000000000D+03 -0.3000000000000D+03 TSRSLT -0.1000000000000D+03 TSRSLT FUNCTION AT X(K) TSRSLT 0.7712112446210D+12 TSRSLT GRADIENT AT X(K) TSRSLT -0.5406000277536D+10 -0.9012019999964D+07 -0.4865400281616D+10 TSRSLT -0.8110989138290D+07 TSNSTP RELATIVE GRADIENT CLOSE TO ZERO TSRSLT ITERATION K = 5 TSRSLT X(K) TSRSLT 0.1000000000045D+01 0.9999999998915D+00 0.9999999999587D+00 TSRSLT 0.1000000000108D+01 TSRSLT FUNCTION AT X(K) TSRSLT 0.3620445253087D-17 TSRSLT GRADIENT AT X(K) TSRSLT 0.3980304115768D-07 -0.1990154450127D-07 -0.3440781877366D-07 TSRSLT 0.1720404671870D-07 SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << \SHAR_EOF > 'driver2.f' program example2 c TENSOLVE finds roots of systems of n nonlinear equations in c n unknowns, or minimizers of the sum of squares of m > n c nonlinear functions in n unknowns, using tensor methods. c c This example illustrates the use of TENSOLVE to solve a c nonlinear least-squares problem defined by subroutines c fwood and jwood (included below). integer maxm, maxn, maxp, m, n, itnlim, jacflg, method, + global, ipr, lunc, lnem, lnen, lin, msg, termcd double precision gradtl, steptl, ftol, stepmx, dlt parameter (maxm = 100, maxn = 30, maxp = 6) parameter (lin = 3, lunc = 14, lnem = 51, lnen = 19) integer iwrkn(maxn,lin) double precision x0(maxn), xp(maxn), fp(maxm), gp(maxn), + typx(maxn), typf(maxm), + wrknen(maxn,lnen), wrkunc(maxp,lunc), + wrknem(maxm,lnem) external fwood, jwood c Set dimensions of the problem. m = 6 n = 4 c Set values for the initial point. x0(1) = -300.0d0 x0(2) = -100.0d0 x0(3) = -300.0d0 x0(4) = -100.0d0 c Set default values for the TENSOLVE parameters. call tsdflt(m , n , itnlim, jacflg, gradtl, steptl, + ftol , method, global, stepmx, dlt , + typx , typf , ipr , msg ) c Alter some of the parameters. gradtl = 1.0d-5 ftol = 1.0d-9 steptl = 1.0d-9 global = 1 c Call TENSOLVE. call tsneci(maxm , maxn , maxp , x0 , m , n , + typx , typf , itnlim, jacflg, gradtl, steptl, + ftol , method, global, stepmx, dlt , ipr , + wrkunc, lunc , wrknem, lnem , wrknen, lnen , + iwrkn , lin , fwood , jwood , msg , + xp , fp , gp , termcd ) c end of example2 main program end subroutine fwood ( x, f, m, n ) integer m, n double precision x(n), f(m) c fwood defines function values for the Wood function. f(1) = 10.0d0 * (x(2) - x(1)**2) f(2) = 1.0d0 - x(1) f(3) = sqrt(90.0d0) * (x(4) - x(3)**2) f(4) = 1.0d0 - x(3) f(5) = (x(2) + x(4) - 2.0d0) * sqrt(10.0d0) f(6) = (x(2) - x(4)) / sqrt(10.0d0) c end of fwood. end subroutine jwood ( x, jac, maxm, m, n ) integer maxm, m, n double precision x(n), jac(maxm,n) c jwood defines Jacobian values for the Wood function. integer i, j double precision tval do 20 j = 1, n do 10 i = 1, m jac(i,j) = 0.0d0 10 continue 20 continue jac(1,1) = -20.0d0 * x(1) jac(1,2) = 10.0d0 jac(2,1) = -1.0d0 tval = sqrt(90.0d0) jac(3,3) = -2.0d0 * tval * x(3) jac(3,4) = tval jac(4,3) = -1.0d0 tval = sqrt(10.0d0) jac(5,2) = tval jac(5,4) = tval tval = 1.0d0/tval jac(6,2) = tval jac(6,4) = -tval c end of jwood. end SHAR_EOF fi # end of overwriting check if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << \SHAR_EOF > 'driver1.f' program example1 c TENSOLVE finds roots of systems of n nonlinear equations in c n unknowns, or minimizers of the sum of squares of m > n c nonlinear functions in n unknowns, using tensor methods. c c This example illustrates the use of TENSOLVE to solve a c nonlinear equation problem defined by the subroutine c frosen (included below). integer maxm, maxn, maxp, m, n, lunc, lnem, lnen, + lin, msg, termcd parameter (maxm = 100, maxn = 30, maxp = 6) parameter (lin = 3, lunc = 14, lnem = 51, lnen = 19) integer iwrkn(maxn,lin) double precision x0(maxn), xp(maxn), fp(maxm), gp(maxn), + wrknen(maxn,lnen), wrkunc(maxp,lunc), + wrknem(maxm,lnem) external frosen c Set dimensions of the problem. m = 2 n = 2 c Set values for the initial point. x0(1) = -1.20d0 x0(2) = 1.0d0 c Call TENSOLVE. call tsnesi(maxm , maxn, maxp , x0 , m , n , + wrkunc, lunc, wrknem, lnem, wrknen, lnen, + iwrkn , lin , frosen, msg , + xp , fp , gp , termcd ) c end of example1 main program. end subroutine frosen ( x, f, m, n ) integer m, n double precision x(n), f(m) c frosen defines function values for the Rosenbrock function. f(1) = 10.0d0 * (x(2) - x(1)**2) f(2) = 1.0d0 - x(1) c end of frosen. end SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'blas.f' then echo shar: will not over-write existing file "'blas.f'" else cat << \SHAR_EOF > 'blas.f' DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*),DTEMP INTEGER I,INCX,M,MP1,N,NINCX C DASUM = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DTEMP = DTEMP + DABS(DX(I)) 10 CONTINUE DASUM = DTEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,6 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 50 CONTINUE 60 DASUM = DTEMP RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*), DY(*), DA INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N .LE. 0) RETURN IF (DA .EQ. 0.0D0) RETURN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1, N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE C RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C CLEAN-UP LOOP C 20 CONTINUE M = MOD(N,4) IF (M .EQ. 0) GO TO 40 DO 30 I = 1, M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF (N .LT. 4) RETURN 40 CONTINUE MP1 = M + 1 DO 50 I = MP1, N, 4 DY(I) = DY(I) + DA*DX(I) DY(I+1) = DY(I+1) + DA*DX(I+1) DY(I+2) = DY(I+2) + DA*DX(I+2) DY(I+3) = DY(I+3) + DA*DX(I+3) 50 CONTINUE C RETURN C END C SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*), DY(*) INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N .LE. 0) RETURN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) GO TO 20 C C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1, N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE C RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 CONTINUE M = MOD(N,7) IF (M .EQ. 0) GO TO 40 DO 30 I = 1, M DY(I) = DX(I) 30 CONTINUE IF (N .LT. 7) RETURN 40 CONTINUE MP1 = M + 1 DO 50 I = MP1, N, 7 DY(I) = DX(I) DY(I+1) = DX(I+1) DY(I+2) = DX(I+2) DY(I+3) = DX(I+3) DY(I+4) = DX(I+4) DY(I+5) = DX(I+5) DY(I+6) = DX(I+6) 50 CONTINUE C RETURN C END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER NEXT,N,INCX,I,J,NN DOUBLE PRECISION DX(*), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DGEMV PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS * * Y := ALPHA*A*X + BETA*Y, OR Y := ALPHA*A'*X + BETA*Y, * * WHERE ALPHA AND BETA ARE SCALARS, X AND Y ARE VECTORS AND A IS AN * M BY N MATRIX. * * PARAMETERS * ========== * * TRANS - CHARACTER*1. * ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS * FOLLOWS: * * TRANS = 'N' OR 'N' Y := ALPHA*A*X + BETA*Y. * * TRANS = 'T' OR 'T' Y := ALPHA*A'*X + BETA*Y. * * TRANS = 'C' OR 'C' Y := ALPHA*A'*X + BETA*Y. * * UNCHANGED ON EXIT. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX A. * M MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX A. * N MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY A MUST * CONTAIN THE MATRIX OF COEFFICIENTS. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDA MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( 1 + ( N - 1 )*ABS( INCX ) ) WHEN TRANS = 'N' OR 'N' * AND AT LEAST * ( 1 + ( M - 1 )*ABS( INCX ) ) OTHERWISE. * BEFORE ENTRY, THE INCREMENTED ARRAY X MUST CONTAIN THE * VECTOR X. * UNCHANGED ON EXIT. * * INCX - INTEGER. * ON ENTRY, INCX SPECIFIES THE INCREMENT FOR THE ELEMENTS OF * X. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS * SUPPLIED AS ZERO THEN Y NEED NOT BE SET ON INPUT. * UNCHANGED ON EXIT. * * Y - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( 1 + ( M - 1 )*ABS( INCY ) ) WHEN TRANS = 'N' OR 'N' * AND AT LEAST * ( 1 + ( N - 1 )*ABS( INCY ) ) OTHERWISE. * BEFORE ENTRY WITH BETA NON-ZERO, THE INCREMENTED ARRAY Y * MUST CONTAIN THE VECTOR Y. ON EXIT, Y IS OVERWRITTEN BY THE * UPDATED VECTOR Y. * * INCY - INTEGER. * ON ENTRY, INCY SPECIFIES THE INCREMENT FOR THE ELEMENTS OF * Y. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * SET LENX AND LENY, THE LENGTHS OF THE VECTORS X AND Y, AND SET * UP THE START POINTS IN X AND Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * START THE OPERATIONS. IN THIS VERSION THE ELEMENTS OF A ARE * ACCESSED SEQUENTIALLY WITH ONE PASS THROUGH A. * * FIRST FORM Y := BETA*Y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * FORM Y := ALPHA*A*X + Y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * FORM Y := ALPHA*A'*X + Y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * END OF DGEMV . * END SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 3/93 TO RETURN IF INCX .LE. 0. C DOUBLE PRECISION DA, DX(*) INTEGER I, INCX, M, MP1, N, NINCX C IF (N .LE. 0 .OR. INCX .LE. 0) RETURN IF (INCX .EQ. 1) GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1, NINCX, INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 CONTINUE M = MOD(N,5) IF (M .EQ. 0) GO TO 40 DO 30 I = 1, M DX(I) = DA*DX(I) 30 CONTINUE IF (N .LT. 5) RETURN 40 CONTINUE MP1 = M + 1 DO 50 I = MP1, N, 5 DX(I) = DA*DX(I) DX(I+1) = DA*DX(I+1) DX(I+2) = DA*DX(I+2) DX(I+3) = DA*DX(I+3) DX(I+4) = DA*DX(I+4) 50 CONTINUE RETURN END SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) C C INTERCHANGES TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*), DY(*), DTEMP INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N .LE. 0) RETURN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1, N DTEMP = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 CONTINUE M = MOD(N,3) IF (M .EQ. 0) GO TO 40 DO 30 I = 1, M DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP 30 CONTINUE IF (N .LT. 3) RETURN 40 CONTINUE MP1 = M + 1 DO 50 I = MP1, N, 3 DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP DTEMP = DX(I+1) DX(I+1) = DY(I+1) DY(I+1) = DTEMP DTEMP = DX(I+2) DX(I+2) = DY(I+2) DY(I+2) = DTEMP 50 CONTINUE RETURN END INTEGER FUNCTION IDAMAX(N,DX,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 3/93 TO RETURN IF INCX .LE. 0. C DOUBLE PRECISION DX(*), DMAX INTEGER I, INCX, IX, N C IDAMAX = 0 IF (N .LT. 1 .OR. INCX .LE. 0) RETURN IDAMAX = 1 IF (N .EQ. 1) RETURN IF (INCX .EQ. 1) GO TO 30 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO 20 I = 2, N IF (DABS(DX(IX)) .LE. DMAX) GO TO 10 IDAMAX = I DMAX = DABS(DX(IX)) 10 CONTINUE IX = IX + INCX 20 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 30 CONTINUE DMAX = DABS(DX(1)) DO 40 I = 2, N IF (DABS(DX(I)) .LE. DMAX) GO TO 40 IDAMAX = I DMAX = DABS(DX(I)) 40 CONTINUE RETURN END LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK AUXILIARY ROUTINE (VERSION 1.1) -- * UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD., * COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY * FEBRUARY 29, 1992 * * .. SCALAR ARGUMENTS .. CHARACTER CA, CB * .. * * PURPOSE * ======= * * LSAME RETURNS .TRUE. IF CA IS THE SAME LETTER AS CB REGARDLESS OF * CASE. * * ARGUMENTS * ========= * * CA (INPUT) CHARACTER*1 * CB (INPUT) CHARACTER*1 * CA AND CB SPECIFY THE SINGLE CHARACTERS TO BE COMPARED. * * .. INTRINSIC FUNCTIONS .. INTRINSIC ICHAR * .. * .. LOCAL SCALARS .. INTEGER INTA, INTB, ZCODE * .. * .. EXECUTABLE STATEMENTS .. * * TEST IF THE CHARACTERS ARE EQUAL * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * NOW TEST FOR EQUIVALENCE IF BOTH CHARACTERS ARE ALPHABETIC. * ZCODE = ICHAR( 'Z' ) * * USE 'Z' RATHER THAN 'A' SO THAT ASCII CAN BE DETECTED ON PRIME * MACHINES, ON WHICH ICHAR RETURNS A VALUE WITH BIT 8 SET. * ICHAR('A') ON PRIME MACHINES RETURNS 193 WHICH IS THE SAME AS * ICHAR('A') ON AN EBCDIC MACHINE. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII IS ASSUMED - ZCODE IS THE ASCII CODE OF EITHER LOWER OR * UPPER CASE 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC IS ASSUMED - ZCODE IS THE EBCDIC CODE OF EITHER LOWER OR * UPPER CASE 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII IS ASSUMED, ON PRIME MACHINES - ZCODE IS THE ASCII CODE * PLUS 128 OF EITHER LOWER OR UPPER CASE 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * END OF LSAME * END SUBROUTINE XERBLA ( SRNAME, INFO ) * .. SCALAR ARGUMENTS .. INTEGER INFO CHARACTER*6 SRNAME * .. * * PURPOSE * ======= * * XERBLA IS AN ERROR HANDLER FOR THE LEVEL 2 BLAS ROUTINES. * * IT IS CALLED BY THE LEVEL 2 BLAS ROUTINES IF AN INPUT PARAMETER IS * INVALID. * * INSTALLERS SHOULD CONSIDER MODIFYING THE STOP STATEMENT IN ORDER TO * CALL SYSTEM-SPECIFIC EXCEPTION-HANDLING FACILITIES. * * PARAMETERS * ========== * * SRNAME - CHARACTER*6. * ON ENTRY, SRNAME SPECIFIES THE NAME OF THE ROUTINE WHICH * CALLED XERBLA. * * INFO - INTEGER. * ON ENTRY, INFO SPECIFIES THE POSITION OF THE INVALID * PARAMETER IN THE PARAMETER-LIST OF THE CALLING ROUTINE. * * * AUXILIARY ROUTINE FOR LEVEL 2 BLAS. * * WRITTEN ON 20-JULY-1986. * * .. EXECUTABLE STATEMENTS .. * WRITE (*,99999) SRNAME, INFO * STOP * 99999 FORMAT ( ' ** ON ENTRY TO ', A6, ' PARAMETER NUMBER ', I2, $ ' HAD AN ILLEGAL VALUE' ) * * END OF XERBLA. * END SHAR_EOF fi # end of overwriting check if test -f 'dpmeps.f' then echo shar: will not over-write existing file "'dpmeps.f'" else cat << \SHAR_EOF > 'dpmeps.f' DOUBLE PRECISION FUNCTION DPMEPS() C ********** C C SUBROUTINE DPMEPS C C THIS SUBROUTINE COMPUTES THE MACHINE PRECISION PARAMETER C DPMEPS AS THE SMALLEST FLOATING POINT NUMBER SUCH THAT C 1 + DPMEPS DIFFERS FROM 1. C C THIS SUBROUTINE IS BASED ON THE SUBROUTINE MACHAR DESCRIBED IN C C W. J. CODY, MACHAR: A SUBROUTINE TO DYNAMICALLY DETERMINE C MACHINE PARAMETERS, ACM TRANS. MATH. SOFTWARE 14 (1988), 303-311. C C THE FUNCTION STATEMENT IS C C DOUBLE PRECISION DPMEPS() C C MINPACK-2 PROJECT. FEBRUARY 1991. C ARGONNE NATIONAL LABORATORY AND UNIVERSITY OF MINNESOTA. C BRETT M. AVERICK. C C ********** INTEGER I, IBETA, IRND, IT, ITEMP, NEGEP DOUBLE PRECISION A, B, BETA, BETAIN, BETAH, TEMP, TEMPA, TEMP1 DOUBLE PRECISION ZERO, ONE, TWO DATA ZERO, ONE, TWO/0.0D0, 1.0D0, 2.0D0/ C DETERMINE IBETA, BETA ALA MALCOLM. A = ONE B = ONE 10 CONTINUE A = A + A TEMP = A + ONE TEMP1 = TEMP - A IF (TEMP1-ONE .EQ. ZERO) GO TO 10 20 CONTINUE B = B + B TEMP = A + B ITEMP = INT(TEMP-A) IF (ITEMP .EQ. 0) GO TO 20 IBETA = ITEMP BETA = DBLE(IBETA) C DETERMINE IT, IRND. IT = 0 B = ONE 30 CONTINUE IT = IT + 1 B = B*BETA TEMP = B + ONE TEMP1 = TEMP - B IF (TEMP1-ONE .EQ. ZERO) GO TO 30 IRND = 0 BETAH = BETA/TWO TEMP = A + BETAH IF (TEMP-A .NE. ZERO) IRND = 1 TEMPA = A + BETA TEMP = TEMPA + BETAH IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 C DETERMINE DPMEPS. NEGEP = IT + 3 BETAIN = ONE/BETA A = ONE DO 40 I = 1, NEGEP A = A*BETAIN 40 CONTINUE 50 CONTINUE TEMP = ONE + A IF (TEMP-ONE .NE. ZERO) GO TO 60 A = A*BETA GO TO 50 60 CONTINUE DPMEPS = A IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 70 A = (A*(ONE+A))/TWO TEMP = ONE + A IF (TEMP-ONE .NE. ZERO) DPMEPS = A 70 CONTINUE END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' c c TENSOLVE: A Software Package for Solving Systems of Nonlinear c Equations and Nonlinear Least Squares Problems Using c Tensor Methods. c c AUTHORS: Ali Bouaricha c Argonne National Laboratory c MCS Division c e-mail: bouarich@mcs.anl.gov c AND c Robert B. Schnabel c University of colorado at Boulder c Department of computer Science c e-mail: bobby@cs.colorado.edu c c DATE: Version of January, 1997 c c Purpose of Tensolve: c c TENSOLVE finds roots of systems of n nonlinear equations in n c unknowns, or minimizers of the sum of squares of m > n nonlinear c equations in n unknowns. It allows the user to choose between c a tensor method based on a quadratic model and a standard method c based on a linear model. Both models calculate an analytic or c finite-difference Jacobian matrix at each iteration. The tensor c method augments the linear model with a low-rank, second-order c term that is chosen so that the model is hardly more expensive c to form, store, or solve than the standard linear model. Either c a line search or trust-region step-selection strategy may be c selected. The tensor method is significantly more efficient c than the standard method in iterations, function evaluations, and c time. It is especially useful on problems where the Jacobian at c the solution is singular. c The software can be called with an interface where the user c supplies only the function, number of nonlinear equations, number c of variables, and starting point; default choices are made for c all the other input parameters. An alternative interface allows c the user to specify any input parameters that are different from c the defaults. c c List of subroutine and function names called by TENSOLVE: c c TS2DTR,TSBSLV,TSCHKI,TSCHKJ,TSCPMU,TSCPSS,TSDLOD,TSD1SV,TSDFCN,TSDFLT, c TSDUMJ,TSFAFA,TSFDFJ,TSFRMT,TSFSCL,TSFSLV,TSJMUV,TSJQTP,TSLMIN,TSLMSP, c TSLSCH,TSMAFA,TSMDLS,TSMFDA,TSMFDV,TSMGSA,TSMSDA,TSMSDV,TSMSLV,TSNECI, c TSNESI,TSNESV,TSNSTP,TSPRMV,TSRSLT,TSQ1P1,TSQFCN,TSQLFC,TSQMLV,TSQMTS, c TSQMUV,TSQRFC,TSRSID,TSSCLF,TSSCLJ,TSSCLX,TSSLCT,TSSMIN,TSSMRD,TSSQP1, c TSSSTP,TSSTMX,TSTRUD,TSUDQV,TSUNSF,TSUNSX,TSUPSF,TSUTMD. c c Packages called by TENSOLVE: c c UNCMIN (R. B. Schnabel, J. E. Koontz, and B. E. Weiss, c "A Modular System of Algorithms of Unconstrained Minimization", c Acm Trans. Math. Softw., 11 (1985), 419-440). c c BLAS called by TENSOLVE: c c LEVEL 1 BLAS: DASUM,DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP,IDAMAX c LEVEL 2 BLAS: DGEMV c c Parameters and Default Values for the interfaces TSNECI and TSNESI: c ================================================================== c c Following each variable name in the list below appears a one- or c two-headed arrow symbol of the form ->, <-, and <-->. c These symbols signify that the variable is for input, output, and c input-output, respectively. c The symbol EPSM in some parts of this section designates the machine c precision. c MAXM->: A positive integer specifying the row dimension of the work c array WRKNEM in the user's calling program. It must satisfy c MAXM >= M+N+2. The provision of MAXM, MAXN, and MAXP allows c the user the flexibility of solving several problems with different c values of M and N one after the other, with the same work arrays. c MAXN->: A positive integer specifying the row dimension of the work c array WRKNEN in the user's calling program. It must satisfy c MAXN >= N+2. c MAXP->: A positive integer specifying the row dimension of the work c array WRKUNC in the user's calling program. It must satisfy c MAXP >= NINT(sqrt(N)), where NINT is a function that rounds to the c nearest integer. c X0->: An array of length N that contains an initial estimate c of the solution x*. c M->: A positive integer specifying the number of nonlinear equations. c N->: A positive integer specifying the number of variables in the c problem. c TYPX->: An array of length N in which the typical size of the C components of X is specified. The typical component sizes should be c positive real scalars. If a negative value is specified, its absolute c value will be used. If 0.0 is specified, 1.0 will be used. This c vector is used to determine the scaling matrix, Dx. Although the c package may work reasonably well in a large number of instances without c scaling, it may fail when the components of x* are of radically c different magnitude and scaling is not invoked. If the sizes c of the parameters are known to differ by many orders of magnitude, then c the scale vector TYPX should definitely be used. For example, if c it is anticipated that the range of values for the iterates xk would be c x1 in [-10e+10,10e+10] c x2 in [-10e+2,10e+4] c x3 in [-6*10e-6,9*10e-6] c then an appropriate choice would be TYPX = (1.0e+10,1.0e+3,7.0e-6). c Module TSDFLT returns TYPX = (1.0,...,1.0). c TYPF->: An array of length M in which the typical size of the components c of F is specified. The typical component sizes should be positive real c scalars. If a negative value is specified, its absolute value will be c used. If 0.0 is specified, 1.0 will be used. This vector is used to c determine the scaling matrix DF. TYPF should be chosen so that all c the components of DF(x) have similar typical magnitudes at points not c too near a root, and should be chosen in conjunction with FTOL. It is c important to supply values of TYPF when the magnitudes of the components c of F(x) are expected to be very different. If the magnitudes of the c components of F(x) are similar, the choice DF = I suffices. Module c TSDFLT returns TYPF = (1.0,...,1.0). c ITNLIM->: Positive integer specifying the maximum number of iterations c to be performed before the program is terminated. Module TSDFLT returns c ITNLIM = 150. If the user specifies ITNLIM <= 0, the module TSCHKI will c supply the value 150. c JACFLG->: Integer designating whether or not an analytic Jacobian has c been supplied by the user. c JACFLG = 0 : No analytic Jacobian supplied. The Jacobian is obtained c by finite differences. c JACFLG = 1 : Analytic Jacobian supplied. c The module TSDFLT returns the value 0. If the user specifies an illegal c value, the module TSCHKI will supply the value 0. c GRADTL->: Positive scalar giving the tolerance at which the scaled c gradient of f(x) = 0.5*F(x)-trans*F(x) is considered close enough to c zero to terminate the algorithm. The scaled gradient is a measure of c the relative change in F in each direction xj divided by the relative c change in xj. The module TSDFLT returns the value EPSM**(1/3). If the c user specifies a negative value, the module TSCHKI will supply c the value EPSM**(1/3). c STEPTL->: A positive scalar providing the minimum allowable relative c step length. STEPTL should be at least as small as 10**(-d), where d c is the number of accurate digits the user desires in the solution x*. c The program may terminate prematurely if STEPTL is too large. Module c TSDFLT returns the value EPSM**(2/3). If the user specifies a negative c value, the module TSCHKI will supply the value EPSM**(2/3). c FTOL->: A positive scalar giving the tolerance at which the scaled c function DF*F(x) is considered close enough to zero to terminate the c algorithm. The program is halted if ||DF*F(x)|| (in the infinity norm) c is <= FTOL. This is the primary stopping condition for nonlinear c equations; the values of TYPF and FTOL should be chosen so that this c test reflects the user's idea of what constitutes a solution to the c problem. The module TSDFLT returns the value EPSM**(2/3). If the c user specifies a negative value, the module TSCHKI will supply the c value EPSM**(2/3). c METHOD->: An integer designating which method to use. c METHOD = 0 : Newton or Gauss-Newton algorithm is used. c METHOD = 1 : Tensor algorithm is used. c Module TSDFLT returns value 1. If the user specifies an illegal value, c module TSCHKI will reset METHOD to 1. c GLOBAL->: An integer designating which global strategy to use. c GLOBAL = 0 : Line search is used. c GLOBAL = 1 : Two-dimensional trust region is used. c Module TSDFLT returns value of 0. If the user specifies an illegal c value, module TSCHKI will reset GLOBAL to 0. c STEPMX->: A positive scalar providing the maximum allowable scaled step c length ||Dx*(x+ - xc)||2, where Dx = diag(1/TYPX_j). STEPMX is used to c prevent steps that would cause the nonlinear equations problem to c overflow, and to prevent the algorithm from leaving the area of c interest in parameter space. STEPMX should be chosen small enough c to prevent these occurrences but should be larger than any anticipated c "reasonable" step. Module TSDFLT returns the value STEPMX = 10e+3. c If the user specifies a nonpositive value, module TSCHKI sets STEPMX c to 10e+3. c DLT->: A positive scalar giving the initial trust region radius. When c the line search strategy is used, this parameter is ignored. For the c trust region algorithm, if DLT is supplied, its value should reflect c what the user considers a maximum reasonable scaled step length at c the first iteration. If DLT = -1.0, the routine uses the length of c the Cauchy step at the initial iterate instead. The module TSDFLT c returns the value -1.0. If the user specifies a nonpositive value, c module TSCHKI sets DLT = -1.0. c IPR->: The unit on which the package outputs information. TSDFLT returns c the value 6. c WRKUNC->: Workspace used by UNCMIN. The user must declare this c array to have dimensions MAXP*LUNC in the calling routine. c LUNC->: A positive integer specifying the column dimension of the work c array WRKUNC in the user's calling program. It must satisfy c LUNC >= 2*NINT(sqrt(N))+4. c WRKNEM->: Workspace used to store the Jacobian matrix, the function c values matrix FV, the tensor matrix ANLS, and working vectors. The c user must declare this array to have dimensions MAXM*LNEM in the c calling routine. c LNEM->: A positive integer specifying the column dimension of the work c array WRKNEM in the user's calling program. It must satisfy c LNEM >= N+2*NINT(sqrt(N))+11. c WRKNEN->: Workspace used to store the matrix S of previous c directions, the matrix SHAT of linearly independent directions, and c working vectors. The user must declare this array to have dimensions c MAXN*LNEN in the calling routine. c LNEN->: A positive integer specifying the column dimension of the work c array WRKNEN in the user's calling program. It must satisfy c LNEN >= 2*NINT(sqrt(N))+9. c IWRKN->: Workspace used to store the integer working vectors. The user c must declare this array to have dimensions at least MAXN*LIN in the c calling routine. c LIN->: A positive integer specifying the column dimension of the work c array IWRKN in the user's calling program. It must satisfy c LIN >= 3. c FVEC->: The name of a user-supplied subroutine that evaluates the c function F at an arbitrary vector X. The subroutine must c be declared EXTERNAL in the user's calling program and must conform c to the usage c CALL FVEC(X, F, M, N), c where X is a vector of length N and F is a vector of length M. The c subroutine must not alter the values of X. c JAC->: The name of a user-supplied subroutine that evaluates the first c derivative (Jacobian) of the function F(x). The subroutine must be c declared EXTERNAL in the user's program and must conform to the usage c CALL JAC(X, AJA, MAXM, M, N) c where X is a vector of length N and the 2-dimensional array AJA of row c dimension MAXM and column dimension N is the analytic Jacobian of F at c X. When using the interface TSNECI, if no analytic Jacobian is supplied c (JACFLG = 0), the user must use the dummy name TSDUMJ as the value of c this parameter. c MSG<-->: An integer variable that the user may set on input to inhibit c certain automatic checks or to override certain default characteristics c of the package. (For the short call it should be set to 0.) There are c four "message" features that can be used individually or in combination c as discussed below. c MSG = 0 : Values of input parameters, final results, and termination code c are printed. c MSG = 2 : Do not check user's analytic Jacobian routine against its c finite difference estimate. This may be necessary if the user knows the c Jacobian is properly coded, but the program aborts because the comparative c tolerance is too tight. Do not use MSG = 2 if the analytic acobian is c not supplied. c MSG = 8 : Suppress printing of the input state, the final results, and c the stopping condition. c MSG = 16 : Print the intermediate results; that is, the input state, c each iteration including the current iterate xk, 0.5*||DF*F(xk)||2**2, c and grad(f(x)) = J(x)-trans*DF**2 F(x), and the final results including c the stopping conditions. c The user may specify a combination of features by setting MSG to c the sum of the individual components. The module TSDFLT returns a value c of 0. On exit, if the program has terminated because of erroneous c input, MSG contains an error code indicating the reason. c MSG = 0 : No error. c MSG = -1 : Illegal dimension, M <= 0. c MSG = -2 : Illegal dimension, N <= 0. c MSG = -3 : Illegal dimension, MAXM < M+N+2. c MSG = -4 : Illegal dimension, MAXN < N+2. c MSG = -5 : Illegal dimension, MAXP < NINT(sqrt(N)). c MSG = -6 : Illegal dimension, LUNC < 2*NINT(sqrt(N))+4. c MSG = -7 : Illegal dimension, LNEM < N+2*NINT(sqrt(N))+11. c MSG = -8 : Illegal dimension, LNEN < 2*NINT(sqrt(N))+9. c MSG = -9 : Illegal dimension, LIN < 3. c MSG = -10 : Program asked to override check of analytic Jacobian c against finite difference estimate, but routine JAC not c supplied (incompatible input). c MSG = -11 : Probable coding error in the user's analytic Jacobian c routine JAC. Analytic and finite difference Jacobian do not agree c within the assigned tolerance. c XP<-: An array of length N containing the best approximation c to the solution x*. (If the algorithm has not converged, the final c iterate is returned). c FP<-: An array of length M containing the function value F(XP). c GP<-: An array of length N containing the gradient of the c function 0.5*||F(x)||2**2 at XP. c TERMCD<-: An integer specifying the reason for termination. c TERMCD = 0 : No termination criterion satisfied (occurs if package c terminates because of illegal input). c TERMCD = 1 : function tolerance reached. The current iteration is c probably a solution. c TERMCD = 2 : gradient tolerance reached. For nonlinear least c squares, the current iteration is probably a solution; for nonlinear c equations, it could be a solution or a local minimizer. c TERMCD = 3 : Successive iterates within step tolerance. The c current iterate may be a solution, or the algorithm is making very slow c progress and is not near a solution. c TERMCD = 4 : Last global step failed to locate a point lower c than XP. It is likely that either XP is an approximate solution c of the problem or STEPTL is too large. c TERMCD = 5 : Iteration limit exceeded. c TERMCD = 6 : Five consecutive steps of length STEPMX have been taken. c c =========================================================================== c Begin TENSOLVE c =========================================================================== SUBROUTINE TS2DTR(AJA,SHAT,ANLS,DT,G,GBAR,XC,METHOD,NWTAKE, + STEPMX,STEPTL,EPSM,MXTAKE,DLT,FQ,MAXM,MAXN, + M,N,P,CURPOS,PIVOT,PBAR,ITN,IERR,FLAG, + DXN,DFN,FVEC,D,XPLSP,ADT,AG,TEMP,VN,VNP,VNS, + WRK1,CONST1,CONST2,FNORM,XPLS,FP,FPLS,RETCD) INTEGER MAXM,MAXN,M,N,P,ITN,METHOD,IERR,FLAG DOUBLE PRECISION STEPMX,STEPTL,EPSM,DLT,FPLS INTEGER CURPOS(N),PIVOT(N),PBAR(N),RETCD DOUBLE PRECISION DT(N),G(N),GBAR(N),XC(N) DOUBLE PRECISION XPLS(N),FP(M),XPLSP(N),AJA(MAXM,N),D(M) DOUBLE PRECISION TEMP(M),SHAT(MAXN,P),ANLS(MAXM,P),VNS(M) DOUBLE PRECISION VN(M),VNP(M),FQ(M),DXN(N),DFN(M) DOUBLE PRECISION ADT(N),AG(N),WRK1(M),CONST1(P),CONST2(P) LOGICAL NWTAKE,MXTAKE C********************************************************************** C THIS ROUTINE FINDS A NEXT ITERATE BY A 2-DIMENSIONAL TRUST REGION. C********************************************************************** C C C INPUT PARAMETERS : C ----------------- C C AJA : JACOBIAN AT THE CURRENT ITERATE C SHAT : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS C AFTER A QL FACTORIZATION C ANLS : TENSOR TERM MATRIX C DT : CURRENT STEP C G : GRADIENT AT CURRENT ITERATE C GBAR : STEEPEST DESCENT DIRECTION (= -G) C XC : CURRENT ITERATE C METHOD : METHOD TO USE C = 0 : STANDARD METHOD USED C = 1 : TENSOR METHOD USED C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C STEPMX : MAXIMUM STEP ALLOWED C STEPTL : STEP TOLERANCE C EPSM : MACHINE PRECISION C MXTAKE : BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY C ORTHOGONAL MATRICES C MAXM : LEADING DIMENSION OF AJA AND ANLS C MAXN : LEADING DIMENSION OF SHAT C M,N : DIMENSIONS OF PROBLEM C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS C CURPOS : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE C JACOBIAN FROM COLUMN 1 TO N-P) C PIVOT : PIVOT VECTOR ( USED DURING THE FACTORIZATION OF THE C JACOBIAN FROM COLUMN N-P+1 TO N) C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE C JACOBIAN IF IT IS SINGULAR C FNORM : 0.5 * || FC ||**2 C ITN : ITERATION NUMBER C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE: C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C FLAG : RETURN CODE WITH THE FOLLOWING MEANINGS : C FLAG = 0 : NO SINGULARITY DETECTED DURING C FACTORIZATION OF THE JACOBIAN FROM C COLUMN 1 TO N C FLAG = 1 : SINGULARITY DETECTED DURING FACTORIZATION C OF THE JACOBIAN FROM COLUMN 1 TO N-P C FLAG = 2 : SINGULARITY DETECTED DURING FACTORIZATION C OF THE JACOBIAN FROM COLUMN N-P+1 TO N C DXN : DIAGONAL SCALING MATRIX FOR X C DFN : DIAGONAL SCALING MATRIX FOR F C FVEC : SUBROUTINE TO EVALUATE THE USER'S FUNCTION C D,XPLSP,ADT,AG,TEMP,VN,VNP,VNS : WORKING VECTORS C WRK1,CONST1,CONST2,X: WORKING VECTORS C C INPUT-OUTPUT PARAMETERS : C ------------------------ C C DLT : INITIAL TRUST RADIUS (= -1.0D0) IF IT IS NOT SUPPLIED C BY THE USER ON ENTRY AND CURRENT TRUST RADIUS ON EXIT C C OUTPUT PARAMETERS : C ------------------- C C XPLS : NEXT ITERATE C FP : FUNCTION VALUE AT NEXT ITERATE C FPLS : 0.5 * || FP ||**2 C RETCD : RETURN CODE FROM SUBROUTINE (SEE SUBROUTINE TSTRUD C FOR MEANING ) C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DAXPY,DCOPY,DDOT,DNRM2,DSCAL C TENSOLVE ... TSPRMV,TSUTMD,TSJMUV,TSUDQV,TSSMIN,TSRSID, C TENSOLVE ... TSTRUD C C*********************************************************************** INTEGER I DOUBLE PRECISION FNORM,RES,ALPH,SUM,RESG,OPTIM DOUBLE PRECISION SCRES,FPLSP,RRES,RRESG DOUBLE PRECISION DNRM2,DDOT DOUBLE PRECISION ZERO,ONE LOGICAL DTAKEN INTRINSIC SQRT EXTERNAL FVEC DATA ZERO,ONE/0.0D0,1.0D0/ DTAKEN = .FALSE. RETCD = 4 IF(DLT.EQ.-ONE) THEN c set DLT to length of Cauchy step ALPH = DNRM2(N,G,1) ALPH = ALPH**2 CALL TSPRMV(VN,G,PIVOT,N,1) IF(IERR.EQ.0) THEN CALL TSUTMD(AJA,VN,MAXM,M,N,VNP) ELSE CALL TSPRMV(VNS,VN,PBAR,N,1) CALL TSUTMD(AJA,VNS,MAXM,M+N,N,VNP) ENDIF DLT = ALPH*SQRT(ALPH)/DNRM2(N,VNP,1)**2 IF(DLT.GT.STEPMX) THEN DLT = STEPMX ENDIF ENDIF c form an orthonormal basis for the two-dimensional subspace CALL DCOPY(N,G,1,GBAR,1) CALL DSCAL(N,-ONE,GBAR,1) RES = DNRM2(N,DT,1) SUM = -DDOT(N,GBAR,1,DT,1)/RES**2 CALL DAXPY(N,SUM,DT,1,GBAR,1) RESG = DNRM2(N,GBAR,1) IF(RESG.GT.ZERO) THEN RRESG = ONE/RESG CALL DSCAL(N,RRESG,GBAR,1) ENDIF RRES = ONE/RES CALL DSCAL(N,RRES,DT,1) c compute Jacobian times DT CALL TSJMUV(ITN,METHOD,DT,CURPOS,PIVOT,PBAR,AJA,SHAT, + FLAG,IERR,MAXM,MAXN,M,N,P,D,TEMP,VN,ADT) c compute Jacobian times GBAR CALL TSJMUV(ITN,METHOD,GBAR,CURPOS,PIVOT,PBAR,AJA,SHAT, + FLAG,IERR,MAXM,MAXN,M,N,P,D,TEMP,VNP,AG) IF(.NOT. NWTAKE) THEN c compute SHAT times VN CALL TSUDQV(SHAT,VN,MAXN,N,P,CONST1) c compute SHAT times VNP CALL TSUDQV(SHAT,VNP,MAXN,N,P,CONST2) ENDIF 70 CONTINUE c normalize DT IF(RES.LE.DLT) THEN DTAKEN = .TRUE. DO 80 I = 1,N D(I) = DT(I)*RES 80 CONTINUE DLT = RES ELSE c find the global minimizer of one-variable function in the c interval (-dlt, dlt) CALL TSSMIN(ANLS,FQ,ADT,AG,CONST1,CONST2,DLT,MAXM,M,N, + P,NWTAKE,IERR,EPSM,VN,VNP,VNS,OPTIM) c compute the global step DO 90 I = 1,N D(I) = OPTIM*DT(I)+SQRT(DLT**2-OPTIM**2)*GBAR(I) 90 CONTINUE ENDIF c compute the tensor model residual CALL TSRSID(ITN,METHOD,FQ,D,CURPOS,PIVOT,PBAR,AJA,ANLS, + SHAT,FLAG,NWTAKE,IERR,MAXM,MAXN,M,N,P,WRK1,VN, + VNP,VNS,SCRES) c check whether the global step is acceptable CALL TSTRUD(M,N,XC,FNORM,G,D,DTAKEN,STEPMX,STEPTL,DLT, + MXTAKE,DXN,DFN,FVEC,SCRES,RETCD,XPLSP,FPLSP, + TEMP,XPLS,FP,FPLS) IF(RETCD.GE.2) GO TO 70 END SUBROUTINE TSBSLV(R,NR,M,N,B,Y) INTEGER NR,M,N DOUBLE PRECISION R(NR,N),B(N),Y(N) C********************************************************************* C THIS ROUTINE DOES A BACKWARD SOLVE. C********************************************************************* C C INPUT PARAMETERS : C ----------------- C C R : UPPER TRIANGULAR MATRIX OBTAINED FROM A QR FACTORIZATION C OF AN M BY N MATRIX A. DIAG(R) IS STORED IN ROW M+2. THIS C IS THE STORAGE SCHEME USED IN STEWART, G. W., III(1973) C "INTRODUCTION TO MATRIX COMPUTATION", ACADEMIC PRESS, C NEW YORK C NR : LEADING DIMENSION OF MATRIX A C M : ROW DIMENSION OF MATRIX A C N : COLUMN DIMENSION OF MATRIX A C B : RIGHT HAND SIDE C C OUTPUT PARAMETERS : C ----------------- C C Y : VECTOR SOLUTION ON EXIT C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DAXPY C TENSOLVE ... TSDLOD C C********************************************************************* INTEGER J DOUBLE PRECISION ZERO DATA ZERO/0.0D0/ c solve R Y = B Y(N) = B(N) / R(M+2,N) IF(N .GT. 2) THEN CALL TSDLOD(N-1,ZERO,Y,1) DO 20 J = N-1,2,-1 CALL DAXPY(J,Y(J+1),R(1,J+1),1,Y,1) Y(J) = (B(J)-Y(J))/R(M+2,J) 20 CONTINUE Y(1) = Y(1) + R(1,2) * Y(2) Y(1) = (B(1) - Y(1)) / R(M+2,1) ELSEIF(N .EQ. 2) THEN Y(1) = (B(1) - (R(1,2) * Y(2))) / R(M+2,1) ENDIF RETURN END SUBROUTINE TSCHKI(MAXM,MAXN,MAXP,M,N,LUNC,LNEM,LNEN,LIN,GRADTL, + STEPTL,FTOL,ITNLIM,JACFLG,METHOD,GLOBAL, + STEPMX,DLT,EPSM,MSG,TYPX,TYPF,DXN,DFN, + SQRN,TERMCD,IPR) INTEGER MAXM,MAXN,MAXP,M,N,LUNC,LNEM,LNEN,LIN,IPR,MSG,JACFLG INTEGER METHOD,GLOBAL,ITNLIM,SQRN,TERMCD DOUBLE PRECISION GRADTL,STEPTL,FTOL,STEPMX,DLT,EPSM DOUBLE PRECISION TYPX(N),TYPF(M),DXN(N),DFN(M) C********************************************************************* C THIS ROUTINE CHECKS INPUT FOR REASONABLENESS. C********************************************************************* C C INPUT PARAMETERS : C ----------------- C C MAXM : LEADING DIMENSION OF WORKSPACE WRKNEM C (SEE TOP OF THIS FILE ) C MAXN : LEADING DIMENSION OF WORKSPACE WRKNEN C (SEE TOP OF THIS FILE ) C MAXP : LEADING DIMENSION OF WORKSPACE WRKUNC C (SEE TOP OF THIS FILE ) C M,N : DIMENSIONS OF PROBLEM C IPR : DEVICE TO WHICH TO SEND OUTPUT C C INPUT-OUTPUT PARAMETERS : C ------------------------ C C GRADTL : TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE ALGORITHM C STEPTL : TOLERANCE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM C FTOL : TOLERANCE AT WHICH FUNCTION VALUE CONSIDERED C CLOSE ENOUGH TO ZERO C ITNLIM : MAXIMUM NUMBER OF ALLOWABLE ITERATIONS C STEPMX : MAXIMUM STEP ALLOWED IN TRUST REGION C DLT : TRUST RADIUS C JACFLG : JACOBIAN FLAG WITH THE FOLLOWING MEANINGS : C JACFLG = 1 : ANALYTIC JACOBIAN SUPPLIED C JACFLG = 0 : ANALYTIC JACOBIAN NOT SUPPLIED C METHOD : METHOD TO USE C METHOD = 0 : STANDARD METHOD IS USED C METHOD = 1 : TENSOR METHOD IS USED C GLOBAL : GLOBAL STRATEGY USED C GLOBAL = 0 : LINE SEARCH USED C GLOBAL = 1 : 2-DIMENSIONAL TRUST REGION USED C TYPX : TYPICAL SIZE FOR EACH COMPONENT OF X C TYPF : TYPICAL SIZE FOR EACH COMPONENT OF F C MSG : MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT C C OUTPUT PARAMETERS : C ------------------- C C TERMCD: TERMINATION CODE C DXN : DIAGONAL SCALING MATRIX FOR X C DFN : DIAGONAL SCALING MATRIX FOR F C SQRN : MAXIMUM COLUMN DIMENSION OF S AND FV C C SUBPROGRAMS CALLED: C C UNCMIN ... DPMEPS C C********************************************************************* INTEGER I,LEN DOUBLE PRECISION DPMEPS,ZERO,ONE,TWO,THREE,THOUS,TEMP INTRINSIC MOD,NINT,SQRT DATA ZERO,ONE,TWO,THREE,THOUS/0.0D0,1.0D0,2.0D0,3.0D0,1000.0D0/ c check that parameters only take on acceptable values c if not, set them to default values c set TERMCD to zero in case we abort prematuraly TERMCD = 0 c compute machine precision EPSM = DPMEPS() c check dimensions of the problem IF(M.LE.0) THEN WRITE(IPR,601) M MSG = -1 RETURN ENDIF IF(N.LE.0) THEN WRITE(IPR,602) N MSG = -2 RETURN ENDIF c check leading dimensions of the problem LEN = M+N+2 IF(MAXM .LT. LEN) THEN WRITE(IPR,603) MAXM,LEN MSG = -3 RETURN ENDIF LEN = N+2 IF(MAXN .LT. LEN) THEN WRITE(IPR,604) MAXN,LEN MSG = -4 RETURN ENDIF TEMP = SQRT(DBLE(N)) SQRN = NINT(TEMP) IF(MAXP .LT. SQRN) THEN WRITE(IPR,605) MAXP,SQRN MSG = -5 RETURN ENDIF c check column dimensions of workspace arrays LEN = 2*SQRN+4 IF(LUNC.LT.LEN) THEN WRITE(IPR,606) LUNC,LEN MSG = -6 RETURN ENDIF LEN = N+2*SQRN+11 IF(LNEM.LT.LEN) THEN WRITE(IPR,607) LNEM,LEN MSG = -7 RETURN ENDIF LEN = 2*SQRN+9 IF(LNEN.LT.LEN) THEN WRITE(IPR,608) LNEN,LEN MSG = -8 RETURN ENDIF IF(LIN.LT.3) THEN WRITE(IPR,609) LIN MSG = -9 RETURN ENDIF c check JACFLG, METHOD, and GLOBAL IF(JACFLG.NE.1) JACFLG = 0 IF(METHOD.NE.0 .AND. METHOD.NE.1) METHOD = 1 IF(GLOBAL.NE.0 .AND. GLOBAL.NE.1) GLOBAL = 0 IF(MOD(MSG/2,2).EQ.1 .AND. JACFLG.EQ.0) THEN WRITE(IPR,610) MSG,JACFLG MSG = -10 RETURN ENDIF c check scale matrices DO 10 I = 1,N IF(TYPX(I).EQ.ZERO) TYPX(I) = ONE IF(TYPX(I).LT.ZERO) TYPX(I) = -TYPX(I) DXN(I) = ONE/TYPX(I) 10 CONTINUE DO 20 I = 1,M IF(TYPF(I).EQ.ZERO) TYPF(I) = ONE IF(TYPF(I).LT.ZERO) TYPF(I) = -TYPF(I) DFN(I) = ONE/TYPF(I) 20 CONTINUE c check gradient, step, and function tolerances TEMP = ONE/THREE IF(GRADTL.LT.ZERO) THEN GRADTL = EPSM**TEMP ENDIF IF(STEPTL.LT.ZERO) THEN STEPTL = EPSM**(TWO*TEMP) ENDIF IF(FTOL.LT.ZERO) THEN FTOL = EPSM**(TWO*TEMP) ENDIF c check iteration limit IF(ITNLIM.LE.0) THEN ITNLIM = 150 ENDIF c check STEPMX and DLT IF(STEPMX.LT.ZERO) STEPMX = THOUS IF(DLT.LE.ZERO) THEN DLT = -ONE IF(DLT.GT.STEPMX) DLT = STEPMX ENDIF 601 FORMAT(' TSCHKI ILLEGAL DIMENSION M =',I5) 602 FORMAT(' TSCHKI ILLEGAL DIMENSION N =',I5) 603 FORMAT(' TSCHKI ILLEGAL DIMENSION MAXM =',I5,' < M+N+2 =',I5) 604 FORMAT(' TSCHKI ILLEGAL DIMENSION MAXN =',I5,' < N+2 =',I5) 605 FORMAT(' TSCHKI ILLEGAL DIMENSION MAXP =',I5,' <', + ' NINT(SQRT (N)) =',I5) 606 FORMAT(' TSCHKI ILLEGAL DIMENSION LUNC =',I5,' <', + ' 2*NINT(SQRT (N))+4 =',I5) 607 FORMAT(' TSCHKI ILLEGAL DIMENSION LNEM =',I5,' <', + ' N+2*NINT(SQRT (N))+11 =',I5) 608 FORMAT(' TSCHKI ILLEGAL DIMENSION LNEN =',I5,' <', + ' 2*NINT(SQRT (N))+9 =',I5) 609 FORMAT(' TSCHKI ILLEGAL DIMENSION LIN =',I5,' < 3') 610 FORMAT(' TSCHKI USER REQUESTS THAT ANALYTIC JACOBIAN BE', + ' ACCEPTED AS PROPERLY CODED (MSG =',I5,')'/ + ' TSCHKI BUT ANALYTIC JACOBIAN NOT SUPPLIED', + ' (JACFLG =',I5,')') END SUBROUTINE TSCHKJ(AJANAL,XC,FC,NR,M,N,EPSM,DFN,DXN, + TYPX,IPR,FHAT,WRK1,FVEC,MSG) INTEGER NR,M,N,IPR,MSG DOUBLE PRECISION AJANAL(NR,N),XC(N),FC(M) DOUBLE PRECISION EPSM,DFN(M),DXN(N),TYPX(N) DOUBLE PRECISION FHAT(M),WRK1(M) EXTERNAL FVEC C********************************************************************* C THIS ROUTINE CHECKS THE ANALYTIC JACOBIAN AGAINST ITS FINITE C DIFFERENCE APPROXIMATION. C********************************************************************* C C INPUT PARAMETERS : C ----------------- C C AJANAL : ANALYTIC JACOBIAN AT XC C XC : CURRENT ITERATE C FC : FUNCTION VALUE AT XC C NR : LEADING DIMENSION OF AJANAL C M,N : DIMENSIONS OF PROBLEM C EPSM : MACHINE PRECISION C DFN : DIAGONAL SCALING MATRIX FOR F C DXN : DIAGONAL SCALING MATRIX FOR X C TYPX : TYPICAL SIZE FOR EACH COMPONENT OF X C IPR : DEVICE TO WHICH TO SEND OUTPUT C FHAT,WRK1 : WORKSPACE C FVEC : SUBROUTINE TO EVALUATE THE USER'S FUNCTION C C INPUT-OUTPUT PARAMETERS : C ------------------------ C C MSG : MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... IDAMAX C TENSOLVE ... TSUNSX,TSUNSF,TSSCLX,TSSCLF C USER ... FVEC C C********************************************************************* INTEGER I,J DOUBLE PRECISION NDIGIT,RNOISE,SQRNS,STEPSZ,XTMPJ,DINF,RSTPSZ DOUBLE PRECISION TOL,QUART,ONE,TEN INTEGER IDAMAX INTRINSIC ABS,MAX,SQRT DATA QUART,ONE,TEN/0.250D0,1.0D0,10.0D0/ c unscale XC and FC CALL TSUNSX(XC,DXN,N) CALL TSUNSF(FC,DFN,M) c compute the finite difference Jacobian and check it against c the analytic one NDIGIT = -LOG10(EPSM) RNOISE = MAX(TEN**(-NDIGIT),EPSM) SQRNS = SQRT(RNOISE) TOL = EPSM**QUART DO 40 J = 1,N STEPSZ = SQRNS*MAX(ABS(XC(J)),ONE) XTMPJ = XC(J) XC(J) = XTMPJ+STEPSZ CALL FVEC(XC,FHAT,M,N) XC(J) = XTMPJ RSTPSZ = ONE/STEPSZ DO 10 I = 1,M WRK1(I) = (FHAT(I)-FC(I))*RSTPSZ 10 CONTINUE DO 20 I = 1,M WRK1(I) = WRK1(I)*DFN(I)*TYPX(J) 20 CONTINUE DINF = ABS(WRK1(IDAMAX(M,WRK1,1))) DO 30 I = 1,M IF(ABS(AJANAL(I,J)-WRK1(I)).GT.TOL*DINF) THEN WRITE(IPR,50) MSG = -11 RETURN ENDIF 30 CONTINUE 40 CONTINUE c scale back XC and FC CALL TSSCLX(XC,DXN,N) CALL TSSCLF(FC,DFN,M) 50 FORMAT(/,' TSCHKJ PROBABLE ERROR IN CODING OF ANALYTIC', + ' JACOBIAN') RETURN END SUBROUTINE TSCPMU(R,NR,N,EPSM,MU) INTEGER NR,N DOUBLE PRECISION R(NR,N),EPSM,MU C********************************************************************* C THIS ROUTINE COMPUTES A SMALL PERTURBATION MU. MU IS USED IN THE C COMPUTATION OF THE LEVENBERG-MARQUARDT STEP. C********************************************************************* C C INPUT PARAMETERS : C ----------------- C C R : UPPER TRIANGULAR MATRIX C NR : LEADING DIMENSION OF R C N : COLUMN DIMENSION OF R C EPSM : MACHINE PRECISION C C OUTPUT PARAMETERS : C ------------------ C C MU = SQRT(L1 NORM OF R * INFINITY NORM OF R * N * EPSM * 100) C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DASUM C C********************************************************************* INTEGER I,J DOUBLE PRECISION AIFNRM,SUM,AL1NRM,ZERO,HUND DOUBLE PRECISION DASUM INTRINSIC ABS,MAX,SQRT DATA ZERO,HUND/0.0D0,100.0D0/ c compute the infinity norm of R AIFNRM = ZERO DO 20 I = 1,N SUM = ZERO DO 10 J = I,N SUM = SUM+ABS(R(I,J)) 10 CONTINUE AIFNRM = MAX(AIFNRM,SUM) 20 CONTINUE c compute the l1 norm of R AL1NRM = ZERO DO 40 J = 1,N SUM = DASUM(J,R(1,J),1) AL1NRM = MAX(AL1NRM,SUM) 40 CONTINUE c compute MU MU = SQRT(AIFNRM*AL1NRM*N*EPSM*HUND) RETURN END SUBROUTINE TSCPSS(S,MAXM,MAXN,M,N,P,METHOD,GLOBAL,EPSM,FCQ, + Y,W,FQT,AL2NRM,QHAT,ANLS,DN,FQQ,PTILDA, + CURPOS,PBAR,ZERO1,IERR,RESNEW,FLAG) INTEGER MAXM,MAXN,M,N,P,FLAG,ZERO1,GLOBAL,IERR DOUBLE PRECISION EPSM,RESNEW INTEGER METHOD,PTILDA(N),CURPOS(N),PBAR(N) DOUBLE PRECISION S(MAXN,P),FCQ(M) DOUBLE PRECISION Y(N),W(M),FQT(M),AL2NRM(N) DOUBLE PRECISION QHAT(MAXM,N),ANLS(MAXM,P) DOUBLE PRECISION DN(N),FQQ(M) C********************************************************************** C THIS ROUTINE COMPUTES THE STANDARD STEP. NOTE THAT AT THIS STAGE C THE JACOBIAN MATRIX (QHAT) HAS ALREADY BEEN FACTORED FROM COLUMNS 1 C TO N-P DURING THE TENSOR STEP COMPUTATION. THIS ROUTINE FACTORS C THE MATRIX QHAT FROM COLUMN N-P+1 TO N, THEREBY OBTAINING A QR C FACTORIZATION OF THE FULL MATRIX QHAT, THEN COMPUTES THE STANDARD C STEP BY PREMULTIPLYING THE RIGH-HAND SIDE FCQ BY AN ORTHOGONAL C MATRIX AND BY PERFORMING A BACKWARD SOLVE. C********************************************************************** C C C INPUT PARAMETERS : C ----------------- C C S : FACTORED MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS C (OBTAINED FROM TSQLFC SUBROUTINE) C MAXM : LEADING DIMENSION OF QHAT AND ANLS C MAXN : LEADING DIMENSION OF S C M,N : DIMENSIONS OF PROBLEM C P : COLUMN DIMENSION OF MATRIX S C METHOD : METHOD USED : C METHOD = 0 : STANDARD METHOD IS USED C METHOD = 1 : TENSOR METHOD IS USED C GLOBAL : GLOBAL STRATEGY USED C GLOBAL = 0 : LINE SEARCH IS USED C GLOBAL = 1 : 2-DIMENSIONAL TRUST REGION IS USED C EPSM : MACHINE PRECISION C FCQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY AN C ORTHOGONAL MATRIX C CURPOS : PIVOT VECTOR FOR THE FACTORIZATION OF QHAT FROM COLUMN C 1 TO N-P C Y,W,FQT,AL2NRM : WORKING VECTORS C C C INPUT-OUTPUT PARAMETERS : C ------------------------ C C QHAT : FACTORED MATRIX FROM COLUMN 1 TO N-P C ON ENTRY AND FACTORED MATRIX FROM 1 TO N ON EXIT C ANLS : TENSOR TERM MATRIX ON ENTRY AND ANLS MULTIPLIED BY C ORTHOGONAL MATRICES ON EXIT (THIS IS PERFORMED IN THE C CASE WHERE THE GLOBAL STRATEGY USED IS THE 2-DIMENSIONAL C TRUST REGION) C C OUTPUT PARAMETERS : C ------------------- C C DN : STANDARD STEP C FQQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY C ORTHOGONAL MATRICES (THIS IS USED IN THE CASE WHERE C THE GLOBAL STRATEGY USED IS THE 2-DIMENSIONAL C TRUST REGION) C PTILDA: PIVOT VECTOR FOR THE FACTORIZATION OF THE C MATRIX QHAT FROM N-P+1 TO N C PBAR : PIVOT VECTOR FOR THE FACTORIZATION OF THE C TRANSFORMED MATRIX QHAT FROM 1 TO N C IN CASE OF SINGULARITY C ZERO1 : FIRST ZERO COLUMN OF MATRIX QHAT IN CASE OF SINGULARITY C IERR : RETURNED CODE WITH THE FOLLOWING MEANING : C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C IERR = 0 : OTHERWISE C RESNEW: RESIDUAL OF THE STANDARD MODEL C FLAG : RETURNED CODE WITH THE FOLLOWING MEANINGS : C FLAG = 0 : NO SINGULARITY DETECTED C FLAG = 1 : SINGULARITY DETECTED DURING QR FACTORIZATION C OF QHAT FROM COLUMN 1 TO N-P C FLAG = 2 : SINGULARITY DETECTED DURING QR FACTORIZATION C OF QHAT FROM COLUMN N-P+1 TO N C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DCOPY,DSCAL C TENSOLVE ... TSQRFC,TSQMUV,TSQMTS,TSSMRD,TSBSLV,TSPRMV C TENSOLVE ... TSDLOD,TSQMLV,TSCPMU C C ********************************************************************** INTEGER ZEROTM,I,J DOUBLE PRECISION MU,ZERO,ONE DATA ZERO,ONE/0.0D0,1.0D0/ FLAG = 0 c initialization CALL TSDLOD (M+N,ZERO,FQQ,1) CALL DCOPY(M,FCQ,1,W,1) CALL DSCAL(M,-ONE,W,1) c if the Jacobian is singular then compute the Levenberg-Marquardt c step (label 20) IF(IERR.EQ.1) THEN FLAG = 1 GO TO 20 ENDIF c factor the matrix QHAT from column n-p+1 to n CALL TSQRFC(QHAT,MAXM,N,M,N-P+1,N,IERR,EPSM,AL2NRM,PTILDA,ZERO1) IF((M.EQ.N).AND.(IERR.EQ.0)) THEN ZEROTM = ZERO1-1 ELSE ZEROTM = ZERO1 ENDIF c premultiply W by the orthogonal matrix resulting from the QR c factorization of QHAT CALL TSQMUV(QHAT,W,FQQ,MAXM,M,N-P+1,ZEROTM,.FALSE.) IF(METHOD.EQ.1 .AND. GLOBAL.EQ.1) THEN c premultiply ANLS by the orthogonal matrix resulting from the QR c factorization of QHAT CALL TSQMTS(ANLS,QHAT,MAXM,M,N,M,P,N-P+1,FCQ,ZEROTM) ENDIF IF(IERR.EQ.1) THEN FLAG = 2 GO TO 20 ENDIF c computate the residual of the standard model CALL TSSMRD(FQQ,RESNEW,DN,MU,IERR,M,N) c if QHAT is nonsingular perform a backward solve to obtain Y CALL TSBSLV(QHAT,MAXM,M,N,FQQ,Y) c pivot Y CALL TSPRMV(DN,Y,PTILDA,N,0) IF(N .NE. 1) THEN CALL TSPRMV(Y,DN,CURPOS,N,0) c premultiply Y by the orthogonal matrix resulting from the QL c factorization of S CALL TSQMLV(MAXN,N,P,S,Y,DN,.TRUE.) ENDIF IF(GLOBAL.EQ.1) THEN IERR = 0 CALL DSCAL(M,-ONE,FQQ,1) ENDIF RETURN 20 CONTINUE c @ SINGULAR CASE @ c solve ( QHAT-trans QHAT + MU I ) DN = -QHAT-trans W c put the diagonal elements stored in row m+2 of QHAT into their c propre positions and zero out the unwanted portions of QHAT DO 30 J = 1, ZERO1-1 QHAT(J,J) = QHAT(M+2,J) CALL TSDLOD (M+N-J,ZERO,QHAT(J+1,J),1) 30 CONTINUE DO 40 J = ZERO1, N CALL TSDLOD (M+N-ZERO1+1,ZERO,QHAT(ZERO1,J),1) 40 CONTINUE c compute a small perturbation MU CALL TSCPMU(QHAT,MAXM,N,EPSM,MU) c form the augmented matrix QHAT by adding an (n*n) diag(MU) in c the bottom DO 50 I = M+1,M+N QHAT(I,I-M) = MU 50 CONTINUE c factor the transformed matrix QHAT from 1 to n CALL TSQRFC(QHAT,MAXM,N,M+N,1,N,IERR,EPSM,AL2NRM,PBAR,ZERO1) IF(METHOD.EQ.1 .AND. GLOBAL.EQ.1) THEN c premultiply ANLS by the orthogonal matrix resulting from the QR c factorization of QHAT CALL TSQMTS(ANLS,QHAT,MAXM,M+N,N,M,P,1,FCQ,ZERO1) ENDIF c compute the Levenberg-Marquardt step and the residual of the c standard model IF(FLAG.EQ.1) THEN CALL TSQMUV(QHAT,W,FQQ,MAXM,M+N,1,N+1,.FALSE.) CALL TSBSLV(QHAT,MAXM,M+N,N,FQQ,Y) CALL TSPRMV(DN,Y,PBAR,N,0) CALL TSPRMV(Y,DN,CURPOS,N,0) CALL TSQMLV(MAXN,N,P,S,Y,DN,.TRUE.) CALL TSSMRD(FQQ,RESNEW,DN,MU,IERR,M,N) IF(GLOBAL.EQ.1) THEN IERR = 1 CALL DSCAL(M+N,-ONE,FQQ,1) ENDIF RETURN ELSE CALL TSQMUV(QHAT,FQQ,FQT,MAXM,M+N,1,N+1,.FALSE.) CALL TSBSLV(QHAT,MAXM,M+N,N,FQT,DN) CALL TSPRMV(Y,DN,PBAR,N,0) CALL TSPRMV(DN,Y,PTILDA,N,0) CALL TSPRMV(Y,DN,CURPOS,N,0) CALL TSQMLV(MAXN,N,P,S,Y,DN,.TRUE.) CALL TSSMRD(FQT,RESNEW,DN,MU,IERR,M,N) IF(GLOBAL.EQ.1) THEN IERR = 1 CALL DCOPY(M+N,FQT,1,FQQ,1) CALL DSCAL(M+N,-ONE,FQQ,1) ENDIF ENDIF END SUBROUTINE TSDLOD ( N, CONST, X, INCX ) DOUBLE PRECISION CONST INTEGER INCX, N DOUBLE PRECISION X(*) C********************************************************************** C THIS ROUTINE LOADS ELEMENTS OF X WITH CONST. C********************************************************************** C C INPUT PARAMETERS : C ---------------- C C N : DIMENSION OF THE VECTOR X C CONST : CONSTANT VALUE C INCX : INCREMENT C C OUTPUT PARAMETERS : C ------------------ C C X : VECTOR WITH ELEMENTS EQUAL TO CONST C C********************************************************************** DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) INTEGER IX IF (N .GT. 0) THEN IF (INCX .EQ. 1 .AND. CONST .EQ. ZERO) THEN DO 10 IX = 1, N X(IX) = ZERO 10 CONTINUE ELSE DO 20 IX = 1, 1 + (N - 1)*INCX, INCX X(IX) = CONST 20 CONTINUE ENDIF ENDIF END SUBROUTINE TSD1SV(AJA,S,ANLS,FN,X,MAXM,MAXN,M,N,P,Q,EPSM, + WRK1,WRK2,WRK3,PIVOT,D1) INTEGER MAXM,MAXN,M,N,P,Q INTEGER PIVOT(N) DOUBLE PRECISION EPSM DOUBLE PRECISION AJA(MAXM,N),S(MAXN,P),ANLS(MAXM,P),FN(M),X(P) DOUBLE PRECISION WRK1(N),WRK2(N),WRK3(N),D1(N) C********************************************************************* C THIS ROUTINE SOLVES THE FIRST N-Q LINEAR EQUATIONS IN N-P UNKNOWNS C OF THE TENSOR MODEL. C********************************************************************* C C INPUT PARAMETERS : C ---------------- C C AJA : JACOBIAN MATRIX AT CURRENT ITERATE C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS C ANLS: TENSOR TERM MATRIX AT CURRENT ITERATE C FN : FUNCTION VALUE AT CURRENT ITERATE C X : SOLUTION OF THE LOWER M-N+Q QUADRATIC EQUATIONS IN P C UNKNOWNS OF THE TENSOR MODEL C MAXM: LEADING DIMENSION OF AJA AND ANLS C MAXN: LEADING DIMENSION OF S C M,N : DIMENSIONS OF PROBLEM C P : COLUMN DIMENSION OF S AND ANLS C Q : NUMERICAL RANK OF JACOBIAN : C Q > P : JACOBIAN IS SINGULAR C Q = P : OTHERWISE C EPSM: MACHINE PRECISION C WRK1,WRK2,WRK3 : WORKSPACE C C C OUTPUT PARAMETERS : C ------------------ C C PIVOT : PIVOT VECTOR C D1 : SOLUTION OF THE N-Q LINEAR EQUATIONS IN N-P UNKNOWNS OF C THE TENSOR MODEL C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DCOPY C LEVEL 2 BLAS ... DGEMV C TENSOLVE ... TSDLOD,TSSTMX,TSBSLV,TSQRFC,TSPRMV C TENSOLVE ... TSFSLV,TSQMUV C C********************************************************************* INTEGER ZERO1,I,J,IERR,ICOL DOUBLE PRECISION EPSM1,ZERO,HALF,ALPHA,ONE DATA ZERO,ALPHA,HALF,ONE/0.0D0,1.0D-4,0.50D0,1.0D0/ c compute the top right (n-q) x p submatrix of AJA times X CALL DGEMV('N',N-Q,P,ONE,AJA(1,N-P+1),MAXM, + X,1,ZERO,D1,1) c compute S-trans times X CALL TSSTMX(S,X,MAXN,N,P,WRK3,WRK2) c compute 0.5 * (S-trans times X)**2 DO 10 I = 1, P WRK1(I) = HALF * WRK2(I)**2 10 CONTINUE c compute 0.5 * (top (n-q) x p submatrix of ANLS) * c (S-trans times X)**2 CALL DGEMV('N',N-Q,P,ONE,ANLS(1,1),MAXM,WRK1,1,ZERO,WRK2,1) DO 20 I = 1,N-Q WRK1(I) = -FN(I)-D1(I)-WRK2(I) 20 CONTINUE c if the Jacobian is nonsingular then solve for the first c n-p components of the tensor step and return IF(P.EQ.Q) THEN CALL TSBSLV(AJA,MAXM,M,N-P,WRK1,D1) RETURN ENDIF CALL TSDLOD(Q-P,ZERO,WRK2(N-Q+1),1) c copy top left (n-q) x (n-p) submatrix of AJA into bottom of AJA DO 30 J = 1,N-P CALL DCOPY(N-Q,AJA(1,J),1,AJA(M+3,J),1) 30 CONTINUE c copy the transpose of the top left (n-q) x (n-p) submatrix of AJA c into top of AJA DO 50 J = 1,N-Q AJA(J,J) = AJA(M+2,J) DO 40 I = J+1,N-P AJA(I,J) = AJA(J,I) 40 CONTINUE 50 CONTINUE c zero out the upper triangular (n-q) x (n-q) triangular part of c the transpose of the top left (n-q) x (n-p) submatrix of AJA DO 60 J = 1,N-Q CALL TSDLOD(J-1,ZERO,AJA(1,J),1) 60 CONTINUE c factor the transpose of the top left (n-q) x (n-p) submatrix of AJA EPSM1 = EPSM*ALPHA CALL TSQRFC(AJA,MAXM,N-Q,N-P,1,N-Q,IERR,EPSM1,WRK3,PIVOT,ZERO1) IF(IERR .EQ. 0) THEN ICOL = N-Q ELSE ICOL = ZERO1-1 ENDIF CALL TSPRMV(D1,WRK1,PIVOT,N-Q,0) c solve for the first n-p components of the tensor step CALL TSFSLV(AJA,D1,MAXM,N-P,ICOL,WRK2) CALL TSQMUV(AJA,WRK2,D1,MAXM,N-P,1,ZERO1,.TRUE.) c copy the (n-q) x (n-p) submatrix of AJA from bottom back to c top of AJA DO 70 J = 1,N-P CALL DCOPY(N-Q,AJA(M+3,J),1,AJA(1,J),1) 70 CONTINUE RETURN END SUBROUTINE TSDFCN(P,X,G,AJA,ANLS,S,FN,WRK1,WRK2, + WRK3,WRK4,WRK5,MAXM,MAXN,M,N,Q) INTEGER P,MAXM,MAXN,M,N,Q DOUBLE PRECISION X(P),G(P),AJA(MAXM,N),ANLS(MAXM,P),S(MAXN,P) DOUBLE PRECISION FN(M),WRK1(M),WRK2(P),WRK3(P),WRK4(M),WRK5(M) C********************************************************************* C THIS ROUTINE COMPUTES THE ANALYTIC GRADIENT OF THE FUNCTION GIVEN C BY SUBROUTINE TSQFCN. C********************************************************************* C C INPUT PARAMETERS : C ----------------- C C P : COLUMN DIMENSION OF ANLS AND S C X : POINT AT WHICH GRADIENT IS EVALUATED C AJA: JACOBIAN MATRIX AT CURRENT ITERATE C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS C FN : FUNCTION VALUE AT CURRENT ITERATE C WRK1,WRK2,WRK3,WRK4,WRK5 : WORKSPACE C MAXM : LEADING DIMENSION OF AJA AND ANLS C MAXN : LEADING DIMENSION OF S C M,N : DIMENSIONS OF PROBLEM C Q : NUMERICAL RANK OF JACOBIAN : C Q > P : JACOBIAN IS SINGULAR C Q = P : OTHERWISE C C C OUTPUT PARAMETERS : C ----------------- C C G : GRADIENT AT X C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DAXPY,DDOT C LEVEL 2 BLAS ... DGEMV C TENSOLVE ... TSSTMX,TSDLOD C C********************************************************************* INTEGER I,J,K,L DOUBLE PRECISION ZERO,HALF,ONE DOUBLE PRECISION DDOT DATA ZERO,HALF,ONE/0.0D0,0.50D0,1.0D0/ c compute the lower right (m-n+q) x p submatrix of AJA times X CALL DGEMV('N',M-N+Q,P,ONE,AJA(N-Q+1,N-P+1),MAXM, + X,1,ZERO,WRK1,1) c compute S-trans times X CALL TSSTMX(S,X,MAXN,N,P,WRK2,WRK3) c compute 0.5 * (S-trans times X)**2 DO 10 I = 1, P WRK2(I) = HALF * WRK3(I)**2 10 CONTINUE c compute 0.5 * (lower (m-n+q) x p submatrix of ANLS) * c (S-trans times X)**2 CALL DGEMV('N',M-N+Q,P,ONE,ANLS(N-Q+1,1),MAXM, + WRK2,1,ZERO,WRK4,1) DO 20 I = 1,M-N+Q WRK4(I) = WRK4(I)+FN(N-Q+I)+WRK1(I) 20 CONTINUE c compute AJA-trans * WRK4 CALL DGEMV('T',M-N+Q,P,ONE,AJA(N-Q+1,N-P+1),MAXM, + WRK4,1,ZERO,WRK1,1) c compute ANLS-trans * WRK4 CALL DGEMV('T',M-N+Q,P,ONE,ANLS(N-Q+1,1),MAXM, + WRK4,1,ZERO,WRK5,1) c compute S * diag(S-trans * WRK3) * WRK5 CALL TSDLOD(P,ZERO,WRK2,1) L = P+1 DO 50 J = 1,P L = L-1 WRK2(L) = S(N+2,L) DO 30 I = L+1,P WRK2(I) = S(N-P+J,I) 30 CONTINUE DO 40 K = 1,P WRK2(K) = WRK2(K)*WRK3(K) 40 CONTINUE G(J) = DDOT(P,WRK2,1,WRK5,1) 50 CONTINUE CALL DAXPY(P,ONE,WRK1,1,G,1) RETURN END SUBROUTINE TSDFLT(M,N,ITNLIM,JACFLG,GRADTL,STEPTL,FTOL,METHOD, + GLOBAL,STEPMX,DLT,TYPX,TYPF,IPR,MSG) INTEGER M,N,ITNLIM,JACFLG,METHOD,GLOBAL,MSG,IPR DOUBLE PRECISION GRADTL,STEPTL,FTOL,STEPMX,DLT DOUBLE PRECISION TYPX(N),TYPF(M) C********************************************************************* C THIS ROUTINE SETS DEFAULT VALUES FOR EACH INPUT VARIABLE TO THE C NONLINEAR EQUATION ALGORITHM. C********************************************************************* C C SUBPROGRAMS CALLED: C C TENSOLVE ... TSDLOD C UNCMIN ... DPMEPS C C********************************************************************** DOUBLE PRECISION EPS,DPMEPS,ONE,TWO,THREE,THOUS DATA ONE,TWO,THREE,THOUS/1.0D0,2.0D0,3.0D0,1000.0D0/ JACFLG = 0 EPS = DPMEPS() GRADTL = EPS**(ONE/THREE) STEPTL = EPS**(TWO/THREE) FTOL = EPS**(TWO/THREE) ITNLIM = 150 METHOD = 1 GLOBAL = 0 STEPMX = THOUS DLT = -ONE MSG = 0 IPR = 6 CALL TSDLOD(N,ONE,TYPX,1) CALL TSDLOD(M,ONE,TYPF,1) RETURN END SUBROUTINE TSDUMJ(X,AJA,NR,M,N) INTEGER NR, M, N DOUBLE PRECISION AJA(NR,N),X(N) C********************************************************************* C THIS IS A DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC C WHEN SPECIFIC ANALYTIC JACOBIAN IS NOT SUPPLIED. C********************************************************************* C C INPUT PARAMETERS: C ----------------- C C X : POINT AT WHICH JACOBIAN IS EVALUATED C AJA : JACOBIAN MATRIX C NR : LEADING DIMENSION OF AJA C M,N : DIMENSIONS OF PROBLEM C C*********************************************************************** RETURN END FUNCTION TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,ALPHA,DLT, + NR,M,N,P,NWTAKE,IERR,VN) INTEGER NR,M,N,P,IERR DOUBLE PRECISION ALPHA,DLT,TSFAFA DOUBLE PRECISION ADT(N),AG(N),CONST1(P),CONST2(P) DOUBLE PRECISION FQ(M),VN(M),ANLS(NR,P) LOGICAL NWTAKE C******************************************************************** C THIS FUNCTION COMPUTES || F + J*D + 0.5*A*D**2 ||**2 IN THE C L2 NORM SENS, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2). C******************************************************************** C C C INPUT PARAMETERS C ---------------- C C ANLS : TENSOR TERM MATRIX C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY C ORTHOGONAL MATRICES C ADT : JACOBIAN MATRIX TIMES DT C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) C CONST1 : SHAT-TRANS TIMES DT C CONST2 : SHAT-TRANS TIMES GBAR C ALPHA : POINT AT WHICH TO EVALUATE THE FUNCTION TSFAFA C DLT : CURRENT TRUST RADIUS C NR : LEADING DIMENSION OF THE JACOBIAN C M,N : DIMENSIONS OF THE PROBLEM C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C C C OUTPUT PARAMETERS C ----------------- C C VN : F + J*D + 0.5*A*D**2 C TSFAFA : || F + J*D + 0.5*A*D**2 ||**2 C WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DDOT C TENSOLVE ... TSMAFA C C******************************************************************** INTEGER LEN DOUBLE PRECISION DDOT DOUBLE PRECISION HALF DATA HALF/0.50D0/ CALL TSMAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,ALPHA,DLT, + NR,M,N,P,NWTAKE,IERR,VN) LEN = M IF(IERR.GT.0) LEN = M + N TSFAFA = HALF*DDOT(LEN,VN,1,VN,1) RETURN END SUBROUTINE TSFDFJ(XC,FC,NR,M,N,EPSM,FVEC,FHAT,AJA) INTEGER NR,M,N DOUBLE PRECISION EPSM DOUBLE PRECISION AJA(NR,N),FHAT(M),XC(N),FC(M) EXTERNAL FVEC C*********************************************************************** C THIS ROUTINE COMPUTES THE FINITE DIFFERENCE JACOBIAN AT THE CURRENT C ITERATE XC. C*********************************************************************** C C INPUT PARAMETERS : C ---------------- C C XC : CURRENT ITERATE C FC : FUNCTION VALUE AT XC C NR : LEADING DIMENSION OF AJA C M,N : DIMENSIONS OF PROBLEM C EPSM : MACHINE PRECISION C FVEC : SUBROUTINE TO EVALUATE THE USER'S FUNCTION C FHAT : WORKSPACE C C OUTPUT PARAMETERS : C -------------------- C C AJA : FINITE DIFFERENCE JACOBIAN AT XC C C SUBPROGRAMS CALLED: C C USER ... FVEC C C*********************************************************************** INTEGER I,J DOUBLE PRECISION NDIGIT,RNOISE,STEPSZ,XTMPJ DOUBLE PRECISION SQRTR,RSTPSZ,ONE,TEN INTRINSIC ABS,MAX,SQRT DATA ONE,TEN/1.0D0,10.0D0/ NDIGIT = -LOG10(EPSM) RNOISE = MAX(TEN**(-NDIGIT),EPSM) SQRTR = SQRT(RNOISE) DO 20 J = 1,N STEPSZ = SQRTR*MAX(ABS(XC(J)),ONE) XTMPJ = XC(J) XC(J) = XTMPJ+STEPSZ CALL FVEC(XC,FHAT,M,N) XC(J) = XTMPJ RSTPSZ = ONE/STEPSZ DO 10 I = 1,M AJA(I,J) = (FHAT(I)-FC(I))*RSTPSZ 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE TSFRMT(SHAT,S,AJA,FV,FN,MAXM,MAXN,MAXP,M,N,P,IDP, + AM,X,B,SCALE,A) INTEGER MAXM,MAXN,MAXP,M,N,P INTEGER IDP(P) DOUBLE PRECISION A(MAXM,P),SHAT(MAXN,P),S(MAXN,P),AJA(MAXM,N) DOUBLE PRECISION FV(MAXM,P),FN(M),AM(MAXP,P),X(P),B(P),SCALE(P) C********************************************************************* C THIS ROUTINE FORM THE TENSOR TERM MATRIX OF THE TENSOR MODEL. C********************************************************************* C C INPUT PARAMETERS : C ---------------- C C SHAT: MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS C S : MATRIX OF PREVIOUS DIRECTIONS C AJA : JACOBIAN MATRIX AT CURRENT ITERATE C FV : MATRIX OF PAST FUNCTION VALUES C FN : FUNCTION VALUE AT CURRENT ITERATE C MAXM: LEADING DIMENSION OF AJA, ANLS, AND FV C MAXN: LEADING DIMENSION OF S AND SHAT C MAXP: LEADING DIMENSION OF AM C M : ROW DIMENSION OF MATRICES A,FV,AND AJA C N : COLUMN DIMENSION OF JACOBIAN MATRIX C P : COLUMN DIMENSION OF MATRIX SHAT C IDP : VECTOR WHICH KEEPS TRACK OF LINEARLY INDEPENDENT C DIRECTION POSITIONS WITHIN THE MATRIX S C AM,X,B,SCALE,: WORKSPACE C C OUTPUT PARAMETERS : C ------------------ C C A : TENSOR TERM MATRIX C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DDOT,DNRM2,DSCAL C UNCMIN ... CHOLDC,LLTSLV C C********************************************************************* INTEGER I,J,JJ DOUBLE PRECISION SUM,SC,TOL,DIAGMX,ADDMAX DOUBLE PRECISION ZERO,ONE,TWO DOUBLE PRECISION DDOT,DNRM2 DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ c scale the matrix SHAT and save scaling in SCALE DO 10 J = 1,P SC = ONE/DNRM2(N,SHAT(1,J),1) CALL DSCAL(N,SC,SHAT(1,J),1) SCALE(J) = SC**2 10 CONTINUE c form the matrix AM = (Si Sj)**2 DO 30 J = 1,P JJ = IDP(J) DO 20 I = 1,P AM(I,J) = DDOT(N,S(1,IDP(I)),1,S(1,JJ),1)**2 20 CONTINUE 30 CONTINUE c scale the matrix AM DO 50 I = 1,P DO 40 J = 1,P AM(I,J) = SCALE(I)*SCALE(J)*AM(I,J) 40 CONTINUE 50 CONTINUE c perform a Cholesky decomposition of AM TOL = ZERO DIAGMX = ZERO CALL CHOLDC(MAXP,P,AM,DIAGMX,TOL,ADDMAX) c form the tensor term matrix A DO 70 I = 1,M DO 60 J = 1,P JJ = IDP(J) SUM = DDOT(N,AJA(I,1),MAXM,S(1,JJ),1) B(J) = TWO*(FV(I,JJ) - FN(I) - SUM) B(J) = SCALE(J)*B(J) 60 CONTINUE c solve AM*X = B CALL LLTSLV(MAXP,P,AM,X,B) c copy X into row i of A CALL DCOPY(P,X,1,A(I,1),MAXM) 70 CONTINUE RETURN END SUBROUTINE TSFSCL(X,DX,DF,M,N,FVEC,F) INTEGER M,N DOUBLE PRECISION X(N),DX(N),F(M),DF(M) EXTERNAL FVEC C******************************************************************** C THIS ROUTINE EVALUATES THE FUNCTION AT THE CURRENT ITERATE X THEN C SCALES ITS VALUE. C******************************************************************** C C INPUT PARAMETERS : C ----------------- C C X : CURRENT ITERATE C DX : DIAGONAL SCALING MATRIX FOR X C DF : DIAGONAL SCALING MATRIX FOR F C M,N : DIMENSIONS OF PROBLEM C FVEC : SUBROUTINE TO EVALUATE FUNCTION C C C OUTPUT PARAMETERS : C ----------------- C C F : SCALED FUNCTION VALUE AT CURRENT ITERATE X C C SUBPROGRAMS CALLED: C C TENSOLVE ... TSUNSX,TSSCLF,TSSCLX C USER ... FVEC C C******************************************************************** CALL TSUNSX(X,DX,N) CALL FVEC(X,F,M,N) CALL TSSCLF(F,DF,M) CALL TSSCLX(X,DX,N) RETURN END SUBROUTINE TSFSLV(L,B,NR,M,N,Y) INTEGER NR,M,N DOUBLE PRECISION B(N),L(NR,N),Y(N) C******************************************************************** C THIS ROUTINE DOES A FORWARD SOLVE. C******************************************************************** C C INPUT PARAMETERS : C ----------------- C C L : THE TRANSPOSE OF THE UPPER TRIANGULAR MATRIX OBTAINED C FROM A QR FACTORIZATION OF AN M BY N MATRIX A. DIAG(L) C IS STORED IN ROW M+2. THIS IS THE STORAGE SCHEME USED C IN STEWART, G. W., III(1973) "INTRODUCTION TO MATRIX C COMPUTATION", ACADEMIC PRESS,NEW YORK C B : RIGHT HAND SIDE C NR : LEADING DIMENSION OF MATRIX A C M : ROW DIMENSION OF MATRIX A C N : COLUMN DIMENSION OF MATRIX A C C OUTPUT PARAMETERS : C ------------------ C C Y : VECTOR SOLUTION ON EXIT C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DDOT C C********************************************************************* INTEGER J DOUBLE PRECISION S DOUBLE PRECISION DDOT c solve L Y = B Y(1) = B(1) / L(M+2,1) IF(N .GT. 1) THEN S = L(1,2) * Y(1) Y(2) = (B(2) - S) / L(M+2,2) DO 10 J = 3,N S = DDOT(J-1,L(1,J),1,Y,1) Y(J) = (B(J) - S) / L(M+2,J) 10 CONTINUE ENDIF RETURN END SUBROUTINE TSJMUV(ITN,METHOD,V,CURPOS,PIVOT,PBAR,AJA,SHAT, + FLAG,IERR,MAXM,MAXN,M,N,P,WRK1,WRK2,VN,AV) INTEGER MAXM,MAXN,M,N,P,IERR,ITN,METHOD,FLAG INTEGER CURPOS(N),PIVOT(N),PBAR(N) DOUBLE PRECISION V(N),WRK1(N),WRK2(N),VN(M),AJA(MAXM,N) DOUBLE PRECISION AV(N),SHAT(MAXN,P) C**************************************************************** C THIS ROUTINE CALCULATES THE PRODUCT JACOBIAN TIMES A VECTOR. C**************************************************************** C C INPUT PARAMETERS C ---------------- C C ITN : CURRENT ITERATION NUMBER C METHOD : METHOD TO BE USED C V : VECTOR TO BE MULTIPLIED BY AJA C CURPOS : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA C FROM COLUMN 1 TO N-P) C PIVOT : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA C FROM COLUMN N-P+1 TO N) C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA C IF IT IS SINGULAR C AJA : JACOBIAN MATRIX AT CURRENT ITERATE C SHAT : MATRIX OF LINEARLY INDEPENDENT DIRECTIONS AFTER C A QL FACTORIZATION C FLAG : RETURN CODE WITH THE FOLLOWING MEANINGS: C FLAG = 0 : NO SINGULARITY DETECTED DURING FACTORIZATION C OF THE JACOBIAN FROM COLUMN 1 TO N C FLAG = 1 : SINGULARITY DETECTED DURING FACTORIZATION C OF THE JACOBIAN FROM COLUMN 1 TO N-P C FLAG = 2 : SINGULARITY DETECTED DURING FACTORIZATION C OF THE JACOBIAN FROM COLUMN N-P+1 TO N C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C MAXM : LEADING DIMENSION OF AJA C MAXN : LEADING DIMENSION OF SHAT C M,N : DIMENSIONS OF THE PROBLEM C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS C WRK1,WRK2,VN : WORKSPACE VECTORS C C OUTPUT PARAMETERS C ----------------- C C AV : JACOBIAN TIMES V C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DCOPY C TENSOLVE ... TSPRMV,TSQMLV,TSUTMD C C ********************************************************************** INTEGER LEN IF(ITN .EQ. 1 .OR. METHOD .EQ. 0) THEN CALL TSPRMV(WRK1,V,PIVOT,N,1) IF(IERR .EQ. 1) THEN CALL TSPRMV(WRK2,WRK1,PBAR,N,1) CALL DCOPY(N,WRK2,1,WRK1,1) ENDIF ELSEIF(N .EQ. 1) THEN VN(1) = V(1) ELSE CALL TSQMLV(MAXN,N,P,SHAT,V,VN,.FALSE.) CALL TSPRMV(WRK2,VN,CURPOS,N,1) IF(FLAG .EQ. 0) THEN CALL TSPRMV(WRK1,WRK2,PIVOT,N,1) ELSEIF(FLAG .EQ. 1) THEN CALL TSPRMV(WRK1,WRK2,PBAR,N,1) ELSEIF(FLAG .EQ. 2 ) THEN CALL TSPRMV(WRK1,WRK2,PIVOT,N,1) CALL TSPRMV(WRK2,WRK1,PBAR,N,1) CALL DCOPY(N,WRK2,1,WRK1,1) ENDIF ENDIF LEN = M IF(IERR .GT. 0) LEN = M + N CALL TSUTMD(AJA,WRK1,MAXM,LEN,N,AV) RETURN END SUBROUTINE TSJQTP(Q,MAXM,MAXN,N,M,P,WRK1,WRK2,AJA) INTEGER MAXM,MAXN,N,M,P DOUBLE PRECISION AJA(MAXM,N),Q(MAXN,P),WRK1(N),WRK2(N) C********************************************************************** C THIS ROUTINE GETS J*(Q-TRANS) BY COMPUTING EACH ROW OF THE C RESULTING MATRIX AS FOLLOWS : (J*Q-TRANS)I-TH ROW<--Q*(J)I-TH ROW. C********************************************************************** C C INPUT PARAMETERS : C ----------------- C C Q : RESULTING MATRIX FROM A QL FACTORIZATION C MAXM : LEADING DIMENSION OF AJA C MAXN : LEADING DIMENSION OF Q C M,N : DIMENSIONS OF PROBLEM C P : COLUMN DIMENSION OF MATRIX Q C WRK1,WRK2: WORKING VECTOR C C INPUT-OUTPUT PARAMETERS : C ------------------------ C C AJA : JACOBIAN MATRIX ON ENTRY AND JACOBIAN MULTIPLIED BY THE C ORTHOGONAL MATRIX Q ON EXIT C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DCOPY C TENSOLVE ... TSQMLV C C********************************************************************** INTEGER I DO 30 I = 1,M c copy the i-th row of AJA into WRK1 CALL DCOPY(N,AJA(I,1),MAXM,WRK1,1) CALL TSQMLV(MAXN,N,P,Q,WRK1,WRK2,.FALSE.) c form the i-th row of AJA*(Q-trans) CALL DCOPY(N,WRK2,1,AJA(I,1),MAXM) 30 CONTINUE RETURN END SUBROUTINE TSLMIN(XC,XP,P1,Q,ANLS,FQ,ADT,AG,CONST1,CONST2, + DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP,VNS,XPLUS) INTEGER NR,M,N,P,IERR DOUBLE PRECISION XC,XP,XPLUS,P1,Q,DLT,TOL DOUBLE PRECISION ADT(N),AG(N),VN(M),VNP(M),VNS(M) DOUBLE PRECISION ANLS(NR,P),FQ(M),CONST1(P),CONST2(P) LOGICAL NWTAKE C*********************************************************************** C THIS ROUTINE FINDS A LOCAL MINIMIZER OF A ONE-VARIABLE FUNCTION IN AN C INTERVAL [XC XP]. C*********************************************************************** C C INPUT PARAMETERS : C ----------------- C C XC,XP : LOWER AND UPPER BOUND OF INTERVAL IN WHICH THE SEARCH C IS PERFORMED C P1,Q : FIRST DERIVATIVES OF THE ONE-VARIABLE FUNCTION C ANLS : TENSOR TERM MATRIX C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY C ORTHOGONAL MATRICES C ADT : JACOBIAN TIMES THE STEP DT (SEE SUBROUTINE TS2DTR) C AG : JACOBIAN TIMES THE GRADIENT G (SEE SUBROUTINE TS2DTR) C CONST1 : SHAT-TRANS * DT (SEE SUBROUTINE TS2DTR) C CONST2 : SHAT-TRANS * GBAR (SEE SUBROUTINE TS2DTR) C DLT : TRUST RADIUS C NR : LEADING DIMENSION OF ANLS MATRIX C M,N : DIMENSIONS OF PROBLEM C P : COLUMN DIMENSION OF MATRIX ANLS C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : OTHERWISE C TOL : SMALL TOLERANCE C VN,VNP,VNS : WORKING VECTORS C C C OUTPUT PARAMETERS : C ----------------- C C XPLUS : LOCAL MINIMIZER OF THE ONE-VARIABLE FUNCTION C C SUBPROGRAMS CALLED : C C TENSOLVE ... TSMSDA,TSFAFA,TSLMSP,TSMFDA C C*********************************************************************** INTEGER ITERCD,RETCD,ITNCNT DOUBLE PRECISION ALEFT,ARIGHT,T,E,TSMSDA,S,SINIT,TSFAFA,TSMFDA DOUBLE PRECISION ZERO,OTT,TWO,SMALL LOGICAL SKIP INTRINSIC ABS,MIN,MAX DATA ZERO,OTT,TWO,SMALL/0.0D0,1.0D-4,2.0D0,2.0D-20/ RETCD = 0 ALEFT = MIN(XC,XP) ARIGHT = MAX(XC,XP) ITNCNT = 0 T = ABS(XC-XP) SKIP = .FALSE. c compute the second derivative value at the current point E = TSMSDA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC,DLT, + NR,M,N,P,NWTAKE,IERR,SKIP,VN,VNP,VNS) 10 IF(E.GT.ZERO) THEN S = -P1/E IF(ABS(S).GT.TWO*T) THEN IF (S.LT.ZERO) THEN S = -TWO*T ELSE S = TWO*T ENDIF ENDIF ELSE IF (P1.GT.ZERO) THEN S = -T ELSE S = T ENDIF ENDIF IF(XC+S.GT.ARIGHT) S = ARIGHT-XC IF(XC+S.LT.ALEFT) S = ALEFT-XC SINIT = ABS(S) 20 CONTINUE c compute a next iterate XPLUS IF (TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC+S,DLT, + NR,M,N,P,NWTAKE,IERR,VN).GT.Q + OTT*S*P1) THEN S = S/2 IF(ABS(S).LT.SMALL*SINIT.OR.S.EQ.ZERO) THEN RETCD = 1 ELSE GO TO 20 ENDIF ENDIF XPLUS = XC+S ITNCNT = ITNCNT+1 c check stopping criteria CALL TSLMSP(XC,XPLUS,ITNCNT,RETCD,ITERCD,ANLS,ADT,AG, + CONST1,CONST2,DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP) IF(ITERCD.GT.0) RETURN c update XC XC = XPLUS c compute function and derivative values at the new point Q = TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC,DLT, + NR,M,N,P,NWTAKE,IERR,VN) P1 = TSMFDA(ANLS,ADT,AG,CONST1,CONST2,XC,DLT, + NR,M,N,P,NWTAKE,IERR,VN,VNP) SKIP = .TRUE. E = TSMSDA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC,DLT, + NR,M,N,P,NWTAKE,IERR,SKIP,VN,VNP,VNS) GO TO 10 END SUBROUTINE TSLMSP(XC,XP,ITNCNT,RETCD,ITERCD,ANLS,ADT,AG,CONST1, + CONST2,DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP) INTEGER NR,M,N,P,IERR,RETCD,ITERCD,ITNCNT DOUBLE PRECISION XC,XP,DLT,TOL DOUBLE PRECISION ADT(N),AG(N),CONST1(P) DOUBLE PRECISION CONST2(P),VN(M),VNP(M),ANLS(NR,P) LOGICAL NWTAKE C*********************************************************************** C THIS ROUTINE CHECKS THE STOPPING CRITERIA FOR A LOCAL MINIMIZER. C*********************************************************************** C C INPUT PARAMETERS : C ----------------- C C XC : CURRENT ITERATE (FROM SEARCH SUBROUTINE) C XP : NEXT ITERATE (FROM SEARCH SUBROUTINE) C ITNCNT : ITERATION LIMIT C RETCD : RETURN CODE FROM LINE SEARCH C DLT : TRUST RADIUS C AJA : JACOBIAN AT THE CURRENT ITERATE C NR : LEADING DIMENSION OF THE JACOBIAN MATRIX C M,N : DIMENSIONS OF THE PROBLEM C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS : C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE : C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : OTHERWISE C TOL : SMALL TOLERANCE C METHOD : METHOD TO USE C = 0 : STANDARD METHOD USED C = 1 : TENSOR METHOD USED C VN,VNP : WORKING VECTORS C C C OUTPUT PARAMETERS : C ------------------ C C ITERCD : RETURN CODE WITH FOLLOWING MEANINGS : C ITERCD = 1 : FIRST DERIVATIVE AT THE CURRENT POINT C CLOSE TO 0 C ITERCD = 2 : SUCCESSIVE ITERATES WITHIN TOLERANCE C ITERCD = 3 : LINE SEARCH FAILED TO LOCATE A POINT C LOWER THAT THE CURRENT POINT C ITERCD = 4 : ITERATION LIMIT EXCEEDED C C*********************************************************************** DOUBLE PRECISION TSMFDA,GRDT,ZERO INTRINSIC ABS,SQRT DATA ZERO/0.0D0/ GRDT = SQRT(TOL) ITERCD = 0 IF(RETCD.EQ.1) THEN ITERCD = 3 ELSEIF(ABS(TSMFDA(ANLS,ADT,AG,CONST1,CONST2,XP,DLT, + NR,M,N,P,NWTAKE,IERR,VN,VNP)) .LT. GRDT) THEN ITERCD = 1 ELSEIF(XP.NE.ZERO .AND. ABS(XP-XC)/ABS(XP).LE.TOL) THEN ITERCD = 2 ELSEIF(ITNCNT.GE.150) THEN ITERCD = 4 ENDIF RETURN END SUBROUTINE TSLSCH(M,N,XC,D,G,STEPTL,DX,DF,FVEC, + MXTAKE,STEPMX,XP,FP,FCNORM,FPNORM,RETCD) INTEGER M,N,RETCD DOUBLE PRECISION STEPTL,FCNORM,FPNORM DOUBLE PRECISION XC(N) DOUBLE PRECISION D(N),G(N),XP(N),FP(M) DOUBLE PRECISION DX(N),DF(M),STEPMX LOGICAL MXTAKE EXTERNAL FVEC C********************************************************************** C THIS ROUTINE FINDS A NEXT ITERATE USING A STANDARD LINE SEARCH METHOD. C********************************************************************** C C INPUT PARAMETERS : C ----------------- C C M,N : DIMENSIONS OF PROBLEM C XC : CURRENT ITERATE C D : SEARCH DIRECTION C G : GRADIENT AT CURRENT ITERATE C STEPTL : RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES C ARE CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM C DX : DIAGONAL SCALING MATRIX FOR X C DF : DIAGONAL SCALING MATRIX FOR F C FVEC: SUBROUTINE TO EVALUATE THE FUNCTION C STEPMX: MAXIMUM ALLOWABLE STEP SIZE C C C OUTPUT PARAMETERS : C ----------------- C C MXTAKE: BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED C XP : NEXT ITARATE C FP : FUNCTION VALUE AT NEXT ITERATE C FCNORM : 0.5 * || F(XC) ||**2 C FPNORM : 0.5 * || F(XP) ||**2 C RETCD : RETURN CODE WITH THE FOLLOWING MEANING : C RETCD = 0 : SATISFACTORY LOCATION OF A NEW ITERATE C RETCD = 1 : NO SATISFACTORY POINT FOUND SUFFICIENTLY C DISTINCT FROM X C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DDOT,DNRM2 C TENSOLVE ... TSFSCL C USER ... FVEC C C********************************************************************** INTEGER I DOUBLE PRECISION ALPHA,SLOPE,RELENG DOUBLE PRECISION TEMP1,TEMP2,ALMDA,TEMP,ALMDAT,ALMDAM DOUBLE PRECISION SLN,SCL DOUBLE PRECISION DDOT,DNRM2 DOUBLE PRECISION ZERO,TENTH,HALF,Z99,ONE,TWO,TEN INTRINSIC ABS PARAMETER (ALPHA = 1.0D-4) DATA ZERO,TENTH,HALF,Z99,ONE,TWO,TEN/0.0D0,0.10D0,0.50D0,0.99D0, + 1.0D0,2.0D0,10.0D0/ MXTAKE = .FALSE. SLN = DNRM2(N,D,1) IF(SLN .GT. STEPMX) THEN c step longer than maximum allowed SCL = STEPMX/SLN CALL DSCAL(N,SCL,D,1) SLN = STEPMX ENDIF c compute SLOPE = G-trans * D SLOPE = DDOT(N,G,1,D,1) c initialization of RETCD RETCD = 0 c compute the smallest value allowable for the damping c parameter ALMDA, i.e, ALMDAM RELENG = ZERO DO 20 I = 1,N TEMP1 = MAX(ABS(XC(I)), ONE) TEMP2 = ABS(D(I))/TEMP1 RELENG = MAX(RELENG,TEMP2) 20 CONTINUE ALMDAM = STEPTL/RELENG ALMDA = ONE 40 CONTINUE c compute the next iterate XP DO 50 I = 1,N XP(I) = XC(I)+ALMDA*D(I) 50 CONTINUE c evaluate the objective function at XP and its residual CALL TSFSCL(XP,DX,DF,M,N,FVEC,FP) FPNORM = HALF*DNRM2(M,FP,1)**2 c test whether the full step produces enough decrease in c the l2 norm of the objective function. If not update ALMDA c and compute a new step IF (FPNORM.GT.(FCNORM + (ALPHA* ALMDA * SLOPE))) THEN ALMDAT = ((-ALMDA**2)*SLOPE)/(TWO*(FPNORM-FCNORM-ALMDA*SLOPE)) TEMP = ALMDA/TEN ALMDA = MAX(TEMP,ALMDAT) IF(ALMDA.LT.ALMDAM) THEN RETCD = 1 RETURN ENDIF GO TO 40 ELSE IF(ALMDA.EQ.TENTH .AND. SLN.GT.Z99*STEPMX) MXTAKE=.TRUE. ENDIF RETURN END SUBROUTINE TSMAFA(ANLS,F,ADT,AG,CONST1,CONST2,ALPHA,DLT, + NR,M,N,P,NWTAKE,IERR,VN) INTEGER NR,M,N,P,IERR DOUBLE PRECISION ALPHA,DLT DOUBLE PRECISION ADT(N),AG(N),CONST1(P) DOUBLE PRECISION CONST2(P),F(M),VN(M),ANLS(NR,P) LOGICAL NWTAKE C*********************************************************************** C THIS ROUTINE COMPUTES THE VECTOR VN = F(XC) + J(XC)*D + 0.5*A*D**2, C WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2). C*********************************************************************** C C C INPUT PARAMETERS : C ----------------- C C ANLS : TENSOR TERM MATRIX C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) C CONST1: SHAT-TRANS * DT (SEE SUBROUTINE TS2DTR) C CONST2: SHAT-TRABS * GBAR (SEE SUBROUTINE TS2DTR) C ALPHA : POINT AT WHICH DERIVATIVE IS EVALUATED C DLT : CURRENT TRUST RADIUS C NR : LEADING DIMENSION OF ANLS C M,N : DIMENSIONS OF THE PROBLEM C P : COLUMN DIMENSION OF THE MATRIX ANLS C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE : C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C C OUTPUT PARAMETERS : C ------------------- C C VN : F + J*D + 0.5*A*D**2, WHERE C D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DAXPY C TENSOLVE ... TSDLOD C C******************************************************************* INTEGER I,J,LEN DOUBLE PRECISION EXPR,CONST,ZERO,HALF INTRINSIC SQRT DATA ZERO,HALF/0.0D0,0.50D0/ EXPR = SQRT(DLT**2 - ALPHA**2) DO 10 I = 1,N VN(I) = ALPHA*ADT(I) + EXPR*AG(I) 10 CONTINUE CALL TSDLOD (M,ZERO,VN(N+1),1) LEN = M IF(IERR .GT. 0) LEN = M + N DO 30 I = 1, LEN VN(I) = VN(I) + F(I) 30 CONTINUE IF(NWTAKE) RETURN DO 70 J = 1,P CONST = HALF*(ALPHA*CONST1(J) + EXPR*CONST2(J))**2 CALL DAXPY(LEN,CONST,ANLS(1,J),1,VN,1) 70 CONTINUE RETURN END SUBROUTINE TSMDLS(AJA,SHAT,ANLS,XC,M,N,MAXM,MAXN,P,DT,G, + DX,DF,FVEC,METHOD,STEPTL,GLOBAL,STEPMX, + EPSM,FQ,WRK1,WRK2,WRK3,WRK4,DN,FQQ,PIVOT, + CURPOS,PBAR,MXTAKE,XP,FP,FCNORM,FPNORM, + ZERO1,RETCD,IERR) INTEGER M,N,MAXM,MAXN,P,METHOD,GLOBAL,ZERO1,RETCD,IERR INTEGER PIVOT(N),PBAR(N),CURPOS(N) DOUBLE PRECISION STEPTL,STEPMX,EPSM,FCNORM,FPNORM DOUBLE PRECISION AJA(MAXM,N),SHAT(MAXN,P),ANLS(MAXM,P) DOUBLE PRECISION XC(N),DT(N),G(N),DX(N),DF(M),FQ(M) DOUBLE PRECISION WRK1(M),WRK2(M),WRK3(M),WRK4(N) DOUBLE PRECISION DN(N),FQQ(M),XP(N),FP(M) LOGICAL MXTAKE EXTERNAL FVEC C********************************************************************** C THIS ROUTINE FINDS A NEXT ITERATE USING A LINE SEARCH METHOD. IT C TRIES THE FULL TENSOR STEP FIRST. IF THIS IS NOT SUCCESSFUL THEN C IT COMPUTES THE STANDARD DIRECTION AND COMPUTES A STEP IN THAT C DIRECTION. NEXT, IF THE TENSOR DIRECTION IS DESCENT, IT COMPUTES C A STEP IN THE TENSOR DIRECTION. THE ITERATE THAT PRODUCES C THE LOWER RESIDUAL IS THE NEXT ITERATE FOR THE NONLINEAR ALGORITHM. C********************************************************************** C C INPUT PARAMETERS C ---------------- C C AJA : JACOBIAN AT CURRENT ITERATE C SHAT : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS C AFTER A QL FACORIZATION C ANLS : TENSOR TERM MATRIX C XC : CURRENT ITERATE C M,N : DIMENSIONS OF THE PROBLEM C MAXM : LEADING DIMENSION OF AJA AND ANLS C MAXN : LEADING DIMENSION OF SHAT C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS C DT : TENSOR STEP C G : GRADIENT AT CURRENT ITERATE C DX : DIAGONAL SCALING MATRIX FOR X C DF : DIAGONAL SCALING MATRIX FOR F C GBAR : STEEPEST DESCENT DIRECTION (= -G) C METHOD : METHOD TO USE C = 0 : STANDARD METHOD USED C = 1 : TENSOR METHOD USED C STEPTL : STEP TOLERANCE C GLOBAL : GLOBAL STRATEGY USED C = 0 : LINE SEARCH IS USED C = 1 : 2-DIMENSIONAL TRUST REGION IS USED C STEPMX : MAXIMUM ALLOWABLE STEP SIZE C EPSM : MACHINE PRECISION C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY AN C ORTHOGOL MATRIX C WRK1,WRK2,WRK3,WRK4 : WORKING VECTORS C C C OUTPUT PARAMETERS C ----------------- C C DN : NEWTON STEP C FQQ : FQ MULTIPLIED BY AN ORTHOGONAL MATRIX C CURPOS : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE C JACOBIAN FROM COLUMN 1 TO N-P) C PIVOT : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE C JACOBIAN FROM COLUMN N-P+1 TO N) C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE C JACOBIAN IF IT IS SINGULAR C MXTAKE : BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED C XP : NEXT ITERATE C FP : FUNCTION VALUE AT NEXT ITERATE C FCNORM : 0.5 * || F(XC) ||**2 C FPNORM : 0.5 * || F(XP) ||**2 C ZERO1 : FIRST ZERO COLUMN OF THE JACOBIAN IN CASE OF C SINGULARITY C RETCD : RETURN CODE WITH THE FOLLOWING MEANING : C RETCD = 0 : SATISFACTORY LOCATION OF A NEW ITERATE C RETCD = 1 : NO SATISFACTORY POINT FOUND SUFFICIENTLY C DISTINCT FROM X C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DCOPY,DDOT,DNRM2 C TENSOLVE ... TSFSCL,TSCPSS,TSLSCH C C*********************************************************************** INTEGER I,FLAG,RETCD1 DOUBLE PRECISION ALPHA,SLOPE,RELENG DOUBLE PRECISION TEMP1,TEMP2,ALMDA,RESNEW,F1N,DTNORM,GNORM DOUBLE PRECISION SLN,SCL DOUBLE PRECISION BETA,TEMP,ALMDAT,ALMDAM DOUBLE PRECISION DDOT,DNRM2 DOUBLE PRECISION ZERO,TENTH,HALF,Z99,ONE,TWO,TEN INTRINSIC ABS PARAMETER (ALPHA = 1.0D-4) DATA ZERO,TENTH,HALF,Z99,ONE,TWO,TEN/0.0D0,0.10D0,0.50D0,0.99D0, + 1.0D0,2.0D0,10.0D0/ MXTAKE = .FALSE. SLN = DNRM2(N,DT,1) IF(SLN .GT. STEPMX) THEN c step longer than maximum allowed SCL = STEPMX/SLN CALL DSCAL(N,SCL,DT,1) SLN = STEPMX ENDIF c compute SLOPE = G-Trans * DT SLOPE = DDOT(N,G,1,DT,1) c initialization of RETCD RETCD = 0 c compute the smallest value allowable for the damping c parameter ALMDA, i.e, ALMDAM RELENG = ZERO DO 20 I = 1,N TEMP1 = MAX(ABS(XC(I)), ONE) TEMP2 = ABS(DT(I))/TEMP1 RELENG = MAX(RELENG, TEMP2) 20 CONTINUE ALMDAM = STEPTL/RELENG ALMDA = ONE c compute the next iterate XP DO 30 I = 1,N XP(I) = XC(I)+ALMDA*DT(I) 30 CONTINUE c evaluate the objective function at XP and its residual CALL TSFSCL(XP,DX,DF,M,N,FVEC,FP) FPNORM = HALF*DNRM2(M,FP,1)**2 c test whether the full tensor step produces enough decrease in the c l2 norm of of the objective function IF (FPNORM.LT.(FCNORM + (ALPHA* ALMDA * SLOPE))) RETURN c compute the standard direction CALL TSCPSS(SHAT,MAXM,MAXN,M,N,P,METHOD,GLOBAL,EPSM,FQ, + WRK1,WRK2,WRK3,WRK4,AJA,ANLS,DN,FQQ,PIVOT, + CURPOS,PBAR,ZERO1,IERR,RESNEW,FLAG) c compute a step in the standard direction CALL TSLSCH(M,N,XC,DN,G,STEPTL,DX,DF,FVEC, + MXTAKE,STEPMX,WRK1,WRK2,FCNORM,F1N,RETCD1) c test whether the tensor direction is descent DTNORM = DNRM2(N,DT,1) GNORM = DNRM2(N,G,1) IF(M.GT.N) THEN BETA = TENTH ELSE BETA = ALPHA ENDIF TEMP1 = -BETA*DTNORM*GNORM c compute a step in the tensor direction IF(SLOPE .LE. TEMP1) THEN 50 CONTINUE ALMDAT = ((-ALMDA**2)*SLOPE)/(TWO*(FPNORM-FCNORM-ALMDA*SLOPE)) TEMP = ALMDA/TEN ALMDA = MAX(TEMP, ALMDAT) IF(ALMDA.LT.ALMDAM) THEN IF(RETCD1. EQ. 1) THEN RETCD = 1 GO TO 70 ENDIF ENDIF DO 60 I = 1,N XP(I) = XC(I)+ALMDA*DT(I) 60 CONTINUE CALL TSFSCL(XP,DX,DF,M,N,FVEC,FP) FPNORM = HALF*DNRM2(M,FP,1)**2 IF (FPNORM .GT.(FCNORM + (ALPHA* ALMDA * SLOPE))) GO TO 50 IF(ALMDA.EQ.TENTH .AND. SLN.GT.Z99*STEPMX) MXTAKE=.TRUE. 70 CONTINUE c select the next iterate that produces the lower function value IF(F1N .LT. FPNORM) THEN CALL DCOPY(N,WRK1,1,XP,1) CALL DCOPY(M,WRK2,1,FP,1) FPNORM = F1N ENDIF ELSE CALL DCOPY(N,WRK1,1,XP,1) CALL DCOPY(M,WRK2,1,FP,1) FPNORM = F1N ENDIF RETURN END FUNCTION TSMFDA(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, + NR,M,N,P,NWTAKE,IERR,VN,VNP) INTEGER NR,M,N,P,IERR DOUBLE PRECISION ALPHA,DLT,TSMFDA DOUBLE PRECISION ADT(N),AG(N),CONST1(P),CONST2(P),VN(M),VNP(M) DOUBLE PRECISION ANLS(NR,P) LOGICAL NWTAKE C*********************************************************************** C THIS FUNCTION COMPUTES THE DERIVATIVE OF || F + J*D + 0.5*A*D**2 ||**2 C IN THE L2 NORM SENS, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2). C*********************************************************************** C C C INPUT PARAMETERS C ---------------- C C ANLS : TENSOR MATRIX C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY C ORTHOGONAL MATRICES C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) C CONST1 : SHAT-TRANS TIMES DT (SEE SUBROUTINE TS2DTR) C CONST2 : SHAT-TRANS TIMES GBAR (SEE SUBROUTINE TS2DTR) C ALPHA : POINT AT WHICH TO EVALUATE THE DERIVATIVE OF FUNCTION C DLT : CURRENT TRUST RADIUS C NR : LEADING DIMENSION OF THE JACOBIAN C M,N : DIMENSIONS OF THE PROBLEM C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: C IERR=0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR=1 : SINGULARITY OF JACOBIAN DETECTED C C C OUTPUT PARAMETERS C ----------------- C C C VN : F + J*D + 0.5*A*D**2 C VNP : DERIVATIVE IN ALPHA OF F + J*D + 0.5*A*D**2 C TSMFDA : DERIVATIVE IN ALPHA OF || F + J*D + 0.5*A*D**2 ||**2 C WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DDOT C TENSOLVE ... TSMFDV C C*********************************************************************** INTEGER LEN DOUBLE PRECISION DDOT CALL TSMFDV(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, + NR,M,N,P,NWTAKE,IERR,VNP) LEN = M IF(IERR.GT.0) LEN = M + N TSMFDA = DDOT(LEN,VNP,1,VN,1) RETURN END SUBROUTINE TSMFDV(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, + NR,M,N,P,NWTAKE,IERR,VNP) INTEGER NR,M,N,P,IERR DOUBLE PRECISION ALPHA,DLT DOUBLE PRECISION ADT(N),AG(N),CONST1(P) DOUBLE PRECISION CONST2(P),VNP(M),ANLS(NR,P) LOGICAL NWTAKE C*********************************************************************** C THIS ROUTINE COMPUTES THE DERIVATIVE IN ALPHA OF THE VECTOR C VN = F(XC) + J(XC)*D + 0.5*A*D**2, WHERE D = ALPHA*DT + C SQRT(DLT**2-ALPHA**2). C*********************************************************************** C C C INPUT PARAMETERS : C ----------------- C C ANLS : TENSOR TERM MATRIX C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) C CONST1: SHAT-TRANS TIMES DT (SEE SUBROUTINE TS2DTR) C CONST2: SHAT-TRANS TIMES GBAR (SEE SUBROUTINE TS2DTR) C ALPHA : POINT AT WHICH DERIVATIVE IS EVALUATED C DLT : CURRENT TRUST RADIUS C NR : LEADING DIMENSION OF ANLS C M,N : DIMENSIONS OF THE PROBLEM C P : COLUMN DIMENSION OF THE MATRIX ANLS C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS : C NWTAKE = .TRUE. : STANDARD STEP TAKEN C NWTAKE = .FALSE. : TENSOR STEP TAKEN C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED C C C OUTPUT PARAMETERS : C ------------------- C C VNP : THE DERIVATIVE IN ALPHA OF VN = F(XC) + J(XC)*D + C 0.5*A*D**2, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) C C SUBPROGRAMS CALLED: C C LEVEL 1 BLAS ... DAXPY C TENSOLVE ... TSDLOD C C**