|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
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: February, 1994
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
ca tensor method based on a quadratic model and a standard method
cbased 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,TSD1SV,TSDFCN,TSDFLT,TSDUMJ,
c TSFAFA,TSFDFJ,TSFRMT,TSFSCL,TSFSLV,TSJMUV,TSJQTP,TSLMIN,TSLMSP,TSLSCH,
c TSMAFA,TSMDLS,TSMFDA,TSMFDV,TSMGSA,TSMSDA,TSMSDV,TSMSLV,TSNECI,TSNESI,
c TSNESV,TSNSTP,TSPRMV,TSRSLT,TSQ1P1,TSQFCN,TSQLFC,TSQMLV,TSQMTS,TSQMUV,
c TSQRFC,TSRSID,TSSCLF,TSSCLJ,TSSCLX,TSSLCT,TSSMIN,TSSMRD,TSSQP1,TSSSTP,
c 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,DLOAD,DSCAL,DSWAP,IDAMAX
c LEVEL 2 BLAS: DGEMV
c
c Parameters and Default Values:
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 PMAX 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 PMAX->: A positive integer specifying the row dimension of the work
c array WRKUNC in the user's calling program. It must satisfy
c PMAX >= 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 PMAX*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 TSANJA->: 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 TSANJA(JAC, X, MAXM, M, N)
c where X is a vector of length N and the 2-dimensional array JAC 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 TSFVEC->: 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 TSFVEC(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 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, PMAX < 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 TSANJA not
c supplied (incompatible input).
c MSG = -11 : Probable coding error in the user's analytic Jacobian
c routine TSANJA. 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.
SUBROUTINE TS2DTR(AJA,SHAT,ANLS,DT,G,GBAR,XC,FC,METHOD,NWTAKE, TS2DTR
+ STEPMX,STEPTL,EPSM,MXTAKE,DLT,FQ,MAXM,MAXN, TS2DTR
+ M,N,P,CURPOS,PIVOT,PBAR,ITN,IERR,FLAG,DXN, TS2DTR
+ DFN,TSFVEC,D,XPLSP,ADT,AG,TEMP,VN,VNP,VNS, TS2DTR
+ WRK1,CONST1,CONST2,FNORM,XPLS,FP,FPLS,RETCD) TS2DTR
TS2DTR
INTEGER MAXM,MAXN,M,N,P,ITN,METHOD,IERR,FLAG TS2DTR
DOUBLE PRECISION STEPMX,STEPTL,EPSM,DLT,FPLS TS2DTR
INTEGER CURPOS(N),PIVOT(N),PBAR(N),RETCD TS2DTR
DOUBLE PRECISION DT(N),G(N),GBAR(N),XC(N),FC(M) TS2DTR
DOUBLE PRECISION XPLS(N),FP(M),XPLSP(N),AJA(MAXM,N),D(M) TS2DTR
DOUBLE PRECISION TEMP(M),SHAT(MAXN,P),ANLS(MAXM,P),VNS(M) TS2DTR
DOUBLE PRECISION VN(M),VNP(M),FQ(M),DXN(N),DFN(M) TS2DTR
DOUBLE PRECISION ADT(N),AG(N),WRK1(M),CONST1(P),CONST2(P) TS2DTR
LOGICAL NWTAKE,MXTAKE TS2DTR
TS2DTR
C********************************************************************** TS2DTR
C THIS ROUTINE FINDS A NEXT ITERATE BY A 2-DIMENSIONAL TRUST REGION. TS2DTR
C********************************************************************** TS2DTR
C TS2DTR
C TS2DTR
C INPUT PARAMETERS : TS2DTR
C ----------------- TS2DTR
C TS2DTR
C AJA : JACOBIAN AT THE CURRENT ITERATE TS2DTR
C SHAT : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TS2DTR
C AFTER A QL FACTORIZATION TS2DTR
C ANLS : TENSOR TERM MATRIX TS2DTR
C DT : CURRENT STEP TS2DTR
C G : GRADIENT AT CURRENT ITERATE TS2DTR
C GBAR : STEEPEST DESCENT DIRECTION (= -G) TS2DTR
C XC : CURRENT ITERATE TS2DTR
C FC : FUNCTION VALUE AT CURRENT ITERATE TS2DTR
C METHOD : METHOD TO USE TS2DTR
C = 0 : STANDARD METHOD USED TS2DTR
C = 1 : TENSOR METHOD USED TS2DTR
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: TS2DTR
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TS2DTR
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TS2DTR
C STEPMX : MAXIMUM STEP ALLOWED TS2DTR
C STEPTL : STEP TOLERANCE TS2DTR
C EPSM : MACHINE PRECISION TS2DTR
C MXTAKE : BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED TS2DTR
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TS2DTR
C ORTHOGONAL MATRICES TS2DTR
C MAXM : LEADING DIMENSION OF AJA AND ANLS TS2DTR
C MAXN : LEADING DIMENSION OF SHAT TS2DTR
C M,N : DIMENSIONS OF PROBLEM TS2DTR
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TS2DTR
C CURPOS : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE TS2DTR
C JACOBIAN FROM COLUMN 1 TO N-P) TS2DTR
C PIVOT : PIVOT VECTOR ( USED DURING THE FACTORIZATION OF THE TS2DTR
C JACOBIAN FROM COLUMN N-P+1 TO N) TS2DTR
C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE TS2DTR
C JACOBIAN IF IT IS SINGULAR TS2DTR
C FNORM : 0.5 * || FC ||**2 TS2DTR
C ITN : ITERATION NUMBER TS2DTR
C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE: TS2DTR
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TS2DTR
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TS2DTR
C FLAG : RETURN CODE WITH THE FOLLOWING MEANINGS : TS2DTR
C FLAG = 0 : NO SINGULARITY DETECTED DURING TS2DTR
C FACTORIZATION OF THE JACOBIAN FROM TS2DTR
C COLUMN 1 TO N TS2DTR
C FLAG = 1 : SINGULARITY DETECTED DURING FACTORIZATION TS2DTR
C OF THE JACOBIAN FROM COLUMN 1 TO N-P TS2DTR
C FLAG = 2 : SINGULARITY DETECTED DURING FACTORIZATION TS2DTR
C OF THE JACOBIAN FROM COLUMN N-P+1 TO N TS2DTR
C DXN : DIAGONAL SCALING MATRIX FOR X TS2DTR
C DFN : DIAGONAL SCALING MATRIX FOR F TS2DTR
C TSFVEC : SUBROUTINE TO EVALUATE THE USER'S FUNCTION TS2DTR
C D,XPLSP,ADT,AG,TEMP,VN,VNP,VNS : WORKING VECTORS TS2DTR
C WRK1,CONST1,CONST2,X: WORKING VECTORS TS2DTR
C TS2DTR
C INPUT-OUTPUT PARAMETERS : TS2DTR
C ------------------------ TS2DTR
C TS2DTR
C DLT : INITIAL TRUST RADIUS (= -1.0D0) IF IT IS NOT SUPPLIED TS2DTR
C BY THE USER ON ENTRY AND CURRENT TRUST RADIUS ON EXIT TS2DTR
C TS2DTR
C OUTPUT PARAMETERS : TS2DTR
C ------------------- TS2DTR
C TS2DTR
C FNORM : NORM OF THE FUNCTION VALUE AT CURRENT ITERATE TS2DTR
C XPLS : NEXT ITERATE TS2DTR
C FP : FUNCTION VALUE AT NEXT ITERATE TS2DTR
C FPLS : 0.5 * || FP ||**2 TS2DTR
C RETCD : RETURN CODE FROM SUBROUTINE (SEE SUBROUTINE TSTRUD TS2DTR
C FOR MEANING ) TS2DTR
C TS2DTR
C SUBPROGRAMS CALLED: TS2DTR
C TS2DTR
C LEVEL 1 BLAS ... DAXPY,DCOPY,DDOT,DNRM2,DSCAL TS2DTR
C TENSOLVE ... TSPRMV,TSUTMD,TSJMUV,TSUDQV,TSSMIN,TSRSID, TS2DTR
C TENSOLVE ... TSTRUD TS2DTR
C TS2DTR
C***********************************************************************TS2DTR
TS2DTR
INTEGER I TS2DTR
DOUBLE PRECISION FNORM,RES,ALPH,SUM,RESG,OPTIM TS2DTR
DOUBLE PRECISION SCRES,FPLSP,RRES,RRESG TS2DTR
DOUBLE PRECISION DNRM2,DDOT TS2DTR
DOUBLE PRECISION ZERO,HALF,ONE TS2DTR
LOGICAL DTAKEN TS2DTR
INTRINSIC SQRT TS2DTR
EXTERNAL TSFVEC TS2DTR
DATA ZERO,HALF,ONE/0.0D0,0.50D0,1.0D0/ TS2DTR
TS2DTR
DTAKEN = .FALSE. TS2DTR
RETCD = 4 TS2DTR
IF(DLT.EQ.-ONE) THEN TS2DTR
TS2DTR
c set DLT to length of Cauchy step TS2DTR
TS2DTR
ALPH = DNRM2(N,G,1) TS2DTR
ALPH = ALPH**2 TS2DTR
CALL TSPRMV(VN,G,PIVOT,N,1) TS2DTR
IF(IERR.EQ.0) THEN TS2DTR
CALL TSUTMD(AJA,VN,MAXM,M,N,VNP) TS2DTR
ELSE TS2DTR
CALL TSPRMV(VNS,VN,PBAR,N,1) TS2DTR
CALL TSUTMD(AJA,VNS,MAXM,M+N,N,VNP) TS2DTR
ENDIF TS2DTR
DLT = ALPH*SQRT(ALPH)/DNRM2(N,VNP,1)**2 TS2DTR
IF(DLT.GT.STEPMX) THEN TS2DTR
DLT = STEPMX TS2DTR
ENDIF TS2DTR
ENDIF TS2DTR
TS2DTR
c form an orthonormal basis for the two-dimensional subspace TS2DTR
TS2DTR
CALL DCOPY(N,G,1,GBAR,1) TS2DTR
CALL DSCAL(N,-ONE,GBAR,1) TS2DTR
RES = DNRM2(N,DT,1) TS2DTR
SUM = -DDOT(N,GBAR,1,DT,1)/RES**2 TS2DTR
CALL DAXPY(N,SUM,DT,1,GBAR,1) TS2DTR
RESG = DNRM2(N,GBAR,1) TS2DTR
IF(RESG.GT.ZERO) THEN TS2DTR
RRESG = ONE/RESG TS2DTR
CALL DSCAL(N,RRESG,GBAR,1) TS2DTR
ENDIF TS2DTR
RRES = ONE/RES TS2DTR
CALL DSCAL(N,RRES,DT,1) TS2DTR
TS2DTR
c compute Jacobian times DT TS2DTR
TS2DTR
CALL TSJMUV(ITN,METHOD,DT,CURPOS,PIVOT,PBAR,AJA,SHAT, TS2DTR
+ FLAG,IERR,MAXM,MAXN,M,N,P,D,TEMP,VN,ADT) TS2DTR
TS2DTR
c compute Jacobian times GBAR TS2DTR
TS2DTR
CALL TSJMUV(ITN,METHOD,GBAR,CURPOS,PIVOT,PBAR,AJA,SHAT, TS2DTR
+ FLAG,IERR,MAXM,MAXN,M,N,P,D,TEMP,VNP,AG) TS2DTR
TS2DTR
IF(.NOT. NWTAKE) THEN TS2DTR
TS2DTR
c compute SHAT times VN TS2DTR
TS2DTR
CALL TSUDQV(SHAT,VN,MAXN,N,P,CONST1) TS2DTR
TS2DTR
c compute SHAT times VNP TS2DTR
TS2DTR
CALL TSUDQV(SHAT,VNP,MAXN,N,P,CONST2) TS2DTR
ENDIF TS2DTR
TS2DTR
TS2DTR
70 CONTINUE TS2DTR
TS2DTR
c normalize DT TS2DTR
TS2DTR
IF(RES.LE.DLT) THEN TS2DTR
DTAKEN = .TRUE. TS2DTR
DO 80 I = 1,N TS2DTR
D(I) = DT(I)*RES TS2DTR
80 CONTINUE TS2DTR
DLT = RES TS2DTR
TS2DTR
ELSE TS2DTR
TS2DTR
c find the global minimizer of one-variable function in the TS2DTR
c interval (-dlt, dlt) TS2DTR
TS2DTR
CALL TSSMIN(ANLS,FQ,ADT,AG,CONST1,CONST2,DLT,MAXM,M,N, TS2DTR
+ P,NWTAKE,IERR,EPSM,VN,VNP,VNS,OPTIM) TS2DTR
TS2DTR
c compute the global step TS2DTR
TS2DTR
DO 90 I = 1,N TS2DTR
D(I) = OPTIM*DT(I)+SQRT(DLT**2-OPTIM**2)*GBAR(I) TS2DTR
90 CONTINUE TS2DTR
TS2DTR
ENDIF TS2DTR
TS2DTR
c compute the tensor model residual TS2DTR
TS2DTR
CALL TSRSID(ITN,METHOD,FQ,D,CURPOS,PIVOT,PBAR,AJA,ANLS, TS2DTR
+ SHAT,FLAG,NWTAKE,IERR,MAXM,MAXN,M,N,P,WRK1,VN, TS2DTR
+ VNP,VNS,SCRES) TS2DTR
TS2DTR
c check whether the global step is acceptable TS2DTR
TS2DTR
CALL TSTRUD(M,N,XC,FNORM,G,D,DTAKEN,STEPMX,STEPTL,DLT, TS2DTR
+ MXTAKE,DXN,DFN,TSFVEC,SCRES,RETCD,XPLSP,FPLSP, TS2DTR
+ TEMP,XPLS,FP,FPLS) TS2DTR
TS2DTR
IF(RETCD.GE.2) GO TO 70 TS2DTR
TS2DTR
END TS2DTR
SUBROUTINE TSBSLV(R,NR,M,N,B,Y) TSBSLV
TSBSLV
INTEGER NR,M,N TSBSLV
DOUBLE PRECISION R(NR,N),B(N),Y(N) TSBSLV
TSBSLV
C********************************************************************* TSBSLV
C THIS ROUTINE DOES A BACKWARD SOLVE. TSBSLV
C********************************************************************* TSBSLV
C TSBSLV
C INPUT PARAMETERS : TSBSLV
C ----------------- TSBSLV
C TSBSLV
C R : UPPER TRIANGULAR MATRIX OBTAINED FROM A QR FACTORIZATION TSBSLV
C OF AN M BY N MATRIX A. DIAG(R) IS STORED IN ROW M+2. THIS TSBSLV
C IS THE STORAGE SCHEME USED IN STEWART, G. W., III(1973) TSBSLV
C "INTRODUCTION TO MATRIX COMPUTATION", ACADEMIC PRESS, TSBSLV
C NEW YORK TSBSLV
C NR : LEADING DIMENSION OF MATRIX A TSBSLV
C M : ROW DIMENSION OF MATRIX A TSBSLV
C N : COLUMN DIMENSION OF MATRIX A TSBSLV
C B : RIGHT HAND SIDE TSBSLV
C TSBSLV
C OUTPUT PARAMETERS : TSBSLV
C ----------------- TSBSLV
C TSBSLV
C Y : VECTOR SOLUTION ON EXIT TSBSLV
C TSBSLV
C SUBPROGRAMS CALLED: TSBSLV
C TSBSLV
C LEVEL 1 BLAS ... DAXPY,DLOAD TSBSLV
C TSBSLV
C********************************************************************* TSBSLV
TSBSLV
INTEGER J TSBSLV
DOUBLE PRECISION ZERO TSBSLV
DATA ZERO/0.0D0/ TSBSLV
TSBSLV
c solve R Y = B TSBSLV
TSBSLV
Y(N) = B(N) / R(M+2,N) TSBSLV
TSBSLV
IF(N .GT. 2) THEN TSBSLV
CALL DLOAD(N-1,ZERO,Y,1) TSBSLV
DO 20 J = N-1,2,-1 TSBSLV
CALL DAXPY(J,Y(J+1),R(1,J+1),1,Y,1) TSBSLV
Y(J) = (B(J)-Y(J))/R(M+2,J) TSBSLV
20 CONTINUE TSBSLV
Y(1) = Y(1) + R(1,2) * Y(2) TSBSLV
Y(1) = (B(1) - Y(1)) / R(M+2,1) TSBSLV
ELSEIF(N .EQ. 2) THEN TSBSLV
Y(1) = (B(1) - (R(1,2) * Y(2))) / R(M+2,1) TSBSLV
ENDIF TSBSLV
TSBSLV
RETURN TSBSLV
END TSBSLV
SUBROUTINE TSCHKI(MAXM,MAXN,MAXP,M,N,LUNC,LNEM,LNEN,LIN,GRADTL, TSCHKI
+ STEPTL,FTOL,ITNLIM,JACFLG,METHOD,GLOBAL, TSCHKI
+ STEPMX,DLT,EPSM,MSG,TYPX,TYPF,DXN,DFN, TSCHKI
+ SQRN,TERMCD,IPR) TSCHKI
TSCHKI
INTEGER MAXM,MAXN,MAXP,M,N,LUNC,LNEM,LNEN,LIN,IPR,MSG,JACFLG TSCHKI
INTEGER METHOD,GLOBAL,ITNLIM,SQRN,TERMCD TSCHKI
DOUBLE PRECISION GRADTL,STEPTL,FTOL,STEPMX,DLT,EPSM TSCHKI
DOUBLE PRECISION TYPX(N),TYPF(M),DXN(N),DFN(M) TSCHKI
TSCHKI
C********************************************************************* TSCHKI
C THIS ROUTINE CHECKS INPUT FOR REASONABLENESS. TSCHKI
C********************************************************************* TSCHKI
C TSCHKI
C INPUT PARAMETERS : TSCHKI
C ----------------- TSCHKI
C TSCHKI
C MAXM : LEADING DIMENSION OF WORKSPACE WRKNEM (SEE TSNECI) TSCHKI
C MAXN : LEADING DIMENSION OF WORKSPACE WRKNEN (SEE TSNECI) TSCHKI
C MAXP : LEADING DIMENSION OF WORKSPACE WRKUNC (SEE TSNECI) TSCHKI
C M,N : DIMENSIONS OF PROBLEM TSCHKI
C IPR : DEVICE TO WHICH TO SEND OUTPUT TSCHKI
C TSCHKI
C INPUT-OUTPUT PARAMETERS : TSCHKI
C ------------------------ TSCHKI
C TSCHKI
C GRADTL : TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE TSCHKI
C ENOUGH TO ZERO TO TERMINATE ALGORITHM TSCHKI
C STEPTL : TOLERANCE AT WHICH SUCCESSIVE ITERATES TSCHKI
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM TSCHKI
C FTOL : TOLERANCE AT WHICH FUNCTION VALUE CONSIDERED TSCHKI
C CLOSE ENOUGH TO ZERO TSCHKI
C ITNLIM : MAXIMUM NUMBER OF ALLOWABLE ITERATIONS TSCHKI
C STEPMX : MAXIMUM STEP ALLOWED IN TRUST REGION TSCHKI
C DLT : TRUST RADIUS TSCHKI
C JACFLG : JACOBIAN FLAG WITH THE FOLLOWING MEANINGS : TSCHKI
C JACFLG = 1 : ANALYTIC JACOBIAN SUPPLIED TSCHKI
C JACFLG = 0 : ANALYTIC JACOBIAN NOT SUPPLIED TSCHKI
C METHOD : METHOD TO USE TSCHKI
C METHOD = 0 : STANDARD METHOD IS USED TSCHKI
C METHOD = 1 : TENSOR METHOD IS USED TSCHKI
C GLOBAL : GLOBAL STRATEGY USED TSCHKI
C GLOBAL = 0 : LINE SEARCH USED TSCHKI
C GLOBAL = 1 : 2-DIMENSIONAL TRUST REGION USED TSCHKI
C TYPX : TYPICAL SIZE FOR EACH COMPONENT OF X TSCHKI
C TYPF : TYPICAL SIZE FOR EACH COMPONENT OF F TSCHKI
C MSG : MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT TSCHKI
C TSCHKI
C OUTPUT PARAMETERS : TSCHKI
C ------------------- TSCHKI
C TSCHKI
C TERMCD: TERMINATION CODE TSCHKI
C DXN : DIAGONAL SCALING MATRIX FOR X TSCHKI
C DFN : DIAGONAL SCALING MATRIX FOR F TSCHKI
C SQRN : MAXIMUM COLUMN DIMENSION OF S AND FV TSCHKI
C TSCHKI
C SUBPROGRAMS CALLED: TSCHKI
C TSCHKI
C UNCMIN ... DPMEPS TSCHKI
C TSCHKI
C********************************************************************* TSCHKI
TSCHKI
INTEGER I,LEN TSCHKI
DOUBLE PRECISION DPMEPS,ZERO,ONE,TWO,THREE,THOUS,TEMP TSCHKI
INTRINSIC MOD,NINT,SQRT TSCHKI
DATA ZERO,ONE,TWO,THREE,THOUS/0.0D0,1.0D0,2.0D0,3.0D0,1000.0D0/ TSCHKI
TSCHKI
c check that parameters only take on acceptable values TSCHKI
c if not, set them to default values TSCHKI
TSCHKI
c set TERMCD to zero in case we abort prematuraly TSCHKI
TSCHKI
TERMCD = 0 TSCHKI
TSCHKI
c compute machine precision TSCHKI
TSCHKI
EPSM = DPMEPS() TSCHKI
TSCHKI
c check dimensions of the problem TSCHKI
TSCHKI
IF(M.LE.0) THEN TSCHKI
WRITE(IPR,601) M TSCHKI
MSG = -1 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
IF(N.LE.0) THEN TSCHKI
WRITE(IPR,602) N TSCHKI
MSG = -2 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
c check leading dimensions of the problem TSCHKI
TSCHKI
LEN = M+N+2 TSCHKI
IF(MAXM .LT. LEN) THEN TSCHKI
WRITE(IPR,603) MAXM,LEN TSCHKI
MSG = -3 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
LEN = N+2 TSCHKI
IF(MAXN .LT. LEN) THEN TSCHKI
WRITE(IPR,604) MAXN,LEN TSCHKI
MSG = -4 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
TEMP = SQRT(DBLE(N)) TSCHKI
SQRN = NINT(TEMP) TSCHKI
TSCHKI
IF(MAXP .LT. SQRN) THEN TSCHKI
WRITE(IPR,605) MAXP,SQRN TSCHKI
MSG = -5 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
c check column dimensions of workspace arrays TSCHKI
TSCHKI
LEN = 2*SQRN+4 TSCHKI
IF(LUNC.LT.LEN) THEN TSCHKI
WRITE(IPR,606) LUNC,LEN TSCHKI
MSG = -6 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
LEN = N+2*SQRN+11 TSCHKI
IF(LNEM.LT.LEN) THEN TSCHKI
WRITE(IPR,607) LNEM,LEN TSCHKI
MSG = -7 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
LEN = 2*SQRN+9 TSCHKI
IF(LNEN.LT.LEN) THEN TSCHKI
WRITE(IPR,608) LNEN,LEN TSCHKI
MSG = -8 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
IF(LIN.LT.3) THEN TSCHKI
WRITE(IPR,609) LIN TSCHKI
MSG = -9 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
c check JACFLG, METHOD, and GLOBAL TSCHKI
TSCHKI
IF(JACFLG.NE.1) JACFLG = 0 TSCHKI
TSCHKI
IF(METHOD.NE.0 .AND. METHOD.NE.1) METHOD = 1 TSCHKI
TSCHKI
IF(GLOBAL.NE.0 .AND. GLOBAL.NE.1) GLOBAL = 0 TSCHKI
TSCHKI
IF(MOD(MSG/2,2).EQ.1 .AND. JACFLG.EQ.0) THEN TSCHKI
WRITE(IPR,610) MSG,JACFLG TSCHKI
MSG = -10 TSCHKI
RETURN TSCHKI
ENDIF TSCHKI
TSCHKI
c check scale matrices TSCHKI
TSCHKI
DO 10 I = 1,N TSCHKI
IF(TYPX(I).EQ.ZERO) TYPX(I) = ONE TSCHKI
IF(TYPX(I).LT.ZERO) TYPX(I) = -TYPX(I) TSCHKI
DXN(I) = ONE/TYPX(I) TSCHKI
10 CONTINUE TSCHKI
TSCHKI
DO 20 I = 1,M TSCHKI
IF(TYPF(I).EQ.ZERO) TYPF(I) = ONE TSCHKI
IF(TYPF(I).LT.ZERO) TYPF(I) = -TYPF(I) TSCHKI
DFN(I) = ONE/TYPF(I) TSCHKI
20 CONTINUE TSCHKI
TSCHKI
c check gradient, step, and function tolerances TSCHKI
TSCHKI
TEMP = ONE/THREE TSCHKI
IF(GRADTL.LT.ZERO) THEN TSCHKI
GRADTL = EPSM**TEMP TSCHKI
ENDIF TSCHKI
TSCHKI
IF(STEPTL.LT.ZERO) THEN TSCHKI
STEPTL = EPSM**(TWO*TEMP) TSCHKI
ENDIF TSCHKI
TSCHKI
IF(FTOL.LT.ZERO) THEN TSCHKI
FTOL = EPSM**(TWO*TEMP) TSCHKI
ENDIF TSCHKI
TSCHKI
c check iteration limit TSCHKI
TSCHKI
IF(ITNLIM.LE.0) THEN TSCHKI
ITNLIM = 150 TSCHKI
ENDIF TSCHKI
TSCHKI
c check STEPMX and DLT TSCHKI
TSCHKI
IF(STEPMX.LT.ZERO) STEPMX = THOUS TSCHKI
TSCHKI
IF(DLT.LE.ZERO) THEN TSCHKI
DLT = -ONE TSCHKI
IF(DLT.GT.STEPMX) DLT = STEPMX TSCHKI
ENDIF TSCHKI
TSCHKI
601 FORMAT(' TSCHKI ILLEGAL DIMENSION M =',I5) TSCHKI
602 FORMAT(' TSCHKI ILLEGAL DIMENSION N =',I5) TSCHKI
603 FORMAT(' TSCHKI ILLEGAL DIMENSION MAXM =',I5,' < M+N+2 =',I5)TSCHKI
604 FORMAT(' TSCHKI ILLEGAL DIMENSION MAXN =',I5,' < N+2 =',I5) TSCHKI
605 FORMAT(' TSCHKI ILLEGAL DIMENSION MAXP =',I5,' <', TSCHKI
+ ' NINT(SQRT (N)) =',I5) TSCHKI
606 FORMAT(' TSCHKI ILLEGAL DIMENSION LUNC =',I5,' <', TSCHKI
+ ' 2*NINT(SQRT (N))+4 =',I5) TSCHKI
607 FORMAT(' TSCHKI ILLEGAL DIMENSION LNEM =',I5,' <', TSCHKI
+ ' N+2*NINT(SQRT (N))+11 =',I5) TSCHKI
608 FORMAT(' TSCHKI ILLEGAL DIMENSION LNEN =',I5,' <', TSCHKI
+ ' 2*NINT(SQRT (N))+9 =',I5) TSCHKI
609 FORMAT(' TSCHKI ILLEGAL DIMENSION LIN =',I5,' < 3') TSCHKI
610 FORMAT(' TSCHKI USER REQUESTS THAT ANALYTIC JACOBIAN BE', TSCHKI
+ ' ACCEPTED AS PROPERLY CODED (MSG =',I5,')'/ TSCHKI
+ ' TSCHKI BUT ANALYTIC JACOBIAN NOT SUPPLIED', TSCHKI
+ ' (JACFLG =',I5,')') TSCHKI
END TSCHKI
SUBROUTINE TSCHKJ(AJANAL,XC,FC,NR,M,N,EPSM,DFN,DXN, TSCHKJ
+ TYPX,IPR,FHAT,WRK1,TSFVEC,MSG) TSCHKJ
TSCHKJ
INTEGER NR,M,N,IPR,MSG TSCHKJ
DOUBLE PRECISION AJANAL(NR,N),XC(N),FC(M) TSCHKJ
DOUBLE PRECISION EPSM,DFN(M),DXN(N),TYPX(N) TSCHKJ
DOUBLE PRECISION FHAT(M),WRK1(M) TSCHKJ
EXTERNAL TSFVEC TSCHKJ
TSCHKJ
C********************************************************************* TSCHKJ
C THIS ROUTINE CHECKS THE ANALYTIC JACOBIAN AGAINST ITS FINITE TSCHKJ
C DIFFERENCE APPROXIMATION. TSCHKJ
C********************************************************************* TSCHKJ
C TSCHKJ
C INPUT PARAMETERS : TSCHKJ
C ----------------- TSCHKJ
C TSCHKJ
C AJANAL : ANALYTIC JACOBIAN AT XC TSCHKJ
C XC : CURRENT ITERATE TSCHKJ
C FC : FUNCTION VALUE AT XC TSCHKJ
C NR : LEADING DIMENSION OF AJANAL TSCHKJ
C M,N : DIMENSIONS OF PROBLEM TSCHKJ
C EPSM : MACHINE PRECISION TSCHKJ
C DFN : DIAGONAL SCALING MATRIX FOR F TSCHKJ
C DXN : DIAGONAL SCALING MATRIX FOR X TSCHKJ
C TYPX : TYPICAL SIZE FOR EACH COMPONENT OF X TSCHKJ
C IPR : DEVICE TO WHICH TO SEND OUTPUT TSCHKJ
C FHAT,WRK1 : WORKSPACE TSCHKJ
C TSFVEC : SUBROUTINE TO EVALUATE THE USER'S FUNCTION TSCHKJ
C TSCHKJ
C INPUT-OUTPUT PARAMETERS : TSCHKJ
C ------------------------ TSCHKJ
C TSCHKJ
C MSG : MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT TSCHKJ
C TSCHKJ
C SUBPROGRAMS CALLED: TSCHKJ
C TSCHKJ
C LEVEL 1 BLAS ... IDAMAX TSCHKJ
C TENSOLVE ... TSUNSX,TSUNSF,TSSCLX,TSSCLF TSCHKJ
C USER ... TSFVEC TSCHKJ
C TSCHKJ
C********************************************************************* TSCHKJ
TSCHKJ
INTEGER I,J TSCHKJ
DOUBLE PRECISION NDIGIT,RNOISE,SQRNS,STEPSZ,XTMPJ,DINF,RSTPSZ TSCHKJ
DOUBLE PRECISION TOL,QUART,ONE,TEN TSCHKJ
INTEGER IDAMAX TSCHKJ
INTRINSIC ABS,MAX,SQRT TSCHKJ
DATA QUART,ONE,TEN/0.250D0,1.0D0,10.0D0/ TSCHKJ
TSCHKJ
c unscale XC and FC TSCHKJ
TSCHKJ
CALL TSUNSX(XC,DXN,N) TSCHKJ
CALL TSUNSF(FC,DFN,M) TSCHKJ
TSCHKJ
c compute the finite difference Jacobian and check it against TSCHKJ
c the analytic one TSCHKJ
TSCHKJ
NDIGIT = -LOG10(EPSM) TSCHKJ
RNOISE = MAX(TEN**(-NDIGIT),EPSM) TSCHKJ
SQRNS = SQRT(RNOISE) TSCHKJ
TOL = EPSM**QUART TSCHKJ
DO 40 J = 1,N TSCHKJ
STEPSZ = SQRNS*MAX(ABS(XC(J)),ONE) TSCHKJ
XTMPJ = XC(J) TSCHKJ
XC(J) = XTMPJ+STEPSZ TSCHKJ
CALL TSFVEC(XC,FHAT,M,N) TSCHKJ
XC(J) = XTMPJ TSCHKJ
RSTPSZ = ONE/STEPSZ TSCHKJ
DO 10 I = 1,M TSCHKJ
WRK1(I) = (FHAT(I)-FC(I))*RSTPSZ TSCHKJ
10 CONTINUE TSCHKJ
DO 20 I = 1,M TSCHKJ
WRK1(I) = WRK1(I)*DFN(I)*TYPX(J) TSCHKJ
20 CONTINUE TSCHKJ
DINF = ABS(WRK1(IDAMAX(M,WRK1,1))) TSCHKJ
DO 30 I = 1,M TSCHKJ
IF(ABS(AJANAL(I,J)-WRK1(I)).GT.TOL*DINF) THEN TSCHKJ
PRINT *, TOL*DINF TSCHKJ
PRINT *, AJANAL(I,J) TSCHKJ
PRINT *, WRK1(I) TSCHKJ
PRINT *, ABS(AJANAL(I,J)-WRK1(I)) TSCHKJ
PRINT *, I TSCHKJ
PRINT *, J TSCHKJ
WRITE(IPR,50) TSCHKJ
MSG = -11 TSCHKJ
RETURN TSCHKJ
ENDIF TSCHKJ
30 CONTINUE TSCHKJ
40 CONTINUE TSCHKJ
TSCHKJ
c scale back XC and FC TSCHKJ
TSCHKJ
CALL TSSCLX(XC,DXN,N) TSCHKJ
CALL TSSCLF(FC,DFN,M) TSCHKJ
TSCHKJ
50 FORMAT(/,' TSCHKJ PROBABLE ERROR IN CODING OF ANALYTIC', TSCHKJ
+ ' JACOBIAN') TSCHKJ
TSCHKJ
RETURN TSCHKJ
END TSCHKJ
SUBROUTINE TSCPMU(R,NR,N,EPSM,MU) TSCPMU
TSCPMU
INTEGER NR,N TSCPMU
DOUBLE PRECISION R(NR,N),EPSM,MU TSCPMU
TSCPMU
C********************************************************************* TSCPMU
C THIS ROUTINE COMPUTES A SMALL PERTURBATION MU. MU IS USED IN THE TSCPMU
C COMPUTATION OF THE LEVENBERG-MARQUARDT STEP. TSCPMU
C********************************************************************* TSCPMU
C TSCPMU
C INPUT PARAMETERS : TSCPMU
C ----------------- TSCPMU
C TSCPMU
C R : UPPER TRIANGULAR MATRIX TSCPMU
C NR : LEADING DIMENSION OF R TSCPMU
C N : COLUMN DIMENSION OF R TSCPMU
C EPSM : MACHINE PRECISION TSCPMU
C TSCPMU
C OUTPUT PARAMETERS : TSCPMU
C ------------------ TSCPMU
C TSCPMU
C MU = SQRT(L1 NORM OF R * INFINITY NORM OF R * N * EPSM * 100) TSCPMU
C TSCPMU
C SUBPROGRAMS CALLED: TSCPMU
C TSCPMU
C LEVEL 1 BLAS ... DASUM TSCPMU
C TSCPMU
C********************************************************************* TSCPMU
TSCPMU
INTEGER I,J TSCPMU
DOUBLE PRECISION AIFNRM,SUM,AL1NRM,ZERO,HUND TSCPMU
DOUBLE PRECISION DASUM TSCPMU
INTRINSIC ABS,MAX,SQRT TSCPMU
DATA ZERO,HUND/0.0D0,100.0D0/ TSCPMU
TSCPMU
c compute the infinity norm of R TSCPMU
TSCPMU
AIFNRM = ZERO TSCPMU
DO 20 I = 1,N TSCPMU
SUM = ZERO TSCPMU
DO 10 J = I,N TSCPMU
SUM = SUM+ABS(R(I,J)) TSCPMU
10 CONTINUE TSCPMU
AIFNRM = MAX(AIFNRM,SUM) TSCPMU
20 CONTINUE TSCPMU
TSCPMU
c compute the l1 norm of R TSCPMU
TSCPMU
AL1NRM = ZERO TSCPMU
DO 40 J = 1,N TSCPMU
SUM = DASUM(J,R(1,J),1) TSCPMU
AL1NRM = MAX(AL1NRM,SUM) TSCPMU
40 CONTINUE TSCPMU
TSCPMU
c compute MU TSCPMU
TSCPMU
MU = SQRT(AIFNRM*AL1NRM*N*EPSM*HUND) TSCPMU
TSCPMU
RETURN TSCPMU
END TSCPMU
SUBROUTINE TSCPSS(S,MAXM,MAXN,M,N,P,METHOD,GLOBAL,EPSM,FCQ, TSCPSS
+ Y,W,FQT,AL2NRM,QHAT,ANLS,DN,FQQ,PTILDA, TSCPSS
+ CURPOS,PBAR,ZERO1,IERR,RESNEW,FLAG) TSCPSS
TSCPSS
INTEGER MAXM,MAXN,M,N,P,FLAG,ZERO1,GLOBAL,IERR TSCPSS
DOUBLE PRECISION EPSM,RESNEW TSCPSS
INTEGER METHOD,PTILDA(N),CURPOS(N),PBAR(N) TSCPSS
DOUBLE PRECISION S(MAXN,P),FCQ(M) TSCPSS
DOUBLE PRECISION Y(N),W(M),FQT(M),AL2NRM(N) TSCPSS
DOUBLE PRECISION QHAT(MAXM,N),ANLS(MAXM,P) TSCPSS
DOUBLE PRECISION DN(N),FQQ(M) TSCPSS
TSCPSS
C********************************************************************** TSCPSS
C THIS ROUTINE COMPUTES THE STANDARD STEP. NOTE THAT AT THIS STAGE TSCPSS
C THE JACOBIAN MATRIX (QHAT) HAS ALREADY BEEN FACTORED FROM COLUMNS 1 TSCPSS
C TO N-P DURING THE TENSOR STEP COMPUTATION. THIS ROUTINE FACTORS TSCPSS
C THE MATRIX QHAT FROM COLUMN N-P+1 TO N, THEREBY OBTAINING A QR TSCPSS
C FACTORIZATION OF THE FULL MATRIX QHAT, THEN COMPUTES THE STANDARD TSCPSS
C STEP BY PREMULTIPLYING THE RIGH-HAND SIDE FCQ BY AN ORTHOGONAL TSCPSS
C MATRIX AND BY PERFORMING A BACKWARD SOLVE. TSCPSS
C********************************************************************** TSCPSS
C TSCPSS
C TSCPSS
C INPUT PARAMETERS : TSCPSS
C ----------------- TSCPSS
C TSCPSS
C S : FACTORED MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSCPSS
C (OBTAINED FROM TSQLFC SUBROUTINE) TSCPSS
C MAXM : LEADING DIMENSION OF QHAT AND ANLS TSCPSS
C MAXN : LEADING DIMENSION OF S TSCPSS
C M,N : DIMENSIONS OF PROBLEM TSCPSS
C P : COLUMN DIMENSION OF MATRIX S TSCPSS
C METHOD : METHOD USED : TSCPSS
C METHOD = 0 : STANDARD METHOD IS USED TSCPSS
C METHOD = 1 : TENSOR METHOD IS USED TSCPSS
C GLOBAL : GLOBAL STRATEGY USED TSCPSS
C GLOBAL = 0 : LINE SEARCH IS USED TSCPSS
C GLOBAL = 1 : 2-DIMENSIONAL TRUST REGION IS USED TSCPSS
C EPSM : MACHINE PRECISION TSCPSS
C FCQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY AN TSCPSS
C ORTHOGONAL MATRIX TSCPSS
C CURPOS : PIVOT VECTOR FOR THE FACTORIZATION OF QHAT FROM COLUMN TSCPSS
C 1 TO N-P TSCPSS
C Y,W,FQT,AL2NRM : WORKING VECTORS TSCPSS
C TSCPSS
C TSCPSS
C INPUT-OUTPUT PARAMETERS : TSCPSS
C ------------------------ TSCPSS
C TSCPSS
C QHAT : FACTORED MATRIX FROM COLUMN 1 TO N-P TSCPSS
C ON ENTRY AND FACTORED MATRIX FROM 1 TO N ON EXIT TSCPSS
C ANLS : TENSOR TERM MATRIX ON ENTRY AND ANLS MULTIPLIED BY TSCPSS
C ORTHOGONAL MATRICES ON EXIT (THIS IS PERFORMED IN THE TSCPSS
C CASE WHERE THE GLOBAL STRATEGY USED IS THE 2-DIMENSIONALTSCPSS
C TRUST REGION) TSCPSS
C TSCPSS
C OUTPUT PARAMETERS : TSCPSS
C ------------------- TSCPSS
C TSCPSS
C DN : STANDARD STEP TSCPSS
C FQQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSCPSS
C ORTHOGONAL MATRICES (THIS IS USED IN THE CASE WHERE TSCPSS
C THE GLOBAL STRATEGY USED IS THE 2-DIMENSIONAL TSCPSS
C TRUST REGION) TSCPSS
C PTILDA: PIVOT VECTOR FOR THE FACTORIZATION OF THE TSCPSS
C MATRIX QHAT FROM N-P+1 TO N TSCPSS
C PBAR : PIVOT VECTOR FOR THE FACTORIZATION OF THE TSCPSS
C TRANSFORMED MATRIX QHAT FROM 1 TO N TSCPSS
C IN CASE OF SINGULARITY TSCPSS
C ZERO1 : FIRST ZERO COLUMN OF MATRIX QHAT IN CASE OF SINGULARITY TSCPSS
C IERR : RETURNED CODE WITH THE FOLLOWING MEANING : TSCPSS
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSCPSS
C IERR = 0 : OTHERWISE TSCPSS
C RESNEW: RESIDUAL OF THE STANDARD MODEL TSCPSS
C FLAG : RETURNED CODE WITH THE FOLLOWING MEANINGS : TSCPSS
C FLAG = 0 : NO SINGULARITY DETECTED TSCPSS
C FLAG = 1 : SINGULARITY DETECTED DURING QR FACTORIZATION TSCPSS
C OF QHAT FROM COLUMN 1 TO N-P TSCPSS
C FLAG = 2 : SINGULARITY DETECTED DURING QR FACTORIZATION TSCPSS
C OF QHAT FROM COLUMN N-P+1 TO N TSCPSS
C TSCPSS
C SUBPROGRAMS CALLED: TSCPSS
C TSCPSS
C LEVEL 1 BLAS ... DCOPY,DLOAD,DSCAL TSCPSS
C TENSOLVE ... TSQRFC,TSQMUV,TSQMTS,TSSMRD,TSBSLV,TSPRMV TSCPSS
C TENSOLVE ... TSQMLV,TSCPMU TSCPSS
C TSCPSS
C **********************************************************************TSCPSS
TSCPSS
INTEGER ZEROTM,I,J TSCPSS
DOUBLE PRECISION MU,ZERO,ONE TSCPSS
DATA ZERO,ONE/0.0D0,1.0D0/ TSCPSS
TSCPSS
FLAG = 0 TSCPSS
TSCPSS
c initialization TSCPSS
TSCPSS
CALL DLOAD (M+N,ZERO,FQQ,1) TSCPSS
TSCPSS
CALL DCOPY(M,FCQ,1,W,1) TSCPSS
CALL DSCAL(M,-ONE,W,1) TSCPSS
TSCPSS
c if the Jacobian is singular then compute the Levenberg-Marquardt TSCPSS
c step (label 20) TSCPSS
TSCPSS
IF(IERR.EQ.1) THEN TSCPSS
FLAG = 1 TSCPSS
GO TO 20 TSCPSS
ENDIF TSCPSS
TSCPSS
c factor the matrix QHAT from column n-p+1 to n TSCPSS
TSCPSS
CALL TSQRFC(QHAT,MAXM,N,M,N-P+1,N,IERR,EPSM,AL2NRM,PTILDA,ZERO1) TSCPSS
TSCPSS
IF((M.EQ.N).AND.(IERR.EQ.0)) THEN TSCPSS
ZEROTM = ZERO1-1 TSCPSS
ELSE TSCPSS
ZEROTM = ZERO1 TSCPSS
ENDIF TSCPSS
TSCPSS
c premultiply W by the orthogonal matrix resulting from the QR TSCPSS
c factorization of QHAT TSCPSS
TSCPSS
CALL TSQMUV(QHAT,W,FQQ,MAXM,M,N-P+1,ZEROTM,.FALSE.) TSCPSS
TSCPSS
IF(METHOD.EQ.1 .AND. GLOBAL.EQ.1) THEN TSCPSS
TSCPSS
c premultiply ANLS by the orthogonal matrix resulting from the QR TSCPSS
c factorization of QHAT TSCPSS
TSCPSS
CALL TSQMTS(ANLS,QHAT,MAXM,M,N,M,P,N-P+1,FCQ,ZEROTM) TSCPSS
ENDIF TSCPSS
TSCPSS
IF(IERR.EQ.1) THEN TSCPSS
FLAG = 2 TSCPSS
GO TO 20 TSCPSS
ENDIF TSCPSS
TSCPSS
c computate the residual of the standard model TSCPSS
TSCPSS
CALL TSSMRD(FQQ,RESNEW,DN,MU,IERR,M,N) TSCPSS
TSCPSS
c if QHAT is nonsingular perform a backward solve to obtain Y TSCPSS
TSCPSS
CALL TSBSLV(QHAT,MAXM,M,N,FQQ,Y) TSCPSS
TSCPSS
c pivot Y TSCPSS
TSCPSS
CALL TSPRMV(DN,Y,PTILDA,N,0) TSCPSS
TSCPSS
IF(N .NE. 1) THEN TSCPSS
TSCPSS
CALL TSPRMV(Y,DN,CURPOS,N,0) TSCPSS
TSCPSS
c premultiply Y by the orthogonal matrix resulting from the QL TSCPSS
c factorization of S TSCPSS
TSCPSS
CALL TSQMLV(MAXN,N,P,S,Y,DN,.TRUE.) TSCPSS
TSCPSS
ENDIF TSCPSS
TSCPSS
IF(GLOBAL.EQ.1) THEN TSCPSS
IERR = 0 TSCPSS
CALL DSCAL(M,-ONE,FQQ,1) TSCPSS
ENDIF TSCPSS
TSCPSS
RETURN TSCPSS
TSCPSS
20 CONTINUE TSCPSS
TSCPSS
c @ SINGULAR CASE @ TSCPSS
TSCPSS
c solve ( QHAT-trans QHAT + MU I ) DN = -QHAT-trans W TSCPSS
TSCPSS
c put the diagonal elements stored in row m+2 of QHAT into their TSCPSS
c propre positions and zero out the unwanted portions of QHAT TSCPSS
TSCPSS
DO 30 J = 1, ZERO1-1 TSCPSS
QHAT(J,J) = QHAT(M+2,J) TSCPSS
CALL DLOAD (M+N-J,ZERO,QHAT(J+1,J),1) TSCPSS
30 CONTINUE TSCPSS
TSCPSS
DO 40 J = ZERO1, N TSCPSS
CALL DLOAD (M+N-ZERO1+1,ZERO,QHAT(ZERO1,J),1) TSCPSS
40 CONTINUE TSCPSS
TSCPSS
c compute a small perturbation MU TSCPSS
TSCPSS
CALL TSCPMU(QHAT,MAXM,N,EPSM,MU) TSCPSS
TSCPSS
c form the augmented matrix QHAT by adding an (n*n) diag(MU) in TSCPSS
c the bottom TSCPSS
TSCPSS
DO 50 I = M+1,M+N TSCPSS
QHAT(I,I-M) = MU TSCPSS
50 CONTINUE TSCPSS
TSCPSS
c factor the transformed matrix QHAT from 1 to n TSCPSS
TSCPSS
CALL TSQRFC(QHAT,MAXM,N,M+N,1,N,IERR,EPSM,AL2NRM,PBAR,ZERO1) TSCPSS
TSCPSS
IF(METHOD.EQ.1 .AND. GLOBAL.EQ.1) THEN TSCPSS
TSCPSS
c premultiply ANLS by the orthogonal matrix resulting from the QR TSCPSS
c factorization of QHAT TSCPSS
TSCPSS
CALL TSQMTS(ANLS,QHAT,MAXM,M+N,N,M,P,1,FCQ,ZERO1) TSCPSS
ENDIF TSCPSS
TSCPSS
c compute the Levenberg-Marquardt step and the residual of the TSCPSS
c standard model TSCPSS
TSCPSS
IF(FLAG.EQ.1) THEN TSCPSS
CALL TSQMUV(QHAT,W,FQQ,MAXM,M+N,1,N+1,.FALSE.) TSCPSS
CALL TSBSLV(QHAT,MAXM,M+N,N,FQQ,Y) TSCPSS
CALL TSPRMV(DN,Y,PBAR,N,0) TSCPSS
CALL TSPRMV(Y,DN,CURPOS,N,0) TSCPSS
CALL TSQMLV(MAXN,N,P,S,Y,DN,.TRUE.) TSCPSS
CALL TSSMRD(FQQ,RESNEW,DN,MU,IERR,M,N) TSCPSS
IF(GLOBAL.EQ.1) THEN TSCPSS
IERR = 1 TSCPSS
CALL DSCAL(M+N,-ONE,FQQ,1) TSCPSS
ENDIF TSCPSS
RETURN TSCPSS
ELSE TSCPSS
CALL TSQMUV(QHAT,FQQ,FQT,MAXM,M+N,1,N+1,.FALSE.) TSCPSS
CALL TSBSLV(QHAT,MAXM,M+N,N,FQT,DN) TSCPSS
CALL TSPRMV(Y,DN,PBAR,N,0) TSCPSS
CALL TSPRMV(DN,Y,PTILDA,N,0) TSCPSS
CALL TSPRMV(Y,DN,CURPOS,N,0) TSCPSS
CALL TSQMLV(MAXN,N,P,S,Y,DN,.TRUE.) TSCPSS
CALL TSSMRD(FQT,RESNEW,DN,MU,IERR,M,N) TSCPSS
IF(GLOBAL.EQ.1) THEN TSCPSS
IERR = 1 TSCPSS
CALL DCOPY(M+N,FQT,1,FQQ,1) TSCPSS
CALL DSCAL(M+N,-ONE,FQQ,1) TSCPSS
ENDIF TSCPSS
ENDIF TSCPSS
TSCPSS
END TSCPSS
SUBROUTINE TSD1SV(AJA,S,ANLS,FN,X,MAXM,MAXN,M,N,P,Q,EPSM, TSD1SV
+ WRK1,WRK2,WRK3,PIVOT,D1) TSD1SV
TSD1SV
INTEGER MAXM,MAXN,M,N,P,Q TSD1SV
INTEGER PIVOT(N) TSD1SV
DOUBLE PRECISION EPSM TSD1SV
DOUBLE PRECISION AJA(MAXM,N),S(MAXN,P),ANLS(MAXM,P),FN(M),X(P) TSD1SV
DOUBLE PRECISION WRK1(N),WRK2(N),WRK3(N),D1(N) TSD1SV
TSD1SV
C********************************************************************* TSD1SV
C THIS ROUTINE SOLVES THE FIRST N-Q LINEAR EQUATIONS IN N-P UNKNOWNS TSD1SV
C OF THE TENSOR MODEL. TSD1SV
C********************************************************************* TSD1SV
C TSD1SV
C INPUT PARAMETERS : TSD1SV
C ---------------- TSD1SV
C TSD1SV
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSD1SV
C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSD1SV
C ANLS: TENSOR TERM MATRIX AT CURRENT ITERATE TSD1SV
C FN : FUNCTION VALUE AT CURRENT ITERATE TSD1SV
C X : SOLUTION OF THE LOWER M-N+Q QUADRATIC EQUATIONS IN P TSD1SV
C UNKNOWNS OF THE TENSOR MODEL TSD1SV
C MAXM: LEADING DIMENSION OF AJA AND ANLS TSD1SV
C MAXN: LEADING DIMENSION OF S TSD1SV
C M,N : DIMENSIONS OF PROBLEM TSD1SV
C P : COLUMN DIMENSION OF S AND ANLS TSD1SV
C Q : NUMERICAL RANK OF JACOBIAN : TSD1SV
C Q > P : JACOBIAN IS SINGULAR TSD1SV
C Q = P : OTHERWISE TSD1SV
C EPSM: MACHINE PRECISION TSD1SV
C WRK1,WRK2,WRK3 : WORKSPACE TSD1SV
C TSD1SV
C TSD1SV
C OUTPUT PARAMETERS : TSD1SV
C ------------------ TSD1SV
CTSD1SV
C PIVOT : PIVOT VECTOR TSD1SV
C D1 : SOLUTION OF THE N-Q LINEAR EQUATIONS IN N-P UNKNOWNS OF TSD1SV
C THE TENSOR MODEL TSD1SV
C TSD1SV
C SUBPROGRAMS CALLED: TSD1SV
C TSD1SV
C LEVEL 1 BLAS ... DCOPY,DLOAD TSD1SV
C LEVEL 2 BLAS ... DGEMV TSD1SV
C TENSOLVE ... TSSTMX,TSBSLV,TSQRFC,TSPRMV,TSFSLV,TSQMUV TSD1SV
C TSD1SV
C********************************************************************* TSD1SV
TSD1SV
INTEGER ZERO1,I,J,IERR,ICOL TSD1SV
DOUBLE PRECISION EPSM1,ZERO,HALF,ALPHA,ONE TSD1SV
DATA ZERO,ALPHA,HALF,ONE/0.0D0,1.0D-4,0.50D0,1.0D0/ TSD1SV
TSD1SV
c compute the top right (n-q) x p submatrix of AJA times X TSD1SV
TSD1SV
CALL DGEMV('N',N-Q,P,ONE,AJA(1,N-P+1),MAXM, TSD1SV
+ X,1,ZERO,D1,1) TSD1SV
TSD1SV
c compute S-trans times X TSD1SV
TSD1SV
CALL TSSTMX(S,X,MAXN,N,P,WRK3,WRK2) TSD1SV
TSD1SV
c compute 0.5 * (S-trans times X)**2 TSD1SV
TSD1SV
DO 10 I = 1, P TSD1SV
WRK1(I) = HALF * WRK2(I)**2 TSD1SV
10 CONTINUE TSD1SV
TSD1SV
c compute 0.5 * (top (n-q) x p submatrix of ANLS) * TSD1SV
c (S-trans times X)**2 TSD1SV
TSD1SV
CALL DGEMV('N',N-Q,P,ONE,ANLS(1,1),MAXM,WRK1,1,ZERO,WRK2,1) TSD1SV
TSD1SV
DO 20 I = 1,N-Q TSD1SV
WRK1(I) = -FN(I)-D1(I)-WRK2(I) TSD1SV
20 CONTINUE TSD1SV
TSD1SV
c if the Jacobian is nonsingular then solve for the first TSD1SV
c n-p components of the tensor step and return TSD1SV
TSD1SV
IF(P.EQ.Q) THEN TSD1SV
CALL TSBSLV(AJA,MAXM,M,N-P,WRK1,D1) TSD1SV
RETURN TSD1SV
ENDIF TSD1SV
TSD1SV
CALL DLOAD(Q-P,ZERO,WRK2(N-Q+1),1) TSD1SV
TSD1SV
c copy top left (n-q) x (n-p) submatrix of AJA into bottom of AJA TSD1SV
TSD1SV
DO 30 J = 1,N-P TSD1SV
CALL DCOPY(N-Q,AJA(1,J),1,AJA(M+3,J),1) TSD1SV
30 CONTINUE TSD1SV
TSD1SV
c copy the transpose of the top left (n-q) x (n-p) submatrix of AJA TSD1SV
c into top of AJA TSD1SV
TSD1SV
DO 50 J = 1,N-Q TSD1SV
AJA(J,J) = AJA(M+2,J) TSD1SV
DO 40 I = J+1,N-P TSD1SV
AJA(I,J) = AJA(J,I) TSD1SV
40 CONTINUE TSD1SV
50 CONTINUE TSD1SV
TSD1SV
c zero out the upper triangular (n-q) x (n-q) triangular part of TSD1SV
c the transpose of the top left (n-q) x (n-p) submatrix of AJA TSD1SV
TSD1SV
DO 60 J = 1,N-Q TSD1SV
CALL DLOAD(J-1,ZERO,AJA(1,J),1) TSD1SV
60 CONTINUE TSD1SV
TSD1SV
c factor the transpose of the top left (n-q) x (n-p) submatrix of AJA TSD1SV
TSD1SV
EPSM1 = EPSM*ALPHA TSD1SV
TSD1SV
CALL TSQRFC(AJA,MAXM,N-Q,N-P,1,N-Q,IERR,EPSM1,WRK3,PIVOT,ZERO1) TSD1SV
TSD1SV
IF(IERR .EQ. 0) THEN TSD1SV
ICOL = N-Q TSD1SV
ELSE TSD1SV
ICOL = ZERO1-1 TSD1SV
ENDIF TSD1SV
TSD1SV
CALL TSPRMV(D1,WRK1,PIVOT,N-Q,0) TSD1SV
TSD1SV
c solve for the first n-p components of the tensor step TSD1SV
TSD1SV
CALL TSFSLV(AJA,D1,MAXM,N-P,ICOL,WRK2) TSD1SV
TSD1SV
CALL TSQMUV(AJA,WRK2,D1,MAXM,N-P,1,ZERO1,.TRUE.) TSD1SV
TSD1SV
c copy the (n-q) x (n-p) submatrix of AJA from bottom back to TSD1SV
c top of AJA TSD1SV
TSD1SV
DO 70 J = 1,N-P TSD1SV
CALL DCOPY(N-Q,AJA(M+3,J),1,AJA(1,J),1) TSD1SV
70 CONTINUE TSD1SV
TSD1SV
RETURN TSD1SV
END TSD1SV
SUBROUTINE TSDFCN(P,X,G,AJA,ANLS,S,FN,WRK1,WRK2, TSDFCN
+ WRK3,WRK4,WRK5,MAXM,MAXN,M,N,Q) TSDFCN
TSDFCN
INTEGER P,MAXM,MAXN,M,N,Q TSDFCN
DOUBLE PRECISION X(P),G(P),AJA(MAXM,N),ANLS(MAXM,P),S(MAXN,P) TSDFCN
DOUBLE PRECISION FN(M),WRK1(M),WRK2(P),WRK3(P),WRK4(M),WRK5(M) TSDFCN
TSDFCN
C********************************************************************* TSDFCN
C THIS ROUTINE COMPUTES THE ANALYTIC GRADIENT OF THE FUNCTION GIVEN TSDFCN
C BY SUBROUTINE TSQFCN. TSDFCN
C********************************************************************* TSDFCN
C TSDFCN
C INPUT PARAMETERS : TSDFCN
C ----------------- TSDFCN
C TSDFCN
C P : COLUMN DIMENSION OF ANLS AND S TSDFCN
C X : POINT AT WHICH GRADIENT IS EVALUATED TSDFCN
C AJA: JACOBIAN MATRIX AT CURRENT ITERATE TSDFCN
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSDFCN
C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSDFCN
C FN : FUNCTION VALUE AT CURRENT ITERATE TSDFCN
C WRK1,WRK2,WRK3,WRK4,WRK5 : WORKSPACE TSDFCN
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSDFCN
C MAXN : LEADING DIMENSION OF S TSDFCN
C M,N : DIMENSIONS OF PROBLEM TSDFCN
C Q : NUMERICAL RANK OF JACOBIAN : TSDFCN
C Q > P : JACOBIAN IS SINGULAR TSDFCN
C Q = P : OTHERWISE TSDFCN
C TSDFCN
C TSDFCN
C OUTPUT PARAMETERS : TSDFCN
C ----------------- TSDFCN
C TSDFCN
C G : GRADIENT AT X TSDFCN
C TSDFCN
C SUBPROGRAMS CALLED: TSDFCN
C TSDFCN
C LEVEL 1 BLAS ... DAXPY,DDOT,DLOAD TSDFCN
C LEVEL 2 BLAS ... DGEMV TSDFCN
C TENSOLVE ... TSSTMX TSDFCN
C TSDFCN
C********************************************************************* TSDFCN
TSDFCN
INTEGER I,J,K,L TSDFCN
DOUBLE PRECISION ZERO,HALF,ONE TSDFCN
DOUBLE PRECISION DDOT TSDFCN
DATA ZERO,HALF,ONE/0.0D0,0.50D0,1.0D0/ TSDFCN
TSDFCN
c compute the lower right (m-n+q) x p submatrix of AJA times X TSDFCN
TSDFCN
CALL DGEMV('N',M-N+Q,P,ONE,AJA(N-Q+1,N-P+1),MAXM, TSDFCN
+ X,1,ZERO,WRK1,1) TSDFCN
TSDFCN
c compute S-trans times X TSDFCN
TSDFCN
CALL TSSTMX(S,X,MAXN,N,P,WRK2,WRK3) TSDFCN
TSDFCN
c compute 0.5 * (S-trans times X)**2 TSDFCN
TSDFCN
DO 10 I = 1, P TSDFCN
WRK2(I) = HALF * WRK3(I)**2 TSDFCN
10 CONTINUE TSDFCN
TSDFCN
c compute 0.5 * (lower (m-n+q) x p submatrix of ANLS) * TSDFCN
c (S-trans times X)**2 TSDFCN
TSDFCN
CALL DGEMV('N',M-N+Q,P,ONE,ANLS(N-Q+1,1),MAXM, TSDFCN
+ WRK2,1,ZERO,WRK4,1) TSDFCN
TSDFCN
DO 20 I = 1,M-N+Q TSDFCN
WRK4(I) = WRK4(I)+FN(N-Q+I)+WRK1(I) TSDFCN
20 CONTINUE TSDFCN
TSDFCN
c compute AJA-trans * WRK4 TSDFCN
TSDFCN
CALL DGEMV('T',M-N+Q,P,ONE,AJA(N-Q+1,N-P+1),MAXM, TSDFCN
+ WRK4,1,ZERO,WRK1,1) TSDFCN
TSDFCN
c compute ANLS-trans * WRK4 TSDFCN
TSDFCN
CALL DGEMV('T',M-N+Q,P,ONE,ANLS(N-Q+1,1),MAXM, TSDFCN
+ WRK4,1,ZERO,WRK5,1) TSDFCN
TSDFCN
c compute S * diag(S-trans * WRK3) * WRK5 TSDFCN
TSDFCN
CALL DLOAD(P,ZERO,WRK2,1) TSDFCN
TSDFCN
L = P+1 TSDFCN
DO 50 J = 1,P TSDFCN
L = L-1 TSDFCN
WRK2(L) = S(N+2,L) TSDFCN
DO 30 I = L+1,P TSDFCN
WRK2(I) = S(N-P+J,I) TSDFCN
30 CONTINUE TSDFCN
DO 40 K = 1,P TSDFCN
WRK2(K) = WRK2(K)*WRK3(K) TSDFCN
40 CONTINUE TSDFCN
G(J) = DDOT(P,WRK2,1,WRK5,1) TSDFCN
50 CONTINUE TSDFCN
TSDFCN
CALL DAXPY(P,ONE,WRK1,1,G,1) TSDFCN
TSDFCN
RETURN TSDFCN
END TSDFCN
SUBROUTINE TSDFLT(M,N,ITNLIM,JACFLG,GRADTL,STEPTL,FTOL,METHOD, TSDFLT
+ GLOBAL,STEPMX,DLT,TYPX,TYPF,IPR,MSG) TSDFLT
TSDFLT
INTEGER M,N,ITNLIM,JACFLG,METHOD,GLOBAL,MSG,IPR TSDFLT
DOUBLE PRECISION GRADTL,STEPTL,FTOL,STEPMX,DLT TSDFLT
DOUBLE PRECISION TYPX(N),TYPF(M) TSDFLT
TSDFLT
C********************************************************************* TSDFLT
C THIS ROUTINE SETS DEFAULT VALUES FOR EACH INPUT VARIABLE TO THE TSDFLT
C NONLINEAR EQUATION ALGORITHM. TSDFLT
C********************************************************************* TSDFLT
C TSDFLT
C SUBPROGRAMS CALLED: TSDFLT
C TSDFLT
C LEVEL 1 BLAS ... DLOAD TSDFLT
C UNCMIN ... DPMEPS TSDFLT
C TSDFLT
C********************************************************************** TSDFLT
TSDFLT
DOUBLE PRECISION EPS,DPMEPS,ONE,TWO,THREE,THOUS TSDFLT
DATA ONE,TWO,THREE,THOUS/1.0D0,2.0D0,3.0D0,1000.0D0/ TSDFLT
TSDFLT
JACFLG = 0 TSDFLT
EPS = DPMEPS() TSDFLT
GRADTL = EPS**(ONE/THREE) TSDFLT
STEPTL = EPS**(TWO/THREE) TSDFLT
FTOL = EPS**(TWO/THREE) TSDFLT
ITNLIM = 150 TSDFLT
METHOD = 1 TSDFLT
GLOBAL = 0 TSDFLT
STEPMX = THOUS TSDFLT
DLT = -ONE TSDFLT
MSG = 0 TSDFLT
IPR = 6 TSDFLT
CALL DLOAD(N,ONE,TYPX,1) TSDFLT
CALL DLOAD(M,ONE,TYPF,1) TSDFLT
TSDFLT
RETURN TSDFLT
END TSDFLT
SUBROUTINE TSDUMJ(AJA,X,NR,N) TSDUMJ
TSDUMJ
INTEGER NR, N TSDUMJ
DOUBLE PRECISION AJA(NR,N),X(N) TSDUMJ
TSDUMJ
C********************************************************************* TSDUMJ
C THIS IS A DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC TSDUMJ
C WHEN SPECIFIC ANALYTIC JACOBIAN IS NOT SUPPLIED. TSDUMJ
C********************************************************************* TSDUMJ
C TSDUMJ
C INPUT PARAMETERS: TSDUMJ
C ----------------- TSDUMJ
C TSDUMJ
C AJA : JACOBIAN MATRIX TSDUMJ
C X : POINT AT WHICH JACOBIAN IS EVALUATED TSDUMJ
C NR : LEADING DIMENSION OF AJA TSDUMJ
C N : DIMENSION OF PROBLEM TSDUMJ
C TSDUMJ
C***********************************************************************TSDUMJ
TSDUMJ
RETURN TSDUMJ
END TSDUMJ
FUNCTION TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSFAFA
+ NR,M,N,P,NWTAKE,IERR,VN) TSFAFA
TSFAFA
INTEGER NR,M,N,P,IERR TSFAFA
DOUBLE PRECISION ALPHA,DLT,TSFAFA TSFAFA
DOUBLE PRECISION ADT(N),AG(N),CONST1(P),CONST2(P) TSFAFA
DOUBLE PRECISION FQ(M),VN(M),ANLS(NR,P) TSFAFA
LOGICAL NWTAKE TSFAFA
TSFAFA
C******************************************************************** TSFAFA
C THIS FUNCTION COMPUTES || F + J*D + 0.5*A*D**2 ||**2 IN THE TSFAFA
C L2 NORM SENS, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2). TSFAFA
C******************************************************************** TSFAFA
C TSFAFA
C TSFAFA
C INPUT PARAMETERS TSFAFA
C ---------------- TSFAFA
C TSFAFA
C ANLS : TENSOR TERM MATRIX TSFAFA
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSFAFA
C ORTHOGONAL MATRICES TSFAFA
C ADT : JACOBIAN MATRIX TIMES DT TSFAFA
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSFAFA
C CONST1 : SHAT-TRANS TIMES DT TSFAFA
C CONST2 : SHAT-TRANS TIMES GBAR TSFAFA
C ALPHA : POINT AT WHICH TO EVALUATE THE FUNCTION TSFAFA TSFAFA
C DLT : CURRENT TRUST RADIUS TSFAFA
C NR : LEADING DIMENSION OF THE JACOBIAN TSFAFA
C M,N : DIMENSIONS OF THE PROBLEM TSFAFA
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSFAFA
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: TSFAFA
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSFAFA
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSFAFA
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: TSFAFA
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSFAFA
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSFAFA
C TSFAFA
C TSFAFA
C OUTPUT PARAMETERS TSFAFA
C ----------------- TSFAFA
C TSFAFA
C VN : F + J*D + 0.5*A*D**2 TSFAFA
C TSFAFA : || F + J*D + 0.5*A*D**2 ||**2 TSFAFA
C WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) TSFAFA
C TSFAFA
C SUBPROGRAMS CALLED: TSFAFA
C TSFAFA
C LEVEL 1 BLAS ... DDOT TSFAFA
C TENSOLVE ... TSMAFA TSFAFA
C TSFAFA
C******************************************************************** TSFAFA
TSFAFA
INTEGER LEN TSFAFA
DOUBLE PRECISION DDOT TSFAFA
DOUBLE PRECISION HALF TSFAFA
DATA HALF/0.50D0/ TSFAFA
TSFAFA
CALL TSMAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSFAFA
+ NR,M,N,P,NWTAKE,IERR,VN) TSFAFA
TSFAFA
LEN = M TSFAFA
IF(IERR.GT.0) LEN = M + N TSFAFA
TSFAFA
TSFAFA = HALF*DDOT(LEN,VN,1,VN,1) TSFAFA
TSFAFA
RETURN TSFAFA
END TSFAFA
SUBROUTINE TSFDFJ(XC,FC,NR,M,N,EPSM,TSFVEC,FHAT,AJA) TSFDFJ
TSFDFJ
INTEGER NR,M,N TSFDFJ
DOUBLE PRECISION EPSM TSFDFJ
DOUBLE PRECISION AJA(NR,N),FHAT(M),XC(N),FC(M) TSFDFJ
EXTERNAL TSFVEC TSFDFJ
TSFDFJ
C***********************************************************************TSFDFJ
C THIS ROUTINE COMPUTES THE FINITE DIFFERENCE JACOBIAN AT THE CURRENT TSFDFJ
C ITERATE XC. TSFDFJ
C***********************************************************************TSFDFJ
C TSFDFJ
C INPUT PARAMETERS : TSFDFJ
C ---------------- TSFDFJ
C TSFDFJ
C XC : CURRENT ITERATE TSFDFJ
C FC : FUNCTION VALUE AT XC TSFDFJ
C NR : LEADING DIMENSION OF AJA TSFDFJ
C M,N : DIMENSIONS OF PROBLEM TSFDFJ
C EPSM : MACHINE PRECISION TSFDFJ
C TSFVEC : SUBROUTINE TO EVALUATE THE USER'S FUNCTION TSFDFJ
C FHAT : WORKSPACE TSFDFJ
C TSFDFJ
C OUTPUT PARAMETERS : TSFDFJ
C -------------------- TSFDFJ
C TSFDFJ
C AJA : FINITE DIFFERENCE JACOBIAN AT XC TSFDFJ
C TSFDFJ
C SUBPROGRAMS CALLED: TSFDFJ
C TSFDFJ
C USER ... TSFVEC TSFDFJ
C TSFDFJ
C***********************************************************************TSFDFJ
TSFDFJ
INTEGER I,J TSFDFJ
DOUBLE PRECISION NDIGIT,RNOISE,STEPSZ,XTMPJ TSFDFJ
DOUBLE PRECISION SQRTR,RSTPSZ,ONE,TEN TSFDFJ
INTRINSIC ABS,MAX,SQRT TSFDFJ
DATA ONE,TEN/1.0D0,10.0D0/ TSFDFJ
TSFDFJ
NDIGIT = -LOG10(EPSM) TSFDFJ
RNOISE = MAX(TEN**(-NDIGIT),EPSM) TSFDFJ
SQRTR = SQRT(RNOISE) TSFDFJ
TSFDFJ
DO 20 J = 1,N TSFDFJ
STEPSZ = SQRTR*MAX(ABS(XC(J)),ONE) TSFDFJ
XTMPJ = XC(J) TSFDFJ
XC(J) = XTMPJ+STEPSZ TSFDFJ
CALL TSFVEC(XC,FHAT,M,N) TSFDFJ
XC(J) = XTMPJ TSFDFJ
RSTPSZ = ONE/STEPSZ TSFDFJ
DO 10 I = 1,M TSFDFJ
AJA(I,J) = (FHAT(I)-FC(I))*RSTPSZ TSFDFJ
10 CONTINUE TSFDFJ
20 CONTINUE TSFDFJ
TSFDFJ
RETURN TSFDFJ
END TSFDFJ
SUBROUTINE TSFRMT(SHAT,S,AJA,FV,FN,MAXM,MAXN,MAXP,M,N,P,IDP, TSFRMT
+ AM,X,B,SCALE,A) TSFRMT
TSFRMT
INTEGER MAXM,MAXN,MAXP,M,N,P TSFRMT
INTEGER IDP(P) TSFRMT
DOUBLE PRECISION A(MAXM,P),SHAT(MAXN,P),S(MAXN,P),AJA(MAXM,N) TSFRMT
DOUBLE PRECISION FV(MAXM,P),FN(M),AM(MAXP,P),X(P),B(P),SCALE(P) TSFRMT
TSFRMT
C********************************************************************* TSFRMT
C THIS ROUTINE FORM THE TENSOR TERM MATRIX OF THE TENSOR MODEL. TSFRMT
C********************************************************************* TSFRMT
C TSFRMT
C INPUT PARAMETERS : TSFRMT
C ---------------- TSFRMT
C TSFRMT
C SHAT: MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSFRMT
C S : MATRIX OF PREVIOUS DIRECTIONS TSFRMT
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSFRMT
C FV : MATRIX OF PAST FUNCTION VALUES TSFRMT
C FN : FUNCTION VALUE AT CURRENT ITERATE TSFRMT
C MAXM: LEADING DIMENSION OF AJA, ANLS, AND FV TSFRMT
C MAXN: LEADING DIMENSION OF S AND SHAT TSFRMT
C MAXP: LEADING DIMENSION OF AM TSFRMT
C M : ROW DIMENSION OF MATRICES A,FV,AND AJA TSFRMT
C N : COLUMN DIMENSION OF JACOBIAN MATRIX TSFRMT
C P : COLUMN DIMENSION OF MATRIX SHAT TSFRMT
C IDP : VECTOR WHICH KEEPS TRACK OF LINEARLY INDEPENDENT TSFRMT
C DIRECTION POSITIONS WITHIN THE MATRIX S TSFRMT
C AM,X,B,SCALE,: WORKSPACE TSFRMT
C TSFRMT
C OUTPUT PARAMETERS : TSFRMT
C ------------------ TSFRMT
C TSFRMT
C A : TENSOR TERM MATRIX TSFRMT
C TSFRMT
C SUBPROGRAMS CALLED: TSFRMT
C TSFRMT
C LEVEL 1 BLAS ... DDOT,DNRM2,DSCAL TSFRMT
C UNCMIN ... CHOLDC,LLTSLV TSFRMT
C TSFRMT
C********************************************************************* TSFRMT
TSFRMT
INTEGER I,J,JJ TSFRMT
DOUBLE PRECISION SUM,SC,TOL,DIAGMX,ADDMAX TSFRMT
DOUBLE PRECISION ZERO,ONE,TWO TSFRMT
DOUBLE PRECISION DDOT,DNRM2 TSFRMT
DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ TSFRMT
TSFRMT
c scale the matrix SHAT and save scaling in SCALE TSFRMT
TSFRMT
DO 10 J = 1,P TSFRMT
SC = ONE/DNRM2(N,SHAT(1,J),1) TSFRMT
CALL DSCAL(N,SC,SHAT(1,J),1) TSFRMT
SCALE(J) = SC**2 TSFRMT
10 CONTINUE TSFRMT
TSFRMT
c form the matrix AM = (Si Sj)**2 TSFRMT
TSFRMT
DO 30 J = 1,P TSFRMT
JJ = IDP(J) TSFRMT
DO 20 I = 1,P TSFRMT
AM(I,J) = DDOT(N,S(1,IDP(I)),1,S(1,JJ),1)**2 TSFRMT
20 CONTINUE TSFRMT
30 CONTINUE TSFRMT
TSFRMT
c scale the matrix AM TSFRMT
TSFRMT
DO 50 I = 1,P TSFRMT
DO 40 J = 1,P TSFRMT
AM(I,J) = SCALE(I)*SCALE(J)*AM(I,J) TSFRMT
40 CONTINUE TSFRMT
50 CONTINUE TSFRMT
TSFRMT
c perform a Cholesky decomposition of AM TSFRMT
TSFRMT
TOL = ZERO TSFRMT
DIAGMX = ZERO TSFRMT
CALL CHOLDC(MAXP,P,AM,DIAGMX,TOL,ADDMAX) TSFRMT
TSFRMT
c form the tensor term matrix A TSFRMT
TSFRMT
DO 70 I = 1,M TSFRMT
DO 60 J = 1,P TSFRMT
JJ = IDP(J) TSFRMT
SUM = DDOT(N,AJA(I,1),MAXM,S(1,JJ),1) TSFRMT
B(J) = TWO*(FV(I,JJ) - FN(I) - SUM) TSFRMT
B(J) = SCALE(J)*B(J) TSFRMT
60 CONTINUE TSFRMT
TSFRMT
c solve AM*X = B TSFRMT
TSFRMT
CALL LLTSLV(MAXP,P,AM,X,B) TSFRMT
TSFRMT
c copy X into row i of A TSFRMT
TSFRMT
CALL DCOPY(P,X,1,A(I,1),MAXM) TSFRMT
TSFRMT
70 CONTINUE TSFRMT
TSFRMT
RETURN TSFRMT
END TSFRMT
SUBROUTINE TSFSCL(X,DX,DF,M,N,TSFVEC,F) TSFSCL
TSFSCL
INTEGER M,N TSFSCL
DOUBLE PRECISION X(N),DX(N),F(M),DF(M) TSFSCL
EXTERNAL TSFVEC TSFSCL
TSFSCL
C******************************************************************** TSFSCL
C THIS ROUTINE EVALUATES THE FUNCTION AT THE CURRENT ITERATE X THEN TSFSCL
C SCALES ITS VALUE. TSFSCL
C******************************************************************** TSFSCL
C TSFSCL
C INPUT PARAMETERS : TSFSCL
C ----------------- TSFSCL
C TSFSCL
C X : CURRENT ITERATE TSFSCL
C DX : DIAGONAL SCALING MATRIX FOR X TSFSCL
C DF : DIAGONAL SCALING MATRIX FOR F TSFSCL
C M,N : DIMENSIONS OF PROBLEM TSFSCL
C TSFVEC : SUBROUTINE TO EVALUATE FUNCTION TSFSCL
C TSFSCL
C TSFSCL
C OUTPUT PARAMETERS : TSFSCL
C ----------------- TSFSCL
C TSFSCL
C F : SCALED FUNCTION VALUE AT CURRENT ITERATE X TSFSCL
C TSFSCL
C SUBPROGRAMS CALLED: TSFSCL
C TSFSCL
C TENSOLVE ... TSUNSX,TSSCLF,TSSCLX TSFSCL
C USER ... TSFVEC TSFSCL
C TSFSCL
C******************************************************************** TSFSCL
TSFSCL
CALL TSUNSX(X,DX,N) TSFSCL
CALL TSFVEC(X,F,M,N) TSFSCL
CALL TSSCLF(F,DF,M) TSFSCL
CALL TSSCLX(X,DX,N) TSFSCL
TSFSCL
RETURN TSFSCL
END TSFSCL
SUBROUTINE TSFSLV(L,B,NR,M,N,Y) TSFSLV
TSFSLV
INTEGER NR,M,N TSFSLV
DOUBLE PRECISION B(N),L(NR,N),Y(N) TSFSLV
TSFSLV
C******************************************************************** TSFSLV
C THIS ROUTINE DOES A FORWARD SOLVE. TSFSLV
C******************************************************************** TSFSLV
C TSFSLV
C INPUT PARAMETERS : TSFSLV
C ----------------- TSFSLV
C TSFSLV
C L : THE TRANSPOSE OF THE UPPER TRIANGULAR MATRIX OBTAINED TSFSLV
C FROM A QR FACTORIZATION OF AN M BY N MATRIX A. DIAG(L) TSFSLV
C IS STORED IN ROW M+2. THIS IS THE STORAGE SCHEME USED TSFSLV
C IN STEWART, G. W., III(1973) "INTRODUCTION TO MATRIX TSFSLV
C COMPUTATION", ACADEMIC PRESS,NEW YORK TSFSLV
C B : RIGHT HAND SIDE TSFSLV
C NR : LEADING DIMENSION OF MATRIX A TSFSLV
C M : ROW DIMENSION OF MATRIX A TSFSLV
C N : COLUMN DIMENSION OF MATRIX A TSFSLV
C TSFSLV
C OUTPUT PARAMETERS : TSFSLV
C ------------------ TSFSLV
C TSFSLV
C Y : VECTOR SOLUTION ON EXIT TSFSLV
C TSFSLV
C SUBPROGRAMS CALLED: TSFSLV
C TSFSLV
C LEVEL 1 BLAS ... DDOT TSFSLV
C TSFSLV
C********************************************************************* TSFSLV
TSFSLV
INTEGER J TSFSLV
DOUBLE PRECISION S TSFSLV
DOUBLE PRECISION DDOT TSFSLV
TSFSLV
c solve L Y = B TSFSLV
TSFSLV
Y(1) = B(1) / L(M+2,1) TSFSLV
IF(N .GT. 1) THEN TSFSLV
S = L(1,2) * Y(1) TSFSLV
Y(2) = (B(2) - S) / L(M+2,2) TSFSLV
DO 10 J = 3,N TSFSLV
S = DDOT(J-1,L(1,J),1,Y,1) TSFSLV
Y(J) = (B(J) - S) / L(M+2,J) TSFSLV
10 CONTINUE TSFSLV
ENDIF TSFSLV
TSFSLV
RETURN TSFSLV
END TSFSLV
SUBROUTINE TSJMUV(ITN,METHOD,V,CURPOS,PIVOT,PBAR,AJA,SHAT, TSJMUV
+ FLAG,IERR,MAXM,MAXN,M,N,P,WRK1,WRK2,VN,AV) TSJMUV
TSJMUV
INTEGER MAXM,MAXN,M,N,P,IERR,ITN,METHOD,FLAG TSJMUV
INTEGER CURPOS(N),PIVOT(N),PBAR(N) TSJMUV
DOUBLE PRECISION V(N),WRK1(N),WRK2(N),VN(M),AJA(MAXM,N) TSJMUV
DOUBLE PRECISION AV(N),SHAT(MAXN,P) TSJMUV
TSJMUV
C**************************************************************** TSJMUV
C THIS ROUTINE CALCULATES THE PRODUCT JACOBIAN TIMES A VECTOR. TSJMUV
C**************************************************************** TSJMUV
C TSJMUV
C INPUT PARAMETERS TSJMUV
C ---------------- TSJMUV
C TSJMUV
C ITN : CURRENT ITERATION NUMBER TSJMUV
C METHOD : METHOD TO BE USED TSJMUV
C V : VECTOR TO BE MULTIPLIED BY AJA TSJMUV
C CURPOS : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA TSJMUV
C FROM COLUMN 1 TO N-P) TSJMUV
C PIVOT : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA TSJMUV
C FROM COLUMN N-P+1 TO N) TSJMUV
C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA TSJMUV
C IF IT IS SINGULAR TSJMUV
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSJMUV
C SHAT : MATRIX OF LINEARLY INDEPENDENT DIRECTIONS AFTER TSJMUV
C A QL FACTORIZATION TSJMUV
C FLAG : RETURN CODE WITH THE FOLLOWING MEANINGS: TSJMUV
C FLAG = 0 : NO SINGULARITY DETECTED DURING FACTORIZATIONTSJMUV
C OF THE JACOBIAN FROM COLUMN 1 TO N TSJMUV
C FLAG = 1 : SINGULARITY DETECTED DURING FACTORIZATION TSJMUV
C OF THE JACOBIAN FROM COLUMN 1 TO N-P TSJMUV
C FLAG = 2 : SINGULARITY DETECTED DURING FACTORIZATION TSJMUV
C OF THE JACOBIAN FROM COLUMN N-P+1 TO N TSJMUV
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: TSJMUV
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSJMUV
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSJMUV
C MAXM : LEADING DIMENSION OF AJA TSJMUV
C MAXN : LEADING DIMENSION OF SHAT TSJMUV
C M,N : DIMENSIONS OF THE PROBLEM TSJMUV
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSJMUV
C WRK1,WRK2,VN : WORKSPACE VECTORS TSJMUV
C TSJMUV
C OUTPUT PARAMETERS TSJMUV
C ----------------- TSJMUV
C TSJMUV
C AV : JACOBIAN TIMES V TSJMUV
C TSJMUV
C SUBPROGRAMS CALLED: TSJMUV
C TSJMUV
C LEVEL 1 BLAS ... DCOPY TSJMUV
C TENSOLVE ... TSPRMV,TSQMLV,TSUTMD TSJMUV
C TSJMUV
C **********************************************************************TSJMUV
TSJMUV
INTEGER LEN TSJMUV
IF(ITN .EQ. 1 .OR. METHOD .EQ. 0) THEN TSJMUV
CALL TSPRMV(WRK1,V,PIVOT,N,1) TSJMUV
IF(IERR .EQ. 1) THEN TSJMUV
CALL TSPRMV(WRK2,WRK1,PBAR,N,1) TSJMUV
CALL DCOPY(N,WRK2,1,WRK1,1) TSJMUV
ENDIF TSJMUV
ELSEIF(N .EQ. 1) THEN TSJMUV
VN(1) = V(1) TSJMUV
ELSE TSJMUV
CALL TSQMLV(MAXN,N,P,SHAT,V,VN,.FALSE.) TSJMUV
CALL TSPRMV(WRK2,VN,CURPOS,N,1) TSJMUV
IF(FLAG .EQ. 0) THEN TSJMUV
CALL TSPRMV(WRK1,WRK2,PIVOT,N,1) TSJMUV
ELSEIF(FLAG .EQ. 1) THEN TSJMUV
CALL TSPRMV(WRK1,WRK2,PBAR,N,1) TSJMUV
ELSEIF(FLAG .EQ. 2 ) THEN TSJMUV
CALL TSPRMV(WRK1,WRK2,PIVOT,N,1) TSJMUV
CALL TSPRMV(WRK2,WRK1,PBAR,N,1) TSJMUV
CALL DCOPY(N,WRK2,1,WRK1,1) TSJMUV
ENDIF TSJMUV
ENDIF TSJMUV
TSJMUV
LEN = M TSJMUV
IF(IERR .GT. 0) LEN = M + N TSJMUV
TSJMUV
CALL TSUTMD(AJA,WRK1,MAXM,LEN,N,AV) TSJMUV
TSJMUV
RETURN TSJMUV
END TSJMUV
SUBROUTINE TSJQTP(Q,MAXM,MAXN,N,M,P,WRK1,WRK2,AJA) TSJQTP
TSJQTP
INTEGER MAXM,MAXN,N,M,P TSJQTP
DOUBLE PRECISION AJA(MAXM,N),Q(MAXN,P),WRK1(N),WRK2(N) TSJQTP
TSJQTP
C********************************************************************** TSJQTP
C THIS ROUTINE GETS J*(Q-TRANS) BY COMPUTING EACH ROW OF THE TSJQTP
C RESULTING MATRIX AS FOLLOWS : (J*Q-TRANS)I-TH ROW<--Q*(J)I-TH ROW. TSJQTP
C********************************************************************** TSJQTP
C TSJQTP
C INPUT PARAMETERS : TSJQTP
C ----------------- TSJQTP
C TSJQTP
C Q : RESULTING MATRIX FROM A QL FACTORIZATION TSJQTP
C MAXM : LEADING DIMENSION OF AJA TSJQTP
C MAXN : LEADING DIMENSION OF Q TSJQTP
C M,N : DIMENSIONS OF PROBLEM TSJQTP
C P : COLUMN DIMENSION OF MATRIX Q TSJQTP
C WRK1,WRK2: WORKING VECTOR TSJQTP
C TSJQTP
C INPUT-OUTPUT PARAMETERS : TSJQTP
C ------------------------ TSJQTP
C TSJQTP
C AJA : JACOBIAN MATRIX ON ENTRY AND JACOBIAN MULTIPLIED BY THE TSJQTP
C ORTHOGONAL MATRIX Q ON EXIT TSJQTP
C TSJQTP
C SUBPROGRAMS CALLED: TSJQTP
C TSJQTP
C LEVEL 1 BLAS ... DCOPY TSJQTP
C TENSOLVE ... TSQMLV TSJQTP
C TSJQTP
C********************************************************************** TSJQTP
TSJQTP
INTEGER I TSJQTP
TSJQTP
DO 30 I = 1,M TSJQTP
TSJQTP
c copy the i-th row of AJA into WRK1 TSJQTP
TSJQTP
CALL DCOPY(N,AJA(I,1),MAXM,WRK1,1) TSJQTP
TSJQTP
CALL TSQMLV(MAXN,N,P,Q,WRK1,WRK2,.FALSE.) TSJQTP
TSJQTP
c form the i-th row of AJA*(Q-trans) TSJQTP
TSJQTP
CALL DCOPY(N,WRK2,1,AJA(I,1),MAXM) TSJQTP
TSJQTP
30 CONTINUE TSJQTP
TSJQTP
RETURN TSJQTP
END TSJQTP
SUBROUTINE TSLMIN(XC,XP,P1,Q,ANLS,FQ,ADT,AG,CONST1,CONST2, TSLMIN
+ DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP,VNS,XPLUS) TSLMIN
TSLMIN
INTEGER NR,M,N,P,IERR TSLMIN
DOUBLE PRECISION XC,XP,XPLUS,P1,Q,DLT,TOL TSLMIN
DOUBLE PRECISION ADT(N),AG(N),VN(M),VNP(M),VNS(M) TSLMIN
DOUBLE PRECISION ANLS(NR,P),FQ(M),CONST1(P),CONST2(P) TSLMIN
LOGICAL NWTAKE TSLMIN
TSLMIN
C***********************************************************************TSLMIN
C THIS ROUTINE FINDS A LOCAL MINIMIZER OF A ONE-VARIABLE FUNCTION IN AN TSLMIN
C INTERVAL [XC XP]. TSLMIN
C***********************************************************************TSLMIN
C TSLMIN
C INPUT PARAMETERS : TSLMIN
C ----------------- TSLMIN
C TSLMIN
C XC,XP : LOWER AND UPPER BOUND OF INTERVAL IN WHICH THE SEARCH TSLMIN
C IS PERFORMED TSLMIN
C P1,Q : FIRST DERIVATIVES OF THE ONE-VARIABLE FUNCTION TSLMIN
C ANLS : TENSOR TERM MATRIX TSLMIN
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSLMIN
C ORTHOGONAL MATRICES TSLMIN
C ADT : JACOBIAN TIMES THE STEP DT (SEE SUBROUTINE TS2DTR) TSLMIN
C AG : JACOBIAN TIMES THE GRADIENT G (SEE SUBROUTINE TS2DTR) TSLMIN
C CONST1 : SHAT-TRANS * DT (SEE SUBROUTINE TS2DTR) TSLMIN
C CONST2 : SHAT-TRANS * GBAR (SEE SUBROUTINE TS2DTR) TSLMIN
C DLT : TRUST RADIUS TSLMIN
C NR : LEADING DIMENSION OF ANLS MATRIX TSLMIN
C M,N : DIMENSIONS OF PROBLEM TSLMIN
C P : COLUMN DIMENSION OF MATRIX ANLS TSLMIN
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: TSLMIN
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSLMIN
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSLMIN
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: TSLMIN
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSLMIN
C IERR = 1 : OTHERWISE TSLMIN
C TOL : SMALL TOLERANCE TSLMIN
C VN,VNP,VNS : WORKING VECTORS TSLMIN
C TSLMIN
C TSLMIN
C OUTPUT PARAMETERS : TSLMIN
C ----------------- TSLMIN
C TSLMIN
C XPLUS : LOCAL MINIMIZER OF THE ONE-VARIABLE FUNCTION TSLMIN
C TSLMIN
C SUBPROGRAMS CALLED : TSLMIN
C TSLMIN
C TENSOLVE ... TSMSDA,TSFAFA,TSLMSP,TSMFDA TSLMIN
C TSLMIN
C***********************************************************************TSLMIN
TSLMIN
INTEGER ITERCD,RETCD,ITNCNT TSLMIN
DOUBLE PRECISION ALEFT,ARIGHT,T,E,TSMSDA,S,SINIT,TSFAFA,TSMFDA TSLMIN
DOUBLE PRECISION ZERO,OTT,TWO,SMALL TSLMIN
LOGICAL SKIP TSLMIN
INTRINSIC ABS,MIN,MAX TSLMIN
DATA ZERO,OTT,TWO,SMALL/0.0D0,1.0D-4,2.0D0,2.0D-20/ TSLMIN
TSLMIN
ALEFT = MIN(XC,XP) TSLMIN
ARIGHT = MAX(XC,XP) TSLMIN
ITNCNT = 0 TSLMIN
T = ABS(XC-XP) TSLMIN
SKIP = .FALSE. TSLMIN
TSLMIN
c compute the second derivative value at the current point TSLMIN
TSLMIN
E = TSMSDA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC,DLT, TSLMIN
+ NR,M,N,P,NWTAKE,IERR,SKIP,VN,VNP,VNS) TSLMIN
TSLMIN
10 IF(E.GT.ZERO) THEN TSLMIN
S = -P1/E TSLMIN
IF(ABS(S).GT.TWO*T) THEN TSLMIN
IF (S.LT.ZERO) THEN TSLMIN
S = -TWO*T TSLMIN
ELSE TSLMIN
S = TWO*T TSLMIN
ENDIF TSLMIN
ENDIF TSLMIN
ELSE TSLMIN
IF (P1.GT.ZERO) THEN TSLMIN
S = -T TSLMIN
ELSE TSLMIN
S = T TSLMIN
ENDIF TSLMIN
ENDIF TSLMIN
TSLMIN
IF(XC+S.GT.ARIGHT) S = ARIGHT-XC TSLMIN
IF(XC+S.LT.ALEFT) S = ALEFT-XC TSLMIN
SINIT = ABS(S) TSLMIN
TSLMIN
20 CONTINUE TSLMIN
TSLMIN
c compute a next iterate XPLUS TSLMIN
TSLMIN
IF (TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC+S,DLT, TSLMIN
+ NR,M,N,P,NWTAKE,IERR,VN).GT.Q + OTT*S*P1) THEN TSLMIN
S = S/2 TSLMIN
IF(ABS(S).LT.SMALL*SINIT.OR.S.EQ.ZERO) THEN TSLMIN
RETCD = 1 TSLMIN
ELSE TSLMIN
GO TO 20 TSLMIN
ENDIF TSLMIN
ENDIF TSLMIN
TSLMIN
XPLUS = XC+S TSLMIN
ITNCNT = ITNCNT+1 TSLMIN
TSLMIN
c check stopping criteria TSLMIN
TSLMIN
CALL TSLMSP(XC,XPLUS,ITNCNT,RETCD,ITERCD,ANLS,ADT,AG, TSLMIN
+ CONST1,CONST2,DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP) TSLMIN
TSLMIN
IF(ITERCD.GT.0) RETURN TSLMIN
TSLMIN
c update XC TSLMIN
TSLMIN
XC = XPLUS TSLMIN
TSLMIN
c compute function and derivative values at the new point TSLMIN
TSLMIN
Q = TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC,DLT, TSLMIN
+ NR,M,N,P,NWTAKE,IERR,VN) TSLMIN
P1 = TSMFDA(ANLS,ADT,AG,CONST1,CONST2,XC,DLT, TSLMIN
+ NR,M,N,P,NWTAKE,IERR,VN,VNP) TSLMIN
SKIP = .TRUE. TSLMIN
E = TSMSDA(ANLS,FQ,ADT,AG,CONST1,CONST2,XC,DLT, TSLMIN
+ NR,M,N,P,NWTAKE,IERR,SKIP,VN,VNP,VNS) TSLMIN
GO TO 10 TSLMIN
TSLMIN
END TSLMIN
SUBROUTINE TSLMSP(XC,XP,ITNCNT,RETCD,ITERCD,ANLS,ADT,AG,CONST1, TSLMSP
+ CONST2,DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP) TSLMSP
TSLMSP
INTEGER NR,M,N,P,IERR,RETCD,ITERCD,ITNCNT TSLMSP
DOUBLE PRECISION XC,XP,DLT,TOL TSLMSP
DOUBLE PRECISION ADT(N),AG(N),CONST1(P) TSLMSP
DOUBLE PRECISION CONST2(P),VN(M),VNP(M),ANLS(NR,P) TSLMSP
LOGICAL NWTAKE TSLMSP
TSLMSP
C***********************************************************************TSLMSP
C THIS ROUTINE CHECKS THE STOPPING CRITERIA FOR A LOCAL MINIMIZER. TSLMSP
C***********************************************************************TSLMSP
C TSLMSP
C INPUT PARAMETERS : TSLMSP
C ----------------- TSLMSP
C TSLMSP
C XC : CURRENT ITERATE (FROM SEARCH SUBROUTINE) TSLMSP
C XP : NEXT ITERATE (FROM SEARCH SUBROUTINE) TSLMSP
C ITNCNT : ITERATION LIMIT TSLMSP
C RETCD : RETURN CODE FROM LINE SEARCH TSLMSP
C DLT : TRUST RADIUS TSLMSP
C AJA : JACOBIAN AT THE CURRENT ITERATE TSLMSP
C NR : LEADING DIMENSION OF THE JACOBIAN MATRIX TSLMSP
C M,N : DIMENSIONS OF THE PROBLEM TSLMSP
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSLMSP
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS : TSLMSP
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSLMSP
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSLMSP
C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE : TSLMSP
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSLMSP
C IERR = 1 : OTHERWISE TSLMSP
C TOL : SMALL TOLERANCE TSLMSP
C METHOD : METHOD TO USE TSLMSP
C = 0 : STANDARD METHOD USED TSLMSP
C = 1 : TENSOR METHOD USED TSLMSP
C VN,VNP : WORKING VECTORS TSLMSP
C TSLMSP
C TSLMSP
C OUTPUT PARAMETERS : TSLMSP
C ------------------ TSLMSP
C TSLMSP
C ITERCD : RETURN CODE WITH FOLLOWING MEANINGS : TSLMSP
C ITERCD = 1 : FIRST DERIVATIVE AT THE CURRENT POINT TSLMSP
C CLOSE TO 0 TSLMSP
C ITERCD = 2 : SUCCESSIVE ITERATES WITHIN TOLERANCE TSLMSP
C ITERCD = 3 : LINE SEARCH FAILED TO LOCATE A POINT TSLMSP
C LOWER THAT THE CURRENT POINT TSLMSP
C ITERCD = 4 : ITERATION LIMIT EXCEEDED TSLMSP
C TSLMSP
C***********************************************************************TSLMSP
TSLMSP
DOUBLE PRECISION TSMFDA,GRDT,ZERO TSLMSP
INTRINSIC DABS,SQRT TSLMSP
DATA ZERO/0.0D0/ TSLMSP
TSLMSP
GRDT = SQRT(TOL) TSLMSP
ITERCD = 0 TSLMSP
TSLMSP
IF(RETCD.EQ.1) THEN TSLMSP
ITERCD = 3 TSLMSP
ELSEIF(DABS(TSMFDA(ANLS,ADT,AG,CONST1,CONST2,XP,DLT, TSLMSP
+ NR,M,N,P,NWTAKE,IERR,VN,VNP)) .LT. GRDT) THEN TSLMSP
ITERCD = 1 TSLMSP
ELSEIF(XP.NE.ZERO .AND. DABS(XP-XC)/DABS(XP).LE.TOL) THEN TSLMSP
ITERCD = 2 TSLMSP
ELSEIF(ITNCNT.GE.150) THEN TSLMSP
ITERCD = 4 TSLMSP
ENDIF TSLMSP
TSLMSP
RETURN TSLMSP
END TSLMSP
SUBROUTINE TSLSCH(M,N,XC,FC,D,G,STEPTL,DX,DF,TSFVEC, TSLSCH
+ MXTAKE,STEPMX,XP,FP,FCNORM,FPNORM,RETCD) TSLSCH
TSLSCH
INTEGER M,N,RETCD TSLSCH
DOUBLE PRECISION STEPTL,FCNORM,FPNORM TSLSCH
DOUBLE PRECISION XC(N),FC(M) TSLSCH
DOUBLE PRECISION D(N),G(N),XP(N),FP(M) TSLSCH
DOUBLE PRECISION DX(N),DF(M),STEPMX TSLSCH
LOGICAL MXTAKE TSLSCH
EXTERNAL TSFVEC TSLSCH
TSLSCH
C********************************************************************** TSLSCH
C THIS ROUTINE FINDS A NEXT ITERATE USING A STANDARD LINE SEARCH METHOD.TSLSCH
C********************************************************************** TSLSCH
C TSLSCH
C INPUT PARAMETERS : TSLSCH
C ----------------- TSLSCH
C TSLSCH
C M,N : DIMENSIONS OF PROBLEM TSLSCH
C XC : CURRENT ITERATE TSLSCH
C FC : FUNCTION VALUE AT CURRENT ITERATE TSLSCH
C D : SEARCH DIRECTION TSLSCH
C G : GRADIENT AT CURRENT ITERATE TSLSCH
C STEPTL : RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES TSLSCH
C ARE CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM TSLSCH
C DX : DIAGONAL SCALING MATRIX FOR X TSLSCH
C DF : DIAGONAL SCALING MATRIX FOR F TSLSCH
C TSFVEC: SUBROUTINE TO EVALUATE THE FUNCTION TSLSCH
C STEPMX: MAXIMUM ALLOWABLE STEP SIZE TSLSCH
C TSLSCH
C TSLSCH
C OUTPUT PARAMETERS : TSLSCH
C ----------------- TSLSCH
C TSLSCH
C MXTAKE: BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED TSLSCH
C XP : NEXT ITARATE TSLSCH
C FP : FUNCTION VALUE AT NEXT ITERATE TSLSCH
C FCNORM : 0.5 * || F(XC) ||**2 TSLSCH
C FPNORM : 0.5 * || F(XP) ||**2 TSLSCH
C RETCD : RETURN CODE WITH THE FOLLOWING MEANING : TSLSCH
C RETCD = 0 : SATISFACTORY LOCATION OF A NEW ITERATE TSLSCH
C RETCD = 1 : NO SATISFACTORY POINT FOUND SUFFICIENTLY TSLSCH
C DISTINCT FROM X TSLSCH
C TSLSCH
C SUBPROGRAMS CALLED: TSLSCH
C TSLSCH
C LEVEL 1 BLAS ... DDOT,DNRM2 TSLSCH
C TENSOLVE ... TSFSCL TSLSCH
C USER ... TSFVEC TSLSCH
C TSLSCH
C********************************************************************** TSLSCH
TSLSCH
INTEGER I,COUNT TSLSCH
DOUBLE PRECISION ALPHA,SLOPE,RELENG TSLSCH
DOUBLE PRECISION TEMP1,TEMP2,ALMDA,TEMP,ALMDAT,ALMDAM TSLSCH
DOUBLE PRECISION SLN,SCL TSLSCH
DOUBLE PRECISION DDOT,DNRM2 TSLSCH
DOUBLE PRECISION ZERO,TENTH,HALF,Z99,ONE,TWO,TEN TSLSCH
INTRINSIC ABS TSLSCH
PARAMETER (ALPHA = 1.0D-4) TSLSCH
DATA ZERO,TENTH,HALF,Z99,ONE,TWO,TEN/0.0D0,0.10D0,0.50D0,0.99D0, TSLSCH
+ 1.0D0,2.0D0,10.0D0/ TSLSCH
TSLSCH
MXTAKE = .FALSE. TSLSCH
SLN = DNRM2(N,D,1) TSLSCH
IF(SLN .GT. STEPMX) THEN TSLSCH
TSLSCH
c step longer than maximum allowed TSLSCH
TSLSCH
SCL = STEPMX/SLN TSLSCH
CALL DSCAL(N,SCL,D,1) TSLSCH
SLN = STEPMX TSLSCH
ENDIF TSLSCH
TSLSCH
c compute SLOPE = G-trans * D TSLSCH
TSLSCH
SLOPE = DDOT(N,G,1,D,1) TSLSCH
TSLSCH
c initialization of RETCD TSLSCH
TSLSCH
RETCD = 0 TSLSCH
TSLSCH
c compute the smallest value allowable for the damping TSLSCH
c parameter ALMDA, i.e, ALMDAM TSLSCH
TSLSCH
RELENG = ZERO TSLSCH
DO 20 I = 1,N TSLSCH
TEMP1 = MAX(ABS(XC(I)), ONE) TSLSCH
TEMP2 = ABS(D(I))/TEMP1 TSLSCH
RELENG = MAX(RELENG,TEMP2) TSLSCH
20 CONTINUE TSLSCH
ALMDAM = STEPTL/RELENG TSLSCH
ALMDA = ONE TSLSCH
TSLSCH
COUNT=0 TSLSCH
TSLSCH
40 CONTINUE TSLSCH
TSLSCH
COUNT=COUNT+1 TSLSCH
TSLSCH
c compute the next iterate XP TSLSCH
TSLSCH
DO 50 I = 1,N TSLSCH
XP(I) = XC(I)+ALMDA*D(I) TSLSCH
50 CONTINUE TSLSCH
TSLSCH
c evaluate the objective function at XP and its residual TSLSCH
TSLSCH
CALL TSFSCL(XP,DX,DF,M,N,TSFVEC,FP) TSLSCH
TSLSCH
FPNORM = HALF*DNRM2(M,FP,1)**2 TSLSCH
TSLSCH
c test whether the full step produces enough decrease in TSLSCH
c the l2 norm of the objective function. If not update ALMDA TSLSCH
c and compute a new step TSLSCH
TSLSCH
IF (FPNORM.GT.(FCNORM + (ALPHA* ALMDA * SLOPE))) THEN TSLSCH
ALMDAT = ((-ALMDA**2)*SLOPE)/(TWO*(FPNORM-FCNORM-ALMDA*SLOPE))TSLSCH
TEMP = ALMDA/TEN TSLSCH
IF (COUNT.GT.2000) THEN TSLSCH
RETCD = 1 TSLSCH
RETURN TSLSCH
ENDIF TSLSCH
IF (ALMDAT.GT.ONE) THEN TSLSCH
ALMDA=TEMP TSLSCH
GO TO 40 TSLSCH
END IF TSLSCH
ALMDA = MAX(TEMP,ALMDAT) TSLSCH
IF(ALMDA.LT.ALMDAM) THEN TSLSCH
RETCD = 1 TSLSCH
RETURN TSLSCH
ENDIF TSLSCH
GO TO 40 TSLSCH
ELSE TSLSCH
IF(ALMDA.EQ.TENTH .AND. SLN.GT.Z99*STEPMX) MXTAKE=.TRUE. TSLSCH
ENDIF TSLSCH
TSLSCH
RETURN TSLSCH
END TSLSCH
SUBROUTINE TSMAFA(ANLS,F,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSMAFA
+ NR,M,N,P,NWTAKE,IERR,VN) TSMAFA
TSMAFA
INTEGER NR,M,N,P,IERR TSMAFA
DOUBLE PRECISION ALPHA,DLT TSMAFA
DOUBLE PRECISION ADT(N),AG(N),CONST1(P) TSMAFA
DOUBLE PRECISION CONST2(P),F(M),VN(M),ANLS(NR,P) TSMAFA
LOGICAL NWTAKE TSMAFA
TSMAFA
C***********************************************************************TSMAFA
C THIS ROUTINE COMPUTES THE VECTOR VN = F(XC) + J(XC)*D + 0.5*A*D**2, TSMAFA
C WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2). TSMAFA
C***********************************************************************TSMAFA
C TSMAFA
C TSMAFA
C INPUT PARAMETERS : TSMAFA
C ----------------- TSMAFA
C TSMAFA
C ANLS : TENSOR TERM MATRIX TSMAFA
C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) TSMAFA
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMAFA
C CONST1: SHAT-TRANS * DT (SEE SUBROUTINE TS2DTR) TSMAFA
C CONST2: SHAT-TRABS * GBAR (SEE SUBROUTINE TS2DTR) TSMAFA
C ALPHA : POINT AT WHICH DERIVATIVE IS EVALUATED TSMAFA
C DLT : CURRENT TRUST RADIUS TSMAFA
C NR : LEADING DIMENSION OF ANLS TSMAFA
C M,N : DIMENSIONS OF THE PROBLEM TSMAFA
C P : COLUMN DIMENSION OF THE MATRIX ANLS TSMAFA
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS TSMAFA
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSMAFA
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSMAFA
C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE : TSMAFA
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSMAFA
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSMAFA
C TSMAFA
C OUTPUT PARAMETERS : TSMAFA
C ------------------- TSMAFA
C TSMAFA
C VN : F + J*D + 0.5*A*D**2, WHERE TSMAFA
C D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) TSMAFA
C TSMAFA
C SUBPROGRAMS CALLED: TSMAFA
C TSMAFA
C LEVEL 1 BLAS ... DLOAD,DAXPY TSMAFA
C TSMAFA
C******************************************************************* TSMAFA
TSMAFA
INTEGER I,J,LEN TSMAFA
DOUBLE PRECISION EXPR,CONST,ZERO,HALF TSMAFA
INTRINSIC SQRT TSMAFA
DATA ZERO,HALF/0.0D0,0.50D0/ TSMAFA
TSMAFA
EXPR = SQRT(DLT**2 - ALPHA**2) TSMAFA
DO 10 I = 1,N TSMAFA
VN(I) = ALPHA*ADT(I) + EXPR*AG(I) TSMAFA
10 CONTINUE TSMAFA
TSMAFA
CALL DLOAD (M,ZERO,VN(N+1),1) TSMAFA
TSMAFA
LEN = M TSMAFA
IF(IERR .GT. 0) LEN = M + N TSMAFA
TSMAFA
DO 30 I = 1, LEN TSMAFA
VN(I) = VN(I) + F(I) TSMAFA
30 CONTINUE TSMAFA
TSMAFA
IF(NWTAKE) RETURN TSMAFA
DO 70 J = 1,P TSMAFA
CONST = HALF*(ALPHA*CONST1(J) + EXPR*CONST2(J))**2 TSMAFA
CALL DAXPY(LEN,CONST,ANLS(1,J),1,VN,1) TSMAFA
70 CONTINUE TSMAFA
TSMAFA
RETURN TSMAFA
END TSMAFA
SUBROUTINE TSMDLS(AJA,SHAT,ANLS,XC,FC,M,N,MAXM,MAXN,P,DT,G, TSMDLS
+ DX,DF,TSFVEC,METHOD,STEPTL,GLOBAL,STEPMX, TSMDLS
+ EPSM,FQ,WRK1,WRK2,WRK3,WRK4,DN,FQQ,PIVOT, TSMDLS
+ CURPOS,PBAR,MXTAKE,XP,FP,FCNORM,FPNORM, TSMDLS
+ ZERO1,RETCD,IERR) TSMDLS
TSMDLS
INTEGER M,N,MAXM,MAXN,P,METHOD,GLOBAL,ZERO1,RETCD,IERR TSMDLS
INTEGER PIVOT(N),PBAR(N),CURPOS(N) TSMDLS
DOUBLE PRECISION STEPTL,STEPMX,EPSM,FCNORM,FPNORM TSMDLS
DOUBLE PRECISION AJA(MAXM,N),SHAT(MAXN,P),ANLS(MAXM,P) TSMDLS
DOUBLE PRECISION XC(N),FC(M),DT(N),G(N),DX(N),DF(M),FQ(M) TSMDLS
DOUBLE PRECISION WRK1(M),WRK2(M),WRK3(M),WRK4(N) TSMDLS
DOUBLE PRECISION DN(N),FQQ(M),XP(N),FP(M) TSMDLS
LOGICAL MXTAKE TSMDLS
EXTERNAL TSFVEC TSMDLS
TSMDLS
C********************************************************************** TSMDLS
C THIS ROUTINE FINDS A NEXT ITERATE USING A LINE SEARCH METHOD. IT TSMDLS
C TRIES THE FULL TENSOR STEP FIRST. IF THIS IS NOT SUCCESSFUL THEN TSMDLS
C IT COMPUTES THE STANDARD DIRECTION AND COMPUTES A STEP IN THAT TSMDLS
C DIRECTION. NEXT, IF THE TENSOR DIRECTION IS DESCENT, IT COMPUTES TSMDLS
C A STEP IN THE TENSOR DIRECTION. THE ITERATE THAT PRODUCES TSMDLS
C THE LOWER RESIDUAL IS THE NEXT ITERATE FOR THE NONLINEAR ALGORITHM. TSMDLS
C********************************************************************** TSMDLS
C TSMDLS
C INPUT PARAMETERS TSMDLS
C ---------------- TSMDLS
C TSMDLS
C AJA : JACOBIAN AT CURRENT ITERATE TSMDLS
C SHAT : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSMDLS
C AFTER A QL FACORIZATION TSMDLS
C ANLS : TENSOR TERM MATRIX TSMDLS
C XC : CURRENT ITERATE TSMDLS
C FC : FUNCTION VALUE AT CURRENT ITERATE TSMDLS
C M,N : DIMENSIONS OF THE PROBLEM TSMDLS
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSMDLS
C MAXN : LEADING DIMENSION OF SHAT TSMDLS
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSMDLS
C DT : TENSOR STEP TSMDLS
C G : GRADIENT AT CURRENT ITERATE TSMDLS
C DX : DIAGONAL SCALING MATRIX FOR X TSMDLS
C DF : DIAGONAL SCALING MATRIX FOR F TSMDLS
C GBAR : STEEPEST DESCENT DIRECTION (= -G) TSMDLS
C METHOD : METHOD TO USE TSMDLS
C = 0 : STANDARD METHOD USED TSMDLS
C = 1 : TENSOR METHOD USED TSMDLS
C STEPTL : STEP TOLERANCE TSMDLS
C GLOBAL : GLOBAL STRATEGY USED TSMDLS
C = 0 : LINE SEARCH IS USED TSMDLS
C = 1 : 2-DIMENSIONAL TRUST REGION IS USED TSMDLS
C STEPMX : MAXIMUM ALLOWABLE STEP SIZE TSMDLS
C EPSM : MACHINE PRECISION TSMDLS
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY AN TSMDLS
C ORTHOGOL MATRIX TSMDLS
C WRK1,WRK2,WRK3,WRK4 : WORKING VECTORS TSMDLS
C TSMDLS
C TSMDLS
C OUTPUT PARAMETERS TSMDLS
C ----------------- TSMDLS
C TSMDLS
C DN : NEWTON STEP TSMDLS
C FQQ : FQ MULTIPLIED BY AN ORTHOGONAL MATRIX TSMDLS
C CURPOS : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE TSMDLS
C JACOBIAN FROM COLUMN 1 TO N-P) TSMDLS
C PIVOT : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE TSMDLS
C JACOBIAN FROM COLUMN N-P+1 TO N) TSMDLS
C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF THE TSMDLS
C JACOBIAN IF IT IS SINGULAR TSMDLS
C MXTAKE : BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED TSMDLS
C XP : NEXT ITERATE TSMDLS
C FP : FUNCTION VALUE AT NEXT ITERATE TSMDLS
C FCNORM : 0.5 * || F(XC) ||**2 TSMDLS
C FPNORM : 0.5 * || F(XP) ||**2 TSMDLS
C ZERO1 : FIRST ZERO COLUMN OF THE JACOBIAN IN CASE OF TSMDLS
C SINGULARITY TSMDLS
C RETCD : RETURN CODE WITH THE FOLLOWING MEANING : TSMDLS
C RETCD = 0 : SATISFACTORY LOCATION OF A NEW ITERATE TSMDLS
C RETCD = 1 : NO SATISFACTORY POINT FOUND SUFFICIENTLY TSMDLS
C DISTINCT FROM X TSMDLS
C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE TSMDLS
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSMDLS
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSMDLS
C TSMDLS
C SUBPROGRAMS CALLED: TSMDLS
C TSMDLS
C LEVEL 1 BLAS ... DCOPY,DDOT,DNRM2 TSMDLS
C TENSOLVE ... TSFSCL,TSCPSS,TSLSCH TSMDLS
C TSMDLS
C***********************************************************************TSMDLS
TSMDLS
INTEGER I,FLAG,RETCD1 TSMDLS
DOUBLE PRECISION ALPHA,SLOPE,RELENG TSMDLS
DOUBLE PRECISION TEMP1,TEMP2,ALMDA,RESNEW,F1N,DTNORM,GNORM TSMDLS
DOUBLE PRECISION SLN,SCL TSMDLS
DOUBLE PRECISION BETA,TEMP,ALMDAT,ALMDAM TSMDLS
DOUBLE PRECISION DDOT,DNRM2 TSMDLS
DOUBLE PRECISION ZERO,TENTH,HALF,Z99,ONE,TWO,TEN TSMDLS
INTRINSIC ABS TSMDLS
PARAMETER (ALPHA = 1.0D-4) TSMDLS
DATA ZERO,TENTH,HALF,Z99,ONE,TWO,TEN/0.0D0,0.10D0,0.50D0,0.99D0,TSMDLS
+ 1.0D0,2.0D0,10.0D0/ TSMDLS
TSMDLS
MXTAKE = .FALSE. TSMDLS
SLN = DNRM2(N,DT,1) TSMDLS
IF(SLN .GT. STEPMX) THEN TSMDLS
TSMDLS
c step longer than maximum allowed TSMDLS
TSMDLS
SCL = STEPMX/SLN TSMDLS
CALL DSCAL(N,SCL,DT,1) TSMDLS
SLN = STEPMX TSMDLS
ENDIF TSMDLS
TSMDLS
c compute SLOPE = G-Trans * DT TSMDLS
TSMDLS
SLOPE = DDOT(N,G,1,DT,1) TSMDLS
TSMDLS
c initialization of RETCD TSMDLS
TSMDLS
RETCD = 0 TSMDLS
TSMDLS
c compute the smallest value allowable for the damping TSMDLS
c parameter ALMDA, i.e, ALMDAM TSMDLS
TSMDLS
RELENG = ZERO TSMDLS
DO 20 I = 1,N TSMDLS
TEMP1 = MAX(ABS(XC(I)), ONE) TSMDLS
TEMP2 = ABS(DT(I))/TEMP1 TSMDLS
RELENG = MAX(RELENG, TEMP2) TSMDLS
20 CONTINUE TSMDLS
ALMDAM = STEPTL/RELENG TSMDLS
ALMDA = ONE TSMDLS
TSMDLS
c compute the next iterate XP TSMDLS
TSMDLS
DO 30 I = 1,N TSMDLS
XP(I) = XC(I)+ALMDA*DT(I) TSMDLS
30 CONTINUE TSMDLS
TSMDLS
c evaluate the objective function at XP and its residual TSMDLS
TSMDLS
CALL TSFSCL(XP,DX,DF,M,N,TSFVEC,FP) TSMDLS
TSMDLS
FPNORM = HALF*DNRM2(M,FP,1)**2 TSMDLS
TSMDLS
c test whether the full tensor step produces enough decrease in the TSMDLS
c l2 norm of of the objective function TSMDLS
TSMDLS
IF (FPNORM.LT.(FCNORM + (ALPHA* ALMDA * SLOPE))) RETURN TSMDLS
TSMDLS
c compute the standard direction TSMDLS
TSMDLS
CALL TSCPSS(SHAT,MAXM,MAXN,M,N,P,METHOD,GLOBAL,EPSM,FQ, TSMDLS
+ WRK1,WRK2,WRK3,WRK4,AJA,ANLS,DN,FQQ,PIVOT, TSMDLS
+ CURPOS,PBAR,ZERO1,IERR,RESNEW,FLAG) TSMDLS
TSMDLS
c compute a step in the standard direction TSMDLS
TSMDLS
CALL TSLSCH(M,N,XC,FC,DN,G,STEPTL,DX,DF,TSFVEC, TSMDLS
+ MXTAKE,STEPMX,WRK1,WRK2,FCNORM,F1N,RETCD1) TSMDLS
TSMDLS
c test whether the tensor direction is descent TSMDLS
TSMDLS
DTNORM = DNRM2(N,DT,1) TSMDLS
GNORM = DNRM2(N,G,1) TSMDLS
IF(M.GT.N) THEN TSMDLS
BETA = TENTH TSMDLS
ELSE TSMDLS
BETA = ALPHA TSMDLS
ENDIF TSMDLS
TEMP1 = -BETA*DTNORM*GNORM TSMDLS
TSMDLS
c compute a step in the tensor direction TSMDLS
TSMDLS
IF(SLOPE .LE. TEMP1) THEN TSMDLS
50 CONTINUE TSMDLS
ALMDAT = ((-ALMDA**2)*SLOPE)/(TWO*(FPNORM-FCNORM-ALMDA*SLOPE)) TSMDLS
TEMP = ALMDA/TEN TSMDLS
ALMDA = MAX(TEMP, ALMDAT) TSMDLS
IF(ALMDA.LT.ALMDAM) THEN TSMDLS
IF(RETCD1. EQ. 1) THEN TSMDLS
RETCD = 1 TSMDLS
GO TO 70 TSMDLS
ENDIF TSMDLS
ENDIF TSMDLS
DO 60 I = 1,N TSMDLS
XP(I) = XC(I)+ALMDA*DT(I) TSMDLS
60 CONTINUE TSMDLS
CALL TSFSCL(XP,DX,DF,M,N,TSFVEC,FP) TSMDLS
FPNORM = HALF*DNRM2(M,FP,1)**2 TSMDLS
IF (FPNORM .GT.(FCNORM + (ALPHA* ALMDA * SLOPE))) GO TO 50 TSMDLS
IF(ALMDA.EQ.TENTH .AND. SLN.GT.Z99*STEPMX) MXTAKE=.TRUE. TSMDLS
70 CONTINUE TSMDLS
TSMDLS
c select the next iterate that produces the lower function value TSMDLS
TSMDLS
IF(F1N .LT. FPNORM) THEN TSMDLS
CALL DCOPY(N,WRK1,1,XP,1) TSMDLS
CALL DCOPY(M,WRK2,1,FP,1) TSMDLS
FPNORM = F1N TSMDLS
ENDIF TSMDLS
ELSE TSMDLS
CALL DCOPY(N,WRK1,1,XP,1) TSMDLS
CALL DCOPY(M,WRK2,1,FP,1) TSMDLS
FPNORM = F1N TSMDLS
ENDIF TSMDLS
TSMDLS
RETURN TSMDLS
END TSMDLS
FUNCTION TSMFDA(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSMFDA
+ NR,M,N,P,NWTAKE,IERR,VN,VNP) TSMFDA
TSMFDA
INTEGER NR,M,N,P,IERR TSMFDA
DOUBLE PRECISION ALPHA,DLT,TSMFDA TSMFDA
DOUBLE PRECISION ADT(N),AG(N),CONST1(P),CONST2(P),VN(M),VNP(M) TSMFDA
DOUBLE PRECISION ANLS(NR,P) TSMFDA
LOGICAL NWTAKE TSMFDA
TSMFDA
C***********************************************************************TSMFDA
C THIS FUNCTION COMPUTES THE DERIVATIVE OF || F + J*D + 0.5*A*D**2 ||**2TSMFDA
C IN THE L2 NORM SENS, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2). TSMFDA
C***********************************************************************TSMFDA
C TSMFDA
C TSMFDA
C INPUT PARAMETERS TSMFDA
C ---------------- TSMFDA
C TSMFDA
C ANLS : TENSOR MATRIX TSMFDA
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSMFDA
C ORTHOGONAL MATRICES TSMFDA
C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) TSMFDA
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMFDA
C CONST1 : SHAT-TRANS TIMES DT (SEE SUBROUTINE TS2DTR) TSMFDA
C CONST2 : SHAT-TRANS TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMFDA
C ALPHA : POINT AT WHICH TO EVALUATE THE DERIVATIVE OF FUNCTION TSMFDA
C DLT : CURRENT TRUST RADIUS TSMFDA
C NR : LEADING DIMENSION OF THE JACOBIAN TSMFDA
C M,N : DIMENSIONS OF THE PROBLEM TSMFDA
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSMFDA
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: TSMFDA
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSMFDA
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSMFDA
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: TSMFDA
C IERR=0 : NO SINGULARITY OF JACOBIAN DETECTED TSMFDA
C IERR=1 : SINGULARITY OF JACOBIAN DETECTED TSMFDA
C TSMFDA
C TSMFDA
C OUTPUT PARAMETERS TSMFDA
C ----------------- TSMFDA
C TSMFDA
C TSMFDA
C VN : F + J*D + 0.5*A*D**2 TSMFDA
C VNP : DERIVATIVE IN ALPHA OF F + J*D + 0.5*A*D**2 TSMFDA
C TSMFDA : DERIVATIVE IN ALPHA OF || F + J*D + 0.5*A*D**2 ||**2 TSMFDA
C WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) TSMFDA
C TSMFDA
C SUBPROGRAMS CALLED: TSMFDA
C TSMFDA
C LEVEL 1 BLAS ... DDOT TSMFDA
C TENSOLVE ... TSMFDV TSMFDA
C TSMFDA
C***********************************************************************TSMFDA
TSMFDA
INTEGER LEN TSMFDA
DOUBLE PRECISION DDOT TSMFDA
TSMFDA
CALL TSMFDV(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSMFDA
+ NR,M,N,P,NWTAKE,IERR,VNP) TSMFDA
TSMFDA
LEN = M TSMFDA
IF(IERR.GT.0) LEN = M + N TSMFDA
TSMFDA
TSMFDA = DDOT(LEN,VNP,1,VN,1) TSMFDA
TSMFDA
RETURN TSMFDA
END TSMFDA
SUBROUTINE TSMFDV(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSMFDV
+ NR,M,N,P,NWTAKE,IERR,VNP) TSMFDV
TSMFDV
INTEGER NR,M,N,P,IERR TSMFDV
DOUBLE PRECISION ALPHA,DLT TSMFDV
DOUBLE PRECISION ADT(N),AG(N),CONST1(P) TSMFDV
DOUBLE PRECISION CONST2(P),VNP(M),ANLS(NR,P) TSMFDV
LOGICAL NWTAKE TSMFDV
TSMFDV
C***********************************************************************TSMFDV
C THIS ROUTINE COMPUTES THE DERIVATIVE IN ALPHA OF THE VECTOR TSMFDV
C VN = F(XC) + J(XC)*D + 0.5*A*D**2, WHERE D = ALPHA*DT + TSMFDV
C SQRT(DLT**2-ALPHA**2). TSMFDV
C***********************************************************************TSMFDV
C TSMFDV
C TSMFDV
C INPUT PARAMETERS : TSMFDV
C ----------------- TSMFDV
C TSMFDV
C ANLS : TENSOR TERM MATRIX TSMFDV
C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) TSMFDV
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMFDV
C CONST1: SHAT-TRANS TIMES DT (SEE SUBROUTINE TS2DTR) TSMFDV
C CONST2: SHAT-TRANS TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMFDV
C ALPHA : POINT AT WHICH DERIVATIVE IS EVALUATED TSMFDV
C DLT : CURRENT TRUST RADIUS TSMFDV
C NR : LEADING DIMENSION OF ANLS TSMFDV
C M,N : DIMENSIONS OF THE PROBLEM TSMFDV
C P : COLUMN DIMENSION OF THE MATRIX ANLS TSMFDV
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS : TSMFDV
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSMFDV
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSMFDV
C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE TSMFDV
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSMFDV
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSMFDV
C TSMFDV
C TSMFDV
C OUTPUT PARAMETERS : TSMFDV
C ------------------- TSMFDV
C TSMFDV
C VNP : THE DERIVATIVE IN ALPHA OF VN = F(XC) + J(XC)*D + TSMFDV
C 0.5*A*D**2, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2) TSMFDV
C TSMFDV
C SUBPROGRAMS CALLED: TSMFDV
C TSMFDV
C LEVEL 1 BLAS ... DLOAD,DAXPY TSMFDV
C TSMFDV
C******************************************************************* TSMFDV
TSMFDV
INTEGER I,J,LEN TSMFDV
DOUBLE PRECISION QUANT1,QUANT2,EXPR,CONST TSMFDV
DOUBLE PRECISION ZERO,HALF,TWO TSMFDV
INTRINSIC SQRT TSMFDV
DATA ZERO,HALF,TWO/0.0D0,0.50D0,2.0D0/ TSMFDV
TSMFDV
QUANT1 = SQRT(DLT**2 - ALPHA**2) TSMFDV
EXPR = - ALPHA/QUANT1 TSMFDV
TSMFDV
DO 10 I = 1,N TSMFDV
VNP(I) = ADT(I) + EXPR*AG(I) TSMFDV
10 CONTINUE TSMFDV
TSMFDV
CALL DLOAD (M,ZERO,VNP(N+1),1) TSMFDV
TSMFDV
IF(NWTAKE) RETURN TSMFDV
TSMFDV
QUANT2 = QUANT1 - ALPHA**2/QUANT1 TSMFDV
TSMFDV
LEN = M TSMFDV
IF(IERR.GT.0) LEN = M + N TSMFDV
TSMFDV
DO 30 J = 1,P TSMFDV
CONST = HALF*(TWO*ALPHA*(CONST1(J)**2 - CONST2(J)**2) TSMFDV
+ +TWO*CONST1(J)*CONST2(J)*QUANT2) TSMFDV
CALL DAXPY(LEN,CONST,ANLS(1,J),1,VNP,1) TSMFDV
30 CONTINUE TSMFDV
TSMFDV
RETURN TSMFDV
END TSMFDV
SUBROUTINE TSMGSA(S,NR,N,SQRN,ITN,SHAT,P,IDP) TSMGSA
TSMGSA
INTEGER NR,N,SQRN,ITN,P TSMGSA
INTEGER IDP(SQRN) TSMGSA
DOUBLE PRECISION S(NR,SQRN),SHAT(NR,SQRN) TSMGSA
TSMGSA
C********************************************************************* TSMGSA
C THIS ROUTINE FINDS A SET OF LINEARLY INDEPENDENT DIRECTIONS USING TSMGSA
C THE MODIFIED GRAM-SCHMIDT ALGORITHM. TSMGSA
C********************************************************************* TSMGSA
C TSMGSA
C INPUT PARAMETERS : TSMGSA
C --------------- TSMGSA
C TSMGSA
C S : MATRIX OF PAST DIRECTIONS TSMGSA
C NR : LEADING DIMENSION OF THE MATRICES S AND SHAT TSMGSA
C N : ROW DIMENSION OF MATRIX S AND SHAT TSMGSA
C SQRN: MAXIMUM COLUMN DIMENSION OF SHAT TSMGSA
C ITN : CURRENT ITERATION NUMBER TSMGSA
C TSMGSA
C OUTPUT PARAMETERS : TSMGSA
C ------------------- TSMGSA
C TSMGSA
C SHAT: MATRIX OF LINEARLY INDEPENDENT DIRECTIONS TSMGSA
C P : COLUMN DIMENSION OF THE MATRIX SHAT TSMGSA
C IDP : VECTOR THAT KEEPS TRACK OF THE INDICES CORRESPONDING TO TSMGSA
C THE LINEARLY INDEPENDENT DIRECTIONS IN THE MATRIX S TSMGSA
C TSMGSA
C SUBPROGRAMS CALLED: TSMGSA
C TSMGSA
C LEVEL 1 BLAS ... DAXPY,DCOPY,DDOT,DNRM2 TSMGSA
C TSMGSA
C********************************************************************* TSMGSA
TSMGSA
INTEGER J,K,L TSMGSA
DOUBLE PRECISION TOL,TJ,SJ,SUM,RTJS,ONE,TWO TSMGSA
DOUBLE PRECISION DNRM2,DDOT TSMGSA
INTRINSIC SQRT TSMGSA
DATA ONE,TWO/1.0D0,2.0D0/ TSMGSA
TSMGSA
IF(SQRN.LT.ITN) THEN TSMGSA
K = SQRN TSMGSA
ELSE TSMGSA
K = ITN-1 TSMGSA
ENDIF TSMGSA
TSMGSA
TOL = SQRT(TWO)/TWO TSMGSA
TSMGSA
DO 10 J = 1,K TSMGSA
CALL DCOPY(N,S(1,J),1,SHAT(1,J),1) TSMGSA
10 CONTINUE TSMGSA
TSMGSA
P = 0 TSMGSA
DO 30 J = 1,K TSMGSA
TJ = DNRM2(N,SHAT(1,J),1) TSMGSA
SJ = DNRM2(N,S(1,J),1) TSMGSA
IF(TJ/SJ.GT.TOL) THEN TSMGSA
P = P+1 TSMGSA
IDP(P) = J TSMGSA
RTJS = ONE/TJ**2 TSMGSA
DO 20 L = J+1,K TSMGSA
SUM = -RTJS*DDOT(N,SHAT(1,L),1,SHAT(1,J),1) TSMGSA
CALL DAXPY(N,SUM,SHAT(1,J),1,SHAT(1,L),1) TSMGSA
20 CONTINUE TSMGSA
ENDIF TSMGSA
30 CONTINUE TSMGSA
TSMGSA
DO 40 J = 1,P TSMGSA
CALL DCOPY(N,S(1,IDP(J)),1,SHAT(1,J),1) TSMGSA
40 CONTINUE TSMGSA
TSMGSA
RETURN TSMGSA
END TSMGSA
FUNCTION TSMSDA(ANLS,FQ,ADT,AG,CONST1,CONST2,ALPHA, TSMSDA
+ DLT,NR,M,N,P,NWTAKE,IERR,SKIP,VN,VNP,VNS) TSMSDA
TSMSDA
INTEGER NR,M,N,P,IERR TSMSDA
DOUBLE PRECISION ALPHA,DLT,TSMSDA TSMSDA
DOUBLE PRECISION ADT(N),AG(N),VN(M),VNP(M) TSMSDA
DOUBLE PRECISION VNS(M),ANLS(NR,P),FQ(M) TSMSDA
DOUBLE PRECISION CONST1(P),CONST2(P) TSMSDA
LOGICAL NWTAKE TSMSDA
TSMSDA
C***********************************************************************TSMSDA
C THIS FUNCTION COMPUTES THE SECOND DERIVATIVE OF || F + J*D + TSMSDA
C 0.5*A*D**2 ||**2 IN THE L2 NORM SENS, WHERE D = ALPHA*DT + TSMSDA
C SQRT(DLT**2-ALPHA**2). TSMSDA
C***********************************************************************TSMSDA
C TSMSDA
C TSMSDA
C INPUT PARAMETERS TSMSDA
C ---------------- TSMSDA
C TSMSDA
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSMSDA
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSMSDA
C ORTHOGONAL MATRICES TSMSDA
C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) TSMSDA
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMSDA
C CONST1 : SHAT-TRANS TIMES DT (SEE SUBROUTINE TS2DTR) TSMSDA
C CONST2 : SHAT-TRANS TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMSDA
C ALPHA : POINT AT WHICH TO EVALUATE THE SECOND DERIVATIVE OF TSMSDA
C FUNCTION TSMSDA
C DLT : CURRENT TRUST RADIUS TSMSDA
C NR : LEADING DIMENSION OF THE JACOBIAN TSMSDA
C M,N : DIMENSIONS OF THE PROBLEM TSMSDA
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSMSDA
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: TSMSDA
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSMSDA
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSMSDA
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE TSMSDA
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSMSDA
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSMSDA
C TSMSDA
C TSMSDA
C OUTPUT PARAMETERS TSMSDA
C ----------------- TSMSDA
C TSMSDA
C VN : F + J*D + 0.5*A*D**2 TSMSDA
C VNP : DERIVATIVE IN ALPHA OF F + J*D + 0.5*A*D**2 TSMSDA
C VNS : SECOND DERIVATIVE IN ALPHA OF F + J*D + 0.5*A*D**2 TSMSDA
C TSMSDA : SECOND DERIVATIVE IN ALPHA OF || F + J*D + TSMSDA
C 0.5*A*D**2 ||**2 TSMSDA
C WHERE D=ALPHA*DT + SQRT(DLT**2-ALPHA**2) TSMSDA
C TSMSDA
C SUBPROGRAMS CALLED: TSMSDA
C TSMSDA
C LEVEL 1 BLAS ... DDOT TSMSDA
C TENSOLVE ... TSMAFA,TSMFDV,TSMSDV TSMSDA
C TSMSDA
C***********************************************************************TSMSDA
TSMSDA
INTEGER LEN TSMSDA
DOUBLE PRECISION DDOT TSMSDA
LOGICAL SKIP TSMSDA
TSMSDA
IF(.NOT. SKIP) THEN TSMSDA
CALL TSMAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSMSDA
+ NR,M,N,P,NWTAKE,IERR,VN) TSMSDA
CALL TSMFDV(ANLS,ADT,AG,CONST1,CONST2,ALPHA,DLT, TSMSDA
+ NR,M,N,P,NWTAKE,IERR,VNP) TSMSDA
ENDIF TSMSDA
TSMSDA
CALL TSMSDV(ANLS,AG,CONST1,CONST2,ALPHA,DLT, TSMSDA
+ NR,M,N,P,NWTAKE,IERR,VNS) TSMSDA
TSMSDA
LEN = M TSMSDA
IF(IERR.GT.0) LEN = M + N TSMSDA
TSMSDA
TSMSDA = DDOT(LEN,VNP,1,VNP,1)+DDOT(M,VNS,1,VN,1) TSMSDA
TSMSDA
RETURN TSMSDA
END TSMSDA
SUBROUTINE TSMSDV(ANLS,AG,CONST1,CONST2,ALPHA,DLT, TSMSDV
+ NR,M,N,P,NWTAKE,IERR,VNS) TSMSDV
TSMSDV
INTEGER NR,M,N,P,IERR TSMSDV
DOUBLE PRECISION ALPHA,DLT TSMSDV
DOUBLE PRECISION AG(N),CONST1(P) TSMSDV
DOUBLE PRECISION CONST2(P),VNS(M),ANLS(NR,P) TSMSDV
LOGICAL NWTAKE TSMSDV
TSMSDV
C***********************************************************************TSMSDV
C THIS ROUTINE COMPUTES THE SECOND DERIVATIVE IN ALPHA OF THE VECTOR TSMSDV
C VN = F(XC) + J(XC)*D + 0.5*A*D**2, WHERE D = ALPHA*DT + TSMSDV
C SQRT(DLT**2-ALPHA**2). TSMSDV
C***********************************************************************TSMSDV
C TSMSDV
C TSMSDV
C INPUT PARAMETERS : TSMSDV
C ----------------- TSMSDV
C TSMSDV
C ANLS : TENSOR TERM MATRIX TSMSDV
C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) TSMSDV
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSMSDV
C CONST1: SHAT-TRANS * DT (SEE SUBROUTINE TS2DTR) TSMSDV
C CONST2: SHAT-TRABS * GBAR (SEE SUBROUTINE TS2DTR) TSMSDV
C ALPHA : POINT AT WHICH DERIVATIVE IS EVALUATED TSMSDV
C DLT : CURRENT TRUST RADIUS TSMSDV
C NR : LEADING DIMENSION OF ANLS TSMSDV
C M,N : DIMENSIONS OF THE PROBLEM TSMSDV
C P : COLUMN DIMENSION OF THE MATRIX ANLS TSMSDV
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS : TSMSDV
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSMSDV
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSMSDV
C IERR : RETURN CODE FROM THE QRP FACTORIZATION ROUTINE : TSMSDV
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSMSDV
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSMSDV
C TSMSDV
C OUTPUT PARAMETERS : TSMSDV
C ------------------- TSMSDV
C TSMSDV
C VNP : THE SECOND DERIVATIVE IN ALPHA OF VN = F(XC) + J(XC)*D TSMSDV
C + 0.5*A*D**2, WHERE D = ALPHA*DT + SQRT(DLT**2-ALPHA**2)TSMSDV
C TSMSDV
C SUBPROGRAMS CALLED: TSMSDV
C TSMSDV
C LEVEL 1 BLAS ... DLOAD,DAXPY TSMSDV
C TSMSDV
C******************************************************************* TSMSDV
TSMSDV
INTEGER I,J,LEN TSMSDV
DOUBLE PRECISION QUANT1,EXPR,CONST,QUANT2 TSMSDV
DOUBLE PRECISION ZERO,HALF,ONEPF,TWO,THREE TSMSDV
INTRINSIC SQRT TSMSDV
DATA ZERO,HALF,ONEPF,TWO,THREE/0.0D0,0.50D0,1.50D0,2.0D0,3.0D0/ TSMSDV
TSMSDV
QUANT1 = DLT**2 - ALPHA**2 TSMSDV
EXPR = -DLT**2 * SQRT(QUANT1) / QUANT1**2 TSMSDV
DO 10 I = 1,N TSMSDV
VNS(I) = EXPR*AG(I) TSMSDV
10 CONTINUE TSMSDV
TSMSDV
CALL DLOAD (M,ZERO,VNS(N+1),1) TSMSDV
TSMSDV
IF(NWTAKE) RETURN TSMSDV
TSMSDV
QUANT2 = -THREE*ALPHA/SQRT(QUANT1)-ALPHA**3/QUANT1**ONEPF TSMSDV
TSMSDV
LEN = M TSMSDV
IF(IERR .GT. 0) LEN = M + N TSMSDV
TSMSDV
DO 30 J = 1,P TSMSDV
CONST = HALF*(TWO*(CONST1(J)**2 - CONST2(J)**2) TSMSDV
+ +TWO*CONST1(J)*CONST2(J)*QUANT2) TSMSDV
CALL DAXPY(LEN,CONST,ANLS(1,J),1,VNS,1) TSMSDV
30 CONTINUE TSMSDV
TSMSDV
RETURN TSMSDV
END TSMSDV
SUBROUTINE TSMSLV(AJA,S,ANLS,FC,P,MAXM,MAXN,SQRN,M,N,EPSM, TSMSLV
+ METHOD,GLOBAL,WRK1,WRK2,WRK3,WRK4,X,TYPXU, TSMSLV
+ XPLS,GPLS,A,WRK,CURPOS,PBAR,PIVOT,FQ,FQQ, TSMSLV
+ DN,DT,RESTNS,RESNEW,ITRMCD,FLAG,ZERO1,IERR) TSMSLV
TSMSLV
INTEGER MAXM,MAXN,M,N,P,GLOBAL,ZERO1,FLAG TSMSLV
INTEGER ITRMCD,IERR,MSG,ITNLIM,IPR,METHOD,SQRN TSMSLV
INTEGER PIVOT(N),PBAR(N),CURPOS(N) TSMSLV
DOUBLE PRECISION EPSM,RESTNS,RESNEW TSMSLV
DOUBLE PRECISION AJA(MAXM,N),S(MAXN,P),ANLS(MAXM,P),FQ(M),FQQ(M) TSMSLV
DOUBLE PRECISION WRK1(M),WRK2(M),WRK3(M),WRK4(M),DN(N),DT(N) TSMSLV
DOUBLE PRECISION FC(M),X(P),TYPXU(P),XPLS(P),GPLS(P),A(SQRN,P) TSMSLV
DOUBLE PRECISION WRK(SQRN,P) TSMSLV
TSMSLV
C********************************************************************** TSMSLV
C THIS ROUTINE FINDS THE TENSOR AND STANDARD STEPS. TSMSLV
C********************************************************************** TSMSLV
C TSMSLV
C INPUT PARAMETERS : TSMSLV
C --------------- TSMSLV
C TSMSLV
C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSMSLV
C P : COLUMN DIMENSION OF MATRICES ANLS AND S TSMSLV
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSMSLV
C MAXN : LEADING DIMENSION OF S TSMSLV
C SQRN : LEADING DIMENSION OF MATRICES A AND WRK TSMSLV
C M,N : DIMENSIONS OF PROBLEM TSMSLV
C EPSM : MACHINE PRECISION TSMSLV
C AJA : JACOBIAN AT CURRENT POINT XC TSMSLV
C ANLS : TENSOR TERM MATRIX AT XC TSMSLV
C FC : FUCTION VALUE XC TSMSLV
C X : ESTIMATE TO A ROOT OF FCN (USED BY UNCMIN) TSMSLV
C TYPXU: TYPICAL SIZE FOR EACH COMPONENT OF X (USED BY UNCMIN) TSMSLV
C A : WORKSPACE FOR HESSIAN (OR ESTIMATE) (USED BY UNCMIN) TSMSLV
C WRK : WORKSPACE (USED BY UNCMIN) TSMSLV
C METHOD : METHOD TO USE TSMSLV
C METHOD = 0 : STANDARD METHOD IS USED TSMSLV
C METHOD = 1 : TENSOR METHOD IS USED TSMSLV
C GLOBAL : GLOBAL STRATEGY USED TSMSLV
C WRK1,WRK2,WRK3,WRK4,FQ,FQQ,WRK3 : WORKSPACE TSMSLV
C TSMSLV
C OUTPUT PARAMETERS : TSMSLV
C ------------------ TSMSLV
C TSMSLV
C DN : STANDARD STEP TSMSLV
C DT : TENSOR STEP TSMSLV
C FLAG : RETURNED CODE WITH THE FOLLOWING MEANING : TSMSLV
C FLAG = 0 : NO SINGULARITY DETECTED WHEN FACTORIZING AJA TSMSLV
C FLAG = 1 : SINGULARITY DETECTED WHEN FACTORIZING AJA TSMSLV
C FROM 1 TO N-P TSMSLV
C FLAG = 2 : SINGULARITY DETECTED WHEN FACTORIZING AJA TSMSLV
C FROM N-P TO N TSMSLV
C IERR : RETURNED CODE WITH THE FOLLOWING MEANING : TSMSLV
C IERR = 0 : NO SINGULARITY DETECTED WHEN FACTORIZING AJATSMSLV
C IERR = 1 : SINGULARITY DETECTED WHEN FACTORIZING AJA TSMSLV
C XPLS : LOCAL MINIMUM OF OPTIMIZATION FUNCTION FCN (USED TSMSLV
C BY UNCMIN) TSMSLV
C FPLS : FUNCTION VALUE AT SOLUTION OF OPTIMIZATION FUNCTION FCN TSMSLV
C (USED IN UNCMIN) TSMSLV
C GPLS : GRADIENT AT SOLUTION XPLS (USED BY UNCMIN) TSMSLV
C CURPOS,PIVOT,PBAR : PIVOT VECTORS TSMSLV
C RESTNS : TENSOR RESIDUAL TSMSLV
C RESNEW : STANDARD RESIDUAL TSMSLV
C ITRMCD : TERMINATION CODE (FOR UNCMIN) TSMSLV
C TSMSLV
C SUBPROGRAMS CALLED: TSMSLV
C TSMSLV
C LEVEL 1 BLAS ... DCOPY,DLOAD TSMSLV
C TENSOLVE ... TSQLFC,QTRNS,TSQRFC,TSQMTS,TSQMUV,TSSQP1 TSMSLV
C TENSOLVE ... TSQ1P1,TSD1SV,TSPRMV,TSQMLV,TSCPSS TSMSLV
C UNCMIN ... DFAUT,OPTIF9 TSMSLV
C TSMSLV
C********************************************************************** TSMSLV
TSMSLV
INTEGER Q,METH,IEXP,NDIGIT,IAGFLG,IAHFLG TSMSLV
DOUBLE PRECISION ROOT,TYPFU,DLT,GRADLT,STEPMX,STEPTL,FPLS TSMSLV
DOUBLE PRECISION ZERO,ONE,TWO TSMSLV
INTRINSIC SQRT TSMSLV
EXTERNAL TSQFCN,TSDFCN,D2FCN TSMSLV
DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ TSMSLV
TSMSLV
IF(N .EQ. 1) THEN TSMSLV
S(2,1) = ONE TSMSLV
S(3,1) = ONE TSMSLV
CURPOS(1) = 1 TSMSLV
CALL DCOPY(M,FC,1,FQ,1) TSMSLV
ELSE TSMSLV
TSMSLV
c perform a QL decomposition of S TSMSLV
TSMSLV
CALL TSQLFC(S,MAXN,N,P,IERR) TSMSLV
TSMSLV
c compute AJA times Q-trans TSMSLV
TSMSLV
CALL TSJQTP(S,MAXM,MAXN,N,M,P,WRK1,FQ,AJA) TSMSLV
TSMSLV
c perform a QR factorization of AJA TSMSLV
TSMSLV
CALL TSQRFC(AJA,MAXM,N,M,1,N-P,IERR,EPSM,WRK1,CURPOS,ZERO1) TSMSLV
TSMSLV
IF(IERR.EQ.1) THEN TSMSLV
Q = N-ZERO1+1 TSMSLV
ELSE TSMSLV
Q = P TSMSLV
ENDIF TSMSLV
CALL TSQMTS(ANLS,AJA,MAXM,M,N,M,P,1,WRK1,ZERO1) TSMSLV
TSMSLV
CALL TSQMUV(AJA,FC,FQ,MAXM,M,1,ZERO1,.FALSE.) TSMSLV
ENDIF TSMSLV
TSMSLV
c minimize the lower m-n+q quadratic equations in p unknowns TSMSLV
c of the tensor model. The minimization is performed analytically TSMSLV
c if p=1,q>1, or p=1,q=1,m>n, or n=1,m>n. Otherwise an unconstrained TSMSLV
c minimization software, UNCMIN, is used. TSMSLV
TSMSLV
IF((P.EQ.1.AND.Q.GT.1).OR.(P.EQ.1 .AND. Q.EQ.1 .AND. M.GT.N) TSMSLV
+ .OR. (N .EQ. 1 .AND. M .GT. N)) THEN TSMSLV
CALL TSSQP1(AJA,ANLS,S,FQ,MAXM,MAXN,M,N,Q,ROOT,RESTNS) TSMSLV
XPLS(1) = ROOT TSMSLV
ELSEIF((M.EQ.N).AND.(P.EQ.1).AND.(Q.EQ.1) .OR. TSMSLV
+ (M.EQ.1.AND.N.EQ.1)) THEN TSMSLV
CALL TSQ1P1(AJA,ANLS,S,FQ,MAXM,MAXN,N,ROOT,RESTNS) TSMSLV
XPLS(1) = ROOT TSMSLV
ELSE TSMSLV
CALL DFAUT(P,TYPXU,TYPFU,METH,IEXP,MSG,NDIGIT,ITNLIM, TSMSLV
+ IAGFLG,IAHFLG,IPR,DLT,GRADLT,STEPMX,STEPTL) TSMSLV
TSMSLV
IAGFLG = 1 TSMSLV
IAHFLG = 0 TSMSLV
IEXP = 0 TSMSLV
METH = 2 TSMSLV
MSG = 9 TSMSLV
TSMSLV
CALL DLOAD (P,ZERO,X,1) TSMSLV
TSMSLV
CALL OPTIF9(SQRN,P,X,TSQFCN,TSDFCN,D2FCN,TYPXU,TYPFU,METH,IEXP,TSMSLV
+ MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,DLT,GRADLT, TSMSLV
+ STEPMX,STEPTL,XPLS,FPLS,GPLS,ITRMCD,A,WRK,AJA, TSMSLV
+ ANLS,S,FQ,WRK1,FQQ,WRK2,WRK3,WRK4,MAXM,MAXN,M,N,Q) TSMSLV
TSMSLV
c compute the tensor residual TSMSLV
TSMSLV
RESTNS = SQRT(TWO*FPLS) TSMSLV
ENDIF TSMSLV
TSMSLV
CALL DCOPY(P,XPLS,1,WRK4(N-P+1),1) TSMSLV
TSMSLV
IF(N .EQ. 1) THEN TSMSLV
DT(1) = WRK4(1) TSMSLV
ELSE TSMSLV
TSMSLV
c compute the first n-p components of the tensor step TSMSLV
TSMSLV
CALL TSD1SV(AJA,S,ANLS,FQ,XPLS,MAXM,MAXN,M,N,P,Q,EPSM, TSMSLV
+ WRK1,FQQ,WRK2,PIVOT,WRK3) TSMSLV
CALL TSPRMV(WRK4,WRK3,CURPOS,N-P,0) TSMSLV
TSMSLV
c premultiply the tensor step by the orthogonal matrix resulting TSMSLV
c from the QL factorization of S TSMSLV
TSMSLV
CALL TSQMLV(MAXN,N,P,S,WRK4,DT,.TRUE.) TSMSLV
ENDIF TSMSLV
TSMSLV
c compute the standard step if needed TSMSLV
TSMSLV
IF(GLOBAL .EQ. 1 .OR. (M .GT. N .AND. GLOBAL .EQ. 0)) THEN TSMSLV
CALL TSCPSS(S,MAXM,MAXN,M,N,P,METHOD,GLOBAL,EPSM,FQ, TSMSLV
+ WRK1,WRK2,WRK3,WRK4,AJA,ANLS,DN,FQQ,PIVOT, TSMSLV
+ CURPOS,PBAR,ZERO1,IERR,RESNEW,FLAG) TSMSLV
ENDIF TSMSLV
TSMSLV
RETURN TSMSLV
END TSMSLV
SUBROUTINE TSNECI(MAXM,MAXN,MAXP,X0,M,N,TYPX,TYPF,ITNLIM, TSNECI
+ JACFLG,GRADTL,STEPTL,FTOL,METHOD,GLOBAL, TSNECI
+ STEPMX,DLT,IPR,WRKUNC,LUNC,WRKNEM,LNEM, TSNECI
+ WRKNEN,LNEN,IWRKN,LIN,TSANJA,TSFVEC,MSG, TSNECI
+ XP,FP,GP,TERMCD) TSNECI
TSNECI
INTEGER MAXM,MAXN,M,N,MAXP,JACFLG,ITNLIM,TERMCD,METHOD TSNECI
INTEGER MSG,GLOBAL,IPR,LUNC,LNEM,LNEN,LIN TSNECI
INTEGER IWRKN(MAXN,LIN) TSNECI
DOUBLE PRECISION STEPTL,GRADTL,FTOL,STEPMX,DLT TSNECI
DOUBLE PRECISION XP(N),FP(M),GP(N),X0(N) TSNECI
DOUBLE PRECISION WRKUNC(MAXP,LUNC) TSNECI
DOUBLE PRECISION WRKNEM(MAXM,LNEM) TSNECI
DOUBLE PRECISION WRKNEN(MAXN,LNEN) TSNECI
DOUBLE PRECISION TYPX(N),TYPF(M) TSNECI
EXTERNAL TSANJA,TSFVEC TSNECI
TSNECI
C TSNECI
C********************************************************************** TSNECI
C THIS ROUTINE PROVIDES A COMPLETE INTERFACE TO THE NONLINEAR EQUATION/ TSNECI
C NONLINEAR LEAST SQUARES PACKAGE. THE USER HAS FULL CONTROL OVER TSNECI
C THE OPTIONS. TSNECI
C********************************************************************** TSNECI
C TSNECI
C SUBPROGRAMS CALLED: TSNECI
C TSNECI
C TENSOLVE ... TSCHKI,TSNESV TSNECI
C TSNECI
C********************************************************************** TSNECI
TSNECI
INTEGER SQRN TSNECI
DOUBLE PRECISION EPSM TSNECI
TSNECI
c check input parameters TSNECI
TSNECI
CALL TSCHKI(MAXM,MAXN,MAXP,M,N,LUNC,LNEM,LNEN,LIN,GRADTL,STEPTL, TSNECI
+ FTOL,ITNLIM,JACFLG,METHOD,GLOBAL,STEPMX,DLT,EPSM, TSNECI
+ MSG,TYPX,TYPF,WRKNEN(1,2),WRKNEM(1,2),SQRN, TSNECI
+ TERMCD,IPR) TSNECI
IF(MSG.LT.0) RETURN TSNECI
TSNECI
c call nonlinear equations/nonlinear least squares solver TSNECI
TSNECI
CALL TSNESV(MAXM,MAXN,MAXP,X0,M,N,TYPX,TYPF,ITNLIM,JACFLG, TSNECI
+ GRADTL,STEPTL,FTOL,METHOD,GLOBAL,STEPMX,DLT,IPR, TSNECI
+ WRKUNC(1,1),WRKUNC(1,2),WRKUNC(1,3),WRKUNC(1,4), TSNECI
+ WRKUNC(1,5),WRKUNC(1,SQRN+5),WRKNEM(1,2),WRKNEM(1,3),TSNECI
+ WRKNEM(1,4),WRKNEM(1,5),WRKNEM(1,6),WRKNEM(1,7), TSNECI
+ WRKNEM(1,8),WRKNEM(1,9),WRKNEM(1,10),WRKNEM(1,11), TSNECI
+ WRKNEM(1,12),WRKNEM(1,SQRN+12),WRKNEM(1,2*SQRN+12), TSNECI
+ WRKNEN(1,2),WRKNEN(1,3),WRKNEN(1,4),WRKNEN(1,5), TSNECI
+ WRKNEN(1,6),WRKNEN(1,7),WRKNEN(1,8),WRKNEN(1,9), TSNECI
+ WRKNEN(1,10),WRKNEN(1,SQRN+10),IWRKN(1,1),IWRKN(1,2),TSNECI
+ IWRKN(1,3),EPSM,SQRN,TSANJA,TSFVEC,MSG,XP,FP,GP, TSNECI
+ TERMCD) TSNECI
TSNECI
RETURN TSNECI
END TSNECI
SUBROUTINE TSNESI(MAXM,MAXN,MAXP,X0,M,N,WRKUNC,LUNC,WRKNEM, TSNESI
+ LNEM,WRKNEN,LNEN,IWRKN,LIN,TSFVEC,MSG,XP, TSNESI
+ FP,GP,TERMCD) TSNESI
TSNESI
INTEGER MAXM,MAXN,M,N,MAXP,JACFLG,ITNLIM,TERMCD,METHOD TSNESI
INTEGER GLOBAL,MSG,IPR,LUNC,LNEM,LNEN,LIN TSNESI
INTEGER IWRKN(MAXN,LIN) TSNESI
DOUBLE PRECISION STEPTL,GRADTL,FTOL,STEPMX,DLT TSNESI
DOUBLE PRECISION XP(N),FP(M),GP(N),X0(N) TSNESI
DOUBLE PRECISION WRKUNC(MAXP,LUNC) TSNESI
DOUBLE PRECISION WRKNEM(MAXM,LNEM) TSNESI
DOUBLE PRECISION WRKNEN(MAXN,LNEN) TSNESI
EXTERNAL TSDUMJ,TSFVEC TSNESI
TSNESI
C********************************************************************** TSNESI
C THIS ROUTINE PROVIDES A SIMPLE INTERFACE TO THE NONLINEAR EQUATION/ TSNESI
C NONLINEAR LEAST SQUARES PROBLEMS PACKAGE. THE USER HAS NO CONTROL TSNESI
C OVER THE PACKAGE OPTIONS. TSNESI
C********************************************************************** TSNESI
C TSNESI
C SUBPROGRAMS CALLED: TSNESI
C TSNESI
C TENSOLVE ... TSDFLT,TSCHKI,TSNESV TSNESI
C TSNESI
C********************************************************************** TSNESI
TSNESI
INTEGER SQRN TSNESI
DOUBLE PRECISION EPSM TSNESI
TSNESI
c set default values for each variable to the nonlinear equations/ TSNESI
c nonlinear least squares solver TSNESI
TSNESI
CALL TSDFLT(M,N,ITNLIM,JACFLG,GRADTL,STEPTL,FTOL,METHOD,GLOBAL, TSNESI
+ STEPMX,DLT,WRKNEN(1,1),WRKNEM(1,1),IPR,MSG) TSNESI
TSNESI
c check input parameters TSNESI
TSNESI
CALL TSCHKI(MAXM,MAXN,MAXP,M,N,LUNC,LNEM,LNEN,LIN,GRADTL,STEPTL, TSNESI
+ FTOL,ITNLIM,JACFLG,METHOD,GLOBAL,STEPMX,DLT,EPSM, TSNESI
+ MSG,WRKNEN(1,1),WRKNEM(1,1),WRKNEN(1,2),WRKNEM(1,2), TSNESI
+ SQRN,TERMCD,IPR) TSNESI
IF(MSG.LT.0) RETURN TSNESI
TSNESI
c call nonlinear equations/nonlinear least squares solver TSNESI
TSNESI
CALL TSNESV(MAXM,MAXN,MAXP,X0,M,N,WRKNEN(1,1),WRKNEM(1,1),ITNLIM,TSNESI
+ JACFLG,GRADTL,STEPTL,FTOL,METHOD,GLOBAL,STEPMX,DLT, TSNESI
+ IPR,WRKUNC(1,1),WRKUNC(1,2),WRKUNC(1,3),WRKUNC(1,4), TSNESI
+ WRKUNC(1,5),WRKUNC(1,SQRN+5), WRKNEM(1,2),WRKNEM(1,3),TSNESI
+ WRKNEM(1,4),WRKNEM(1,5),WRKNEM(1,6),WRKNEM(1,7), TSNESI
+ WRKNEM(1,8),WRKNEM(1,9),WRKNEM(1,10),WRKNEM(1,11), TSNESI
+ WRKNEM(1,12),WRKNEM(1,SQRN+12),WRKNEM(1,2*SQRN+12), TSNESI
+ WRKNEN(1,2),WRKNEN(1,3),WRKNEN(1,4),WRKNEN(1,5), TSNESI
+ WRKNEN(1,6),WRKNEN(1,7),WRKNEN(1,8),WRKNEN(1,9), TSNESI
+ WRKNEN(1,10),WRKNEN(1,SQRN+10),IWRKN(1,1),IWRKN(1,2), TSNESI
+ IWRKN(1,3),EPSM,SQRN,TSDUMJ,TSFVEC,MSG,XP,FP,GP, TSNESI
+ TERMCD) TSNESI
TSNESI
RETURN TSNESI
END TSNESI
SUBROUTINE TSNESV(MAXM,MAXN,MAXP,XC,M,N,TYPX,TYPF,ITNLIM, TSNESV
+ JACFLG,GRADTL,STEPTL,FTOL,METHOD,GLOBAL, TSNESV
+ STEPMX,DLT,IPR,X,TYPXU,XPLS,GPLS,A,WRK,DFN, TSNESV
+ WRK1,WRK2,WRK3,WRK4,WRK5,FQ,FQQ,FC,FHAT, TSNESV
+ ANLS,FV,AJA,DXN,DN,DT,DF,D,GBAR,DBAR,DBARP, TSNESV
+ S,SHAT,CURPOS,PIVOT,PBAR,EPSM,SQRN,TSANJA, TSNESV
+ TSFVEC,MSG,XP,FP,GP,TERMCD) TSNESV
TSNESV
INTEGER MAXM,MAXN,MAXP,M,N,SQRN,TERMCD TSNESV
INTEGER ITNLIM,JACFLG,METHOD,GLOBAL,MSG,IPR TSNESV
INTEGER PBAR(N),CURPOS(N),PIVOT(N) TSNESV
DOUBLE PRECISION GRADTL,STEPTL,FTOL,STEPMX,DLT,FPLS,EPSM TSNESV
DOUBLE PRECISION TYPXU(SQRN),XPLS(SQRN),GPLS(SQRN),A(MAXP,SQRN) TSNESV
DOUBLE PRECISION WRK(MAXP,SQRN),X(SQRN),AJA(MAXM,N),S(MAXN,SQRN)TSNESV
DOUBLE PRECISION ANLS(MAXM,SQRN),SHAT(MAXN,SQRN),FV(MAXM,SQRN) TSNESV
DOUBLE PRECISION XC(N),FC(M),XP(N),FP(M),DN(N),DT(N),DF(N) TSNESV
DOUBLE PRECISION D(N),WRK1(M),WRK2(M),WRK3(M),WRK4(M) TSNESV
DOUBLE PRECISION WRK5(M),FQQ(M),FQ(M),GP(N),FHAT(M) TSNESV
DOUBLE PRECISION GBAR(N),DBAR(N),DBARP(N) TSNESV
DOUBLE PRECISION TYPX(N),TYPF(M),DXN(N),DFN(M) TSNESV
EXTERNAL TSANJA,TSFVEC TSNESV
TSNESV
C********************************************************************** TSNESV
C THIS IS THE DRIVER FOR NONLINEAR EQUATIONS/NONLINEAR LEAST SQUARES TSNESV
C PROBLEMS. TSNESV
C********************************************************************** TSNESV
C TSNESV
C INPUT PARAMETERS : TSNESV
C ----------------- TSNESV
C TSNESV
C MAXM : LEADING DIMENSION OF AJA, ANLS, AND FV TSNESV
C MAXN : LEADING DIMENSION OF S AND SHAT TSNESV
C XC : INITIAL ESTIMATE OF SOLUTION TSNESV
C M,N : DIMENSIONS OF PROBLEM TSNESV
C TYPX : TYPICAL SIZE FOR EACH COMPONENT OF X TSNESV
C TYPF : TYPICAL SIZE FOR EACH COMPONENT OF F TSNESV
C ITNLIM : MAXIMUM NUMBER OF ALLOWABLE ITERATIONS TSNESV
C JACFLG : JACOBIAN FLAG WITH THE FOLLOWING MEANINGS: TSNESV
C JACFLG = 1 IF ANALYTIC JACOBIAN SUPPLIED TSNESV
C JACFLG = 0 IF ANALYTIC JACOBIAN NOT SUPPLIED TSNESV
C GRADTL : TOLERANCE AT WHICH GRADIENT IS CONSIDERED CLOSE ENOUGH TSNESV
C TO ZERO TO TERMINATE ALGORITHM TSNESV
C STEPTL : TOLERANCE AT WHICH SUCCESSIVE ITERATES ARE CONSIDERED TSNESV
C CLOSE ENOUGH TO TERMINATE ALGORITHM TSNESV
C FTOL : TOLERANCE AT WHICH FUNCTION VALUE IS CONSIDERED CLOSE TSNESV
C ENOUGH TO ZERO TSNESV
C METHOD : METHOD TO USE TSNESV
C METHOD = 0 : STANDARD METHOD IS USED TSNESV
C METHOD = 1 : TENSOR METHOD IS USED TSNESV
C GLOBAL : GLOBAL STRATEGY TO USE TSNESV
C GLOBAL = 0 : LINE SEARCH TSNESV
C GLOBAL = 1 : 2-DIMENSIONAL TRUST REGION TSNESV
C STEPMX : MAXIMUM ALLOWABLE STEP SIZE TSNESV
C DLT : TRUST REGION RADIUS TSNESV
C IPR : DEVICE TO WHICH TO SEND OUTPUT TSNESV
C X : ESTIMATE TO A ROOT OF FCN ( USED BY UNCMIN) TSNESV
C TYPXU : TYPICAL SIZE FOR EACH COMPONENT OF X (USED BY UNCMIN) TSNESV
C XPLS : LOCAL MINIMUM OF OPTIMIZATION FUNCTION FCN USED BY TSNESV
C UNCMIN TSNESV
C GPLS : GRADIENT AT SOLUTION XPLS (USED BY UNCMIN) TSNESV
C A : WORKSPACE FOR HESSIAN (OR ESTIMATE) (USED BY UNCMIN) TSNESV
C WRK : WORKSPACE (USED BY UNCMIN) TSNESV
C WRK1,WRK2,WRK3,WRK4,WRK5,FQ,FQQ: WORKSPACE TSNESV
C FC : FUNCTION VALUE AT CURRENT ITERATE TSNESV
C FHAT : WORKSPACE TSNESV
C DFN : DIAGONAL SCALING MATRIX FOR F TSNESV
C ANLS : TENSOR TERM MATRIX TSNESV
C FV : WORKSPACE USED TO STORE PAST FUNCTION VALUES TSNESV
C AJA : JACOBIAN MATRIX TSNESV
C DN : STANDARD STEP TSNESV
C DT : TENSOR STEP TSNESV
C DF,D,GBAR,DBAR,DBARP : WORKSPACE TSNESV
C DXN : DIAGONAL SCALING MATRIX FOR X TSNESV
C S : MATRIX OF PREVIOUS DIRECTIONS TSNESV
C SHAT : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSNESV
C CURPOS,PIVOT,PBAR : PIVOT VECTORS TSNESV
C SQRN : MAXIMUM COLUMN DIMENSION OF ANLS, S, AND SHAT TSNESV
C EPSM : MACHINE PRECISION TSNESV
C TSANJA : (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE JACOBIAN. TSNESV
C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE TSNESV
C TSFVEC : NAME OF SUBROUTINE TO EVALUATE FUNCTION TSNESV
C TSNESV
C TSNESV
C INPUT-OUTPUT PARAMETERS : TSNESV
C ------------------------ TSNESV
C TSNESV
C MSG : MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT TSNESV
C TSNESV
C TSNESV
C OUTPUT PARAMETERS : TSNESV
C ----------------- TSNESV
C TSNESV
C XP : SOLUTION TO THE SYSTEM OF NONLINEAR EQUATIONS TSNESV
C FP : FUNCTION VALUE AT THE SOLUTION TSNESV
C GP : GRADIENT AT THE SOLUTION TSNESV
C TERMCD : TERMINATION CODE TSNESV
C TSNESV
C SUBPROGRAMS CALLED: TSNESV
C TSNESV
C LEVEL 1 BLAS ... DCOPY,DNRM2 TSNESV
C LEVEL 2 BLAS ... DGEMV TSNESV
C TENSOLVE ... TSSCLX,TSFSCL,TSSCLJ,TSCHKJ,TSNSTP,TSSSTP, TSNESV
C TENSOLVE ... TSLSCH,TS2DTR,TSRSLT,TSMGSA,TSFRMT,TSMSLV, TSNESV
C TENSOLVE ... TSSLCT,TSMDLS,TSUPSF TSNESV
C TSNESV
C********************************************************************* TSNESV
TSNESV
INTEGER P,ITN,I,FLAG,RETCD,ZERO1,IERR,ITRMCD,ICSCMX TSNESV
DOUBLE PRECISION FNORM,RESTNS,RESNEW TSNESV
DOUBLE PRECISION ZERO,HALF,ONE TSNESV
DOUBLE PRECISION DNRM2 TSNESV
LOGICAL NWTAKE,MXTAKE TSNESV
DATA ZERO,HALF,ONE/0.0D0,0.50D0,1.0D0/ TSNESV
TSNESV
c----------------- TSNESV
c initialization TSNESV
c----------------- TSNESV
TSNESV
ITN = 0 TSNESV
NWTAKE = .TRUE. TSNESV
TSNESV
CALL TSSCLX(XC,DXN,N) TSNESV
TSNESV
IF(MOD(MSG/8,2).NE.1) THEN TSNESV
WRITE(IPR,896) TSNESV
WRITE(IPR,900) (TYPX(I),I = 1,N) TSNESV
WRITE(IPR,897) TSNESV
WRITE(IPR,900) (DXN(I),I = 1,N) TSNESV
WRITE(IPR,898) TSNESV
WRITE(IPR,900) (TYPF(I),I = 1,M) TSNESV
WRITE(IPR,899) TSNESV
WRITE(IPR,900) (DFN(I),I = 1,M) TSNESV
WRITE(IPR,901) JACFLG TSNESV
WRITE(IPR,902) METHOD TSNESV
WRITE(IPR,903) GLOBAL TSNESV
WRITE(IPR,904) ITNLIM TSNESV
WRITE(IPR,905) EPSM TSNESV
WRITE(IPR,906) STEPTL TSNESV
WRITE(IPR,907) GRADTL TSNESV
WRITE(IPR,908) FTOL TSNESV
WRITE(IPR,909) STEPMX TSNESV
WRITE(IPR,910) DLT TSNESV
ENDIF TSNESV
TSNESV
c evaluate analytic or finite difference Jacobian and check analytic TSNESV
c Jacobian, if requested TSNESV
TSNESV
CALL TSFSCL(XC,DXN,DFN,M,N,TSFVEC,FC) TSNESV
CALL TSSCLJ(XC,DXN,TYPX,FC,DFN,FHAT,MAXM,M,N, TSNESV
+ EPSM,JACFLG,TSFVEC,TSANJA,AJA) TSNESV
IF(JACFLG.EQ.1) THEN TSNESV
IF(MOD(MSG/2,2).EQ.0) THEN TSNESV
CALL TSCHKJ(AJA,XC,FC,MAXM,M,N,EPSM,DFN,DXN,TYPX, TSNESV
+ IPR,FHAT,WRK1,TSFVEC,MSG) TSNESV
IF(MSG.LT.0) RETURN TSNESV
ENDIF TSNESV
ENDIF TSNESV
TSNESV
c compute the gradient at the current iterate XC TSNESV
TSNESV
CALL DGEMV('T',M,N,ONE,AJA,MAXM,FC,1,ZERO,GP,1) TSNESV
TSNESV
c compute the residual of FC TSNESV
TSNESV
FNORM = HALF*DNRM2(M,FC,1)**2 TSNESV
TSNESV
c check stopping criteria for input XC TSNESV
TSNESV
CALL TSNSTP(GP,XC,FC,XC,STEPTL,GRADTL,RETCD,FTOL,ITN, TSNESV
+ ITNLIM,ICSCMX,MXTAKE,M,N,MSG,IPR,FNORM,TERMCD) TSNESV
TSNESV
IF(TERMCD.GT.0) THEN TSNESV
FPLS = FNORM TSNESV
GO TO 120 TSNESV
ENDIF TSNESV
TSNESV
c--------------- TSNESV
c iteration 1 TSNESV
c--------------- TSNESV
TSNESV
ITN = 1 TSNESV
TSNESV
c compute the standard step TSNESV
TSNESV
CALL DCOPY(M,FC,1,FHAT,1) TSNESV
TSNESV
CALL TSSSTP(AJA,FHAT,M,N,MAXM,EPSM,GLOBAL,WRK1,WRK2,WRK3, TSNESV
+ DN,FQQ,PIVOT,PBAR,IERR) TSNESV
TSNESV
c choose next iterate XP by a global strategy TSNESV
TSNESV
IF(GLOBAL.EQ.0) THEN TSNESV
CALL TSLSCH(M,N,XC,FC,DN,GP,STEPTL,DXN,DFN,TSFVEC, TSNESV
+ MXTAKE,STEPMX,XP,FP,FNORM,FPLS,RETCD) TSNESV
ELSE TSNESV
CALL TS2DTR(AJA,SHAT,ANLS,DN,GP,GBAR,XC,FC,METHOD, TSNESV
+ NWTAKE,STEPMX,STEPTL,EPSM,MXTAKE,DLT, TSNESV
+ FQQ,MAXM,MAXN,M,N,P,CURPOS,PIVOT,PBAR, TSNESV
+ ITN,IERR,FLAG,DXN,DFN,TSFVEC,DBAR,DBARP, TSNESV
+ D,FHAT,WRK1,WRK2,WRK3,WRK4,WRK5, TSNESV
+ XPLS,GPLS,FNORM,XP,FP,FPLS,RETCD) TSNESV
ENDIF TSNESV
TSNESV
IF(MOD(MSG/8,2).EQ.0) CALL TSRSLT(N,XC,FNORM,GP,0,IPR) TSNESV
TSNESV
c evaluate the Jacobian at the new iterate XP TSNESV
TSNESV
CALL TSSCLJ(XP,DXN,TYPX,FP,DFN,FHAT,MAXM,M,N,EPSM,JACFLG, TSNESV
+ TSFVEC,TSANJA,AJA) TSNESV
TSNESV
c compute the gradient at the new iterate XP TSNESV
TSNESV
CALL DGEMV('T',M,N,ONE,AJA,MAXM,FP,1,ZERO,GP,1) TSNESV
TSNESV
c check stopping criteria for the new iterate XP TSNESV
TSNESV
CALL TSNSTP(GP,XP,FP,XC,STEPTL,GRADTL,RETCD,FTOL,ITN, TSNESV
+ ITNLIM,ICSCMX,MXTAKE,M,N,MSG,IPR,FPLS,TERMCD) TSNESV
TSNESV
IF(TERMCD.GT.0) GO TO 120 TSNESV
IF(MOD(MSG/16,2).EQ.1) CALL TSRSLT(N,XP,FPLS,GP,ITN,IPR) TSNESV
TSNESV
c update S and FV TSNESV
TSNESV
DO 40 I = 1,N TSNESV
S(I,1) = XC(I)-XP(I) TSNESV
40 CONTINUE TSNESV
CALL DCOPY(M,FC,1,FV(1,1),1) TSNESV
TSNESV
c update XC and FC TSNESV
TSNESV
CALL DCOPY(N,XP,1,XC,1) TSNESV
CALL DCOPY(M,FP,1,FC,1) TSNESV
FNORM = FPLS TSNESV
TSNESV
c--------------- TSNESV
c iteration > 1 TSNESV
c--------------- TSNESV
TSNESV
80 ITN = ITN+1 TSNESV
TSNESV
c if the standard method is selected then compute the standard step TSNESV
TSNESV
IF(METHOD.EQ.0) THEN TSNESV
CALL DCOPY(M,FC,1,FHAT,1) TSNESV
CALL TSSSTP(AJA,FHAT,M,N,MAXM,EPSM,GLOBAL,WRK1,WRK2, TSNESV
+ WRK3,DF,FQQ,PIVOT,PBAR,IERR) TSNESV
ENDIF TSNESV
TSNESV
c if the tensor method is selected then form the tensor model TSNESV
TSNESV
IF(METHOD.EQ.1) THEN TSNESV
TSNESV
c select the past linearly independent directions TSNESV
TSNESV
CALL TSMGSA(S,MAXN,N,SQRN,ITN,SHAT,P,CURPOS) TSNESV
TSNESV
c form the tensor term TSNESV
TSNESV
CALL TSFRMT(SHAT,S,AJA,FV,FC,MAXM,MAXN,MAXP,M,N,P, TSNESV
+ CURPOS,A,X,XPLS,GPLS,ANLS) TSNESV
TSNESV
c solve the tensor model for the tensor step DT and compute DN TSNESV
c as a by-product if the global strategy selected is the TSNESV
c two-dimensional trust region or M > N TSNESV
TSNESV
CALL TSMSLV(AJA,SHAT,ANLS,FC,P,MAXM,MAXN,SQRN,M,N, TSNESV
+ EPSM,METHOD,GLOBAL,WRK1,WRK2,WRK3,WRK4, TSNESV
+ X,TYPXU,XPLS,GPLS,A,WRK,CURPOS,PBAR,PIVOT, TSNESV
+ FQ,FQQ,DN,DT,RESTNS,RESNEW,ITRMCD,FLAG, TSNESV
+ ZERO1,IERR) TSNESV
TSNESV
c decide which step to use (DN or DT) TSNESV
TSNESV
IF(GLOBAL.EQ.1 .OR. (M.GT.N .AND. GLOBAL.EQ.0)) THEN TSNESV
CALL TSSLCT(RESTNS,RESNEW,ITRMCD,FC,M,N,DN,DT,GP, TSNESV
+ DF,NWTAKE) TSNESV
ENDIF TSNESV
TSNESV
ENDIF TSNESV
TSNESV
c choose the next iterate XP by a global strategy TSNESV
TSNESV
IF(GLOBAL.EQ.0) THEN TSNESV
IF(METHOD.EQ.0) THEN TSNESV
CALL TSLSCH(M,N,XC,FC,DF,GP,STEPTL,DXN,DFN,TSFVEC, TSNESV
+ MXTAKE,STEPMX,XP,FP,FNORM,FPLS,RETCD) TSNESV
ELSEIF(M.EQ.N) THEN TSNESV
CALL TSMDLS(AJA,SHAT,ANLS,XC,FC,M,N,MAXM,MAXN,P,DT,GP, TSNESV
+ DXN,DFN,TSFVEC,METHOD,STEPTL,GLOBAL,STEPMX, TSNESV
+ EPSM,FQ,WRK1,WRK2,WRK3,WRK4,DN,FQQ,PIVOT, TSNESV
+ CURPOS,PBAR,MXTAKE,XP,FP,FNORM,FPLS, TSNESV
+ ZERO1,RETCD,IERR) TSNESV
ELSE TSNESV
CALL TSLSCH(M,N,XC,FC,DF,GP,STEPTL,DXN,DFN,TSFVEC, TSNESV
+ MXTAKE,STEPMX,XP,FP,FNORM,FPLS,RETCD) TSNESV
ENDIF TSNESV
ELSE TSNESV
CALL TS2DTR(AJA,SHAT,ANLS,DF,GP,GBAR,XC,FC, TSNESV
+ METHOD,NWTAKE,STEPMX,STEPTL,EPSM,MXTAKE, TSNESV
+ DLT,FQQ,MAXM,MAXN,M,N,P,CURPOS,PIVOT, TSNESV
+ PBAR,ITN,IERR,FLAG,DXN,DFN,TSFVEC,DBAR, TSNESV
+ DBARP,D,FHAT,WRK1,WRK2,WRK3,WRK4,WRK5, TSNESV
+ XPLS,GPLS,FNORM,XP,FP,FPLS,RETCD) TSNESV
ENDIF TSNESV
TSNESV
c evaluate the Jacobian at the new iterate XP TSNESV
TSNESV
CALL TSSCLJ(XP,DXN,TYPX,FP,DFN,FHAT,MAXM,M,N,EPSM, TSNESV
+ JACFLG,TSFVEC,TSANJA,AJA) TSNESV
TSNESV
c evaluate the gradient at the new iterate XP TSNESV
TSNESV
CALL DGEMV('T',M,N,ONE,AJA,MAXM,FP,1,ZERO,GP,1) TSNESV
TSNESV
c check stopping criteria for the new iterate XP TSNESV
TSNESV
CALL TSNSTP(GP,XP,FP,XC,STEPTL,GRADTL,RETCD,FTOL,ITN, TSNESV
+ ITNLIM,ICSCMX,MXTAKE,M,N,MSG,IPR,FPLS,TERMCD) TSNESV
TSNESV
IF(TERMCD.GT.0) GO TO 120 TSNESV
IF(MOD(MSG/16,2).EQ.1) CALL TSRSLT(N,XP,FPLS,GP,ITN,IPR) TSNESV
TSNESV
c if tensor method is selected then update the matrices S and FV TSNESV
TSNESV
IF(METHOD.EQ.1) THEN TSNESV
CALL TSUPSF(FC,XC,XP,SQRN,ITN,MAXM,MAXN,M,N,WRK1,S,FV) TSNESV
ENDIF TSNESV
TSNESV
c update XC, FC, and FNORM TSNESV
TSNESV
CALL DCOPY(N,XP,1,XC,1) TSNESV
CALL DCOPY(M,FP,1,FC,1) TSNESV
FNORM = FPLS TSNESV
GO TO 80 TSNESV
TSNESV
c termination TSNESV
TSNESV
120 IF(MOD(MSG/8,2).EQ.0) THEN TSNESV
IF(ITN.NE.0) THEN TSNESV
CALL TSRSLT(N,XP,FPLS,GP,ITN,IPR) TSNESV
ELSE TSNESV
FPLS = HALF*DNRM2(M,FC,1)**2 TSNESV
CALL TSRSLT(N,XC,FPLS,GP,ITN,IPR) TSNESV
ENDIF TSNESV
ENDIF TSNESV
IF (ITN.EQ.0) THEN TSNESV
CALL DCOPY(N,XC,1,XP,1) TSNESV
CALL DCOPY(M,FC,1,FP,1) TSNESV
FNORM = FPLS TSNESV
ENDIF TSNESV
CALL TSUNSX(XP,DXN,N) TSNESV
TSNESV
896 FORMAT(' TSNESV TYPICAL X') TSNESV
897 FORMAT(' TSNESV DIAGONAL SCALING MATRIX FOR X') TSNESV
898 FORMAT(' TSNESV TYPICAL F') TSNESV
899 FORMAT(' TSNESV DIAGONAL SCALING MATRIX FOR F') TSNESV
900 FORMAT(100(' TSNESV ',3(D20.13,3X),/)) TSNESV
901 FORMAT(' TSNESV JACOBIAN FLAG = ',I1) TSNESV
902 FORMAT(' TSNESV METHOD USED = ',I1) TSNESV
903 FORMAT(' TSNESV GLOBAL STRATEGY = ',I1) TSNESV
904 FORMAT(' TSNESV ITERATION LIMIT =',I5) TSNESV
905 FORMAT(' TSNESV MACHINE EPSILON =',D20.13) TSNESV
906 FORMAT(' TSNESV STEP TOLERANCE =',D20.13) TSNESV
907 FORMAT(' TSNESV GRADIENT TOLERANCE =',D20.13) TSNESV
908 FORMAT(' TSNESV FUNCTION TOLERANCE =',D20.13) TSNESV
909 FORMAT(' TSNESV MAXIMUM STEP SIZE =',D20.13) TSNESV
910 FORMAT(' TSNESV TRUST REG RADIUS =',D20.13) TSNESV
END TSNESV
SUBROUTINE TSNSTP(G,XPLUS,FPLUS,XC,STEPTL,GRADTL,RETCD,FTOL,ITN, TSNSTP
+ ITNLIM,ICSCMX,MXTAKE,M,N,MSG,IPR,FNORM,TERMCD) TSNSTP
TSNSTP
INTEGER M,N,ITN,ITNLIM,MSG,IPR,TERMCD,RETCD,ICSCMX TSNSTP
DOUBLE PRECISION STEPTL,GRADTL,FTOL,FNORM TSNSTP
DOUBLE PRECISION XPLUS(N),FPLUS(M),G(N),XC(N) TSNSTP
LOGICAL MXTAKE TSNSTP
TSNSTP
C********************************************************************** TSNSTP
C THIS ROUTINE DECIDES WHETHER TO TERMINATE THE NONLINEAR ALGORITHM. TSNSTP
C********************************************************************** TSNSTP
C TSNSTP
C INPUT PARAMETERS : TSNSTP
C ------------------ TSNSTP
C TSNSTP
C G : GRADIENT AT XC TSNSTP
C XPLUS : NEW ITERATE TSNSTP
C FPLUS : FUNCTION VALUE AT XPLUS TSNSTP
C XC : CURRENT ITERATE TSNSTP
C STEPTL: STEP TOLERANCE TSNSTP
C GRADTL: GRADIENT TOLERANCE TSNSTP
C RETCD : RETURN CODE WITH THE FOLLOWING MEANINGS : TSNSTP
C RETCD = 0 : SUCCESSFUL GLOBAL STRATEGY TSNSTP
C RETCD = 1 : UNSUCCESSFUL GLOBAL STRATEGY TSNSTP
C FTOL : FUNCTION TOLERANCE TSNSTP
C ITN : ITERATION NUMBER TSNSTP
C ITNLIM: ITERATION NUMBER LIMIT TSNSTP
C ICSCMX: NUMBER CONSECUTIVE STEPS .GE. STEPMX TSNSTP
C MXTAKE: BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH TSNSTP
C M : DIMENSION OF FPLUS TSNSTP
C N : DIMENSION OF G, XC, AND XPLUS TSNSTP
C MSG : MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT TSNSTP
C IPR : DEVICE TO WHICH TO SEND OUTPUT TSNSTP
C TSNSTP
C TSNSTP
C OUTPUT PARAMETERS : TSNSTP
C ------------------ TSNSTP
C TSNSTP
C TERMCD: RETURN CODE WITH THE FOLLOWING MEANINGS : TSNSTP
C TERMCD = 0 NO TERMINATION CRITERION SATISFIED TSNSTP
C TSNSTP
C TERMCD > 0 : SOME TERMINATION CRITERION SATISFIED TSNSTP
C TERMCD = 1 : NORM OF SCALED FUNCTION VALUE IS LESS THAN TSNSTP
C FTOL TSNSTP
C TERMCD = 2 : GRADIENT TOLERANCE REACHED TSNSTP
C TERMCD = 3 : SCALED DISTANCE BETWEEN LAST TWO STEPS TSNSTP
C LESS THAN STEPTL TSNSTP
C TERMCD = 4 : UNSUCCESSFUL GLOBAL STRATEGY TSNSTP
C TERMCD = 5 : ITERATION LIMIT EXCEEDED TSNSTP
C TERMCD = 6 : FIVE CONSECUTIVE STEPS OF LENGTH STEPMX TSNSTP
C HAVE BEEN TAKEN TSNSTP
C TSNSTP
C SUBPROGRAMS CALLED: TSNSTP
C TSNSTP
C LEVEL 1 BLAS ... IDAMAX TSNSTP
C TSNSTP
C********************************************************************** TSNSTP
TSNSTP
INTEGER I TSNSTP
DOUBLE PRECISION RES,D,RGX,RELGRD,RSX,RELSTP,ZERO,ONE TSNSTP
INTEGER IDAMAX TSNSTP
INTRINSIC ABS,MAX TSNSTP
DATA ZERO,ONE/0.0D0,1.0D0/ TSNSTP
TSNSTP
c check whether scaled function is within tolerance TSNSTP
TSNSTP
RES = ABS(FPLUS(IDAMAX(M,FPLUS,1))) TSNSTP
IF(RES.LE.FTOL) THEN TSNSTP
TERMCD = 1 TSNSTP
IF(MOD(MSG/8,2).EQ.0) THEN TSNSTP
WRITE(IPR,701) TSNSTP
ENDIF TSNSTP
RETURN TSNSTP
ENDIF TSNSTP
TSNSTP
c check whether scaled gradient is within tolerance TSNSTP
TSNSTP
D = ONE/MAX(FNORM, DBLE(N/2)) TSNSTP
RGX = ZERO TSNSTP
DO 200 I = 1,N TSNSTP
RELGRD = ABS(G(I)) * MAX(ABS(XPLUS(I)), ONE)*D TSNSTP
RGX = MAX(RGX,RELGRD) TSNSTP
200 CONTINUE TSNSTP
IF(RGX.LE.GRADTL) THEN TSNSTP
TERMCD = 2 TSNSTP
IF(MOD(MSG/8,2).EQ.0) THEN TSNSTP
WRITE(IPR,702) TSNSTP
ENDIF TSNSTP
RETURN TSNSTP
ENDIF TSNSTP
TSNSTP
IF(ITN.EQ.0) RETURN TSNSTP
TSNSTP
IF(RETCD.EQ.1) THEN TSNSTP
TERMCD = 4 TSNSTP
IF(MOD(MSG/8,2).EQ.0) THEN TSNSTP
WRITE(IPR,703) TSNSTP
ENDIF TSNSTP
RETURN TSNSTP
ENDIF TSNSTP
TSNSTP
c check whether relative step length is within tolerance TSNSTP
TSNSTP
RSX = ZERO TSNSTP
DO 300 I = 1,N TSNSTP
RELSTP = ABS(XPLUS(I) - XC(I))/MAX(XPLUS(I), ONE) TSNSTP
RSX = MAX(RSX, RELSTP) TSNSTP
300 CONTINUE TSNSTP
IF(RSX.LE.STEPTL) THEN TSNSTP
TERMCD = 3 TSNSTP
IF(MOD(MSG/8,2).EQ.0) THEN TSNSTP
WRITE(IPR,704) TSNSTP
ENDIF TSNSTP
RETURN TSNSTP
ENDIF TSNSTP
TSNSTP
c check iteration limit TSNSTP
TSNSTP
IF(ITN.GE.ITNLIM) THEN TSNSTP
TERMCD = 5 TSNSTP
IF(MOD(MSG/8,2).EQ.0) THEN TSNSTP
WRITE(IPR,705) TSNSTP
ENDIF TSNSTP
ENDIF TSNSTP
TSNSTP
c check number of consecutive steps .ge. stepmx TSNSTP
TSNSTP
IF(MXTAKE) THEN TSNSTP
ICSCMX = ICSCMX+1 TSNSTP
IF(ICSCMX.GE.5) THEN TSNSTP
TERMCD = 6 TSNSTP
IF(MOD(MSG/8,2).EQ.0) THEN TSNSTP
WRITE(IPR,706) TSNSTP
ENDIF TSNSTP
ENDIF
ELSE
ICSCMX=0
ENDIF
TSNSTP
701 FORMAT(/,' TSNSTP FUNCTION VALUE CLOSE TO ZERO') TSNSTP
702 FORMAT(/,' TSNSTP RELATIVE GRADIENT CLOSE TO ZERO') TSNSTP
703 FORMAT(/,' TSNSTP LAST GLOBAL STEP FAILED TO LOCATE A',/ TSNSTP
+ ' TSNSTP POINT LOWER THAN THE CURRENT ITERATE') TSNSTP
704 FORMAT(/,' TSNSTP SUCCESSIVE ITERATES WITHIN TOLERANCE',/ TSNSTP
+ ' TSNSTP CURRENT ITERATE IS PROBABLY SOLUTION') TSNSTP
705 FORMAT(/,' TSNSTP ITERATION LIMIT EXCEEDED',/ TSNSTP
+ ' TSNSTP ALGORITHM FAILED') TSNSTP
706 FORMAT(/,' TSNSTP MAXIMUM STEP SIZE EXCEEDED 5', TSNSTP
+ ' CONSECUTIVE TIMES',/ TSNSTP
+ ' TSNSTP EITHER THE FUNCTION IS UNBOUNDED BELOW',/ TSNSTP
+ ' TSNSTP BECOMES ASYMPTOTIC TO A FINITE VALUE',/ TSNSTP
+ ' TSNSTP FROM ABOVE IN SOME DIRECTION',/ TSNSTP
+ ' TSNSTP OR STEPMX IS TOO SMALL') TSNSTP
TSNSTP
RETURN TSNSTP
END TSNSTP
SUBROUTINE TSPRMV(X,Y,PIVOT,N,JOB) TSPRMV
TSPRMV
INTEGER N,JOB TSPRMV
INTEGER PIVOT(N) TSPRMV
DOUBLE PRECISION X(N),Y(N) TSPRMV
TSPRMV
C********************************************************************** TSPRMV
C THIS SUBROUTINE PERFORMS A VECTOR PERMUTATION. TSPRMV
C********************************************************************** TSPRMV
C TSPRMV
C INPUT PARAMETERS : TSPRMV
C ----------------- TSPRMV
C TSPRMV
C Y : VECTOR TO TSPRMV TSPRMV
C PIVOT : PIVOT VECTOR TSPRMV
C N : DIMENSION OF THE VECTORS Y AND PIVOT TSPRMV
C TSPRMV
C OUTPUT PARAMETERS : TSPRMV
C ------------------- TSPRMV
C TSPRMV
C X : PIVOTED VECTOR TSPRMV
C TSPRMV
C********************************************************************** TSPRMV
TSPRMV
INTEGER I TSPRMV
TSPRMV
IF(JOB .EQ. 0) THEN TSPRMV
TSPRMV
c permute Y TSPRMV
TSPRMV
DO 20 I = 1,N TSPRMV
X(PIVOT(I)) = Y(I) TSPRMV
20 CONTINUE TSPRMV
ELSE TSPRMV
TSPRMV
c reverse permute of Y TSPRMV
TSPRMV
DO 30 I = 1,N TSPRMV
X(I) = Y(PIVOT(I)) TSPRMV
30 CONTINUE TSPRMV
TSPRMV
ENDIF TSPRMV
TSPRMV
RETURN TSPRMV
END TSPRMV
SUBROUTINE TSRSLT(N,XP,FVAL,GP,ITN,IPR) TSRSLT
TSRSLT
INTEGER N,ITN,IPR TSRSLT
DOUBLE PRECISION FVAL TSRSLT
DOUBLE PRECISION XP(N),GP(N) TSRSLT
TSRSLT
C********************************************************************** TSRSLT
C THIS ROUTINE PRINTS INFORMATION. TSRSLT
C********************************************************************** TSRSLT
C TSRSLT
C INPUT PARAMETERS : TSRSLT
C ----------------- TSRSLT
C TSRSLT
C M,N : DIMENSIONS OF PROBLEM TSRSLT
C XP : NEXT ITERATE TSRSLT
C FVAL : SUM OF SQUARES OF F(XP) TSRSLT
C GP : GRADIENT AT XP TSRSLT
C ITN : ITERATION NUMBER TSRSLT
C IPR : DEVICE TO WHICH TO SEND OUTPUT TSRSLT
C TSRSLT
C********************************************************************** TSRSLT
TSRSLT
INTEGER I TSRSLT
TSRSLT
WRITE(IPR,801) ITN TSRSLT
WRITE(IPR,802) TSRSLT
WRITE(IPR,803) (XP(I),I = 1,N) TSRSLT
WRITE(IPR,804) TSRSLT
WRITE(IPR,805) FVAL TSRSLT
WRITE(IPR,806) TSRSLT
WRITE(IPR,807) (GP(I),I = 1,N) TSRSLT
TSRSLT
801 FORMAT(/,' TSRSLT ITERATION K =',I5) TSRSLT
802 FORMAT(' TSRSLT X(K)') TSRSLT
803 FORMAT(100(' TSRSLT ',3(D20.13,3X),/)) TSRSLT
804 FORMAT(' TSRSLT FUNCTION AT X(K)') TSRSLT
805 FORMAT(' TSRSLT ',D20.13) TSRSLT
806 FORMAT(' TSRSLT GRADIENT AT X(K)') TSRSLT
807 FORMAT(100(' TSRSLT ',3(D20.13,3X),/)) TSRSLT
TSRSLT
RETURN TSRSLT
END TSRSLT
SUBROUTINE TSQ1P1(AJA,ANLS,S,F,MAXM,MAXN,N,ROOT,RESTNS) TSQ1P1
TSQ1P1
INTEGER MAXM,MAXN,N TSQ1P1
DOUBLE PRECISION ROOT,RESTNS TSQ1P1
DOUBLE PRECISION AJA(MAXM,N),S(MAXN,*),F(N),ANLS(MAXM,*) TSQ1P1
TSQ1P1
C********************************************************************** TSQ1P1
C THIS ROUTINE SOLVES THE LOWER M-N+Q QUADRATIC EQUATIONS IN P UNKNOWNS TSQ1P1
C OF THE TENSOR MODEL WHEN Q = 1 AND P = 1. TSQ1P1
C********************************************************************** TSQ1P1
C TSQ1P1
C INPUT PARAMETERS : TSQ1P1
C ----------------- TSQ1P1
C TSQ1P1
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSQ1P1
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSQ1P1
C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSQ1P1
C F : FUNCTION VALUE AT CURRENT ITERATE MULTIPIED BY AN TSQ1P1
C ORTHOGONAL MATRIX TSQ1P1
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSQ1P1
C MAXN : LEADING DIMENSION OF S TSQ1P1
C N : COLUMN DIMENSION OF AJA TSQ1P1
C TSQ1P1
C OUTPUT PARAMETERS : TSQ1P1
C ----------------- TSQ1P1
C TSQ1P1
C ROOT : SOLUTION TO THE SYSTEM TSQ1P1
C RESTNS : TENSOR RESIDUAL TSQ1P1
C TSQ1P1
C********************************************************************** TSQ1P1
TSQ1P1
DOUBLE PRECISION DELTA,T1,T2,ZERO,HALF,ONE,TWO TSQ1P1
INTRINSIC ABS,SQRT TSQ1P1
DATA ZERO,HALF,ONE,TWO/0.0D0,0.50D0,1.0D0,2.0D0/ TSQ1P1
TSQ1P1
c find the roots of the equation: TSQ1P1
c F(N) + AJA(N,N)*D + 0.5*ANLS(N,1)*(S(N+2,1)*D)**2 TSQ1P1
TSQ1P1
T1 = AJA(N,N) TSQ1P1
T2 = ANLS(N,1) * S(N+2,1)**2 TSQ1P1
IF(ANLS(N,1).EQ.ZERO) THEN TSQ1P1
ROOT = -F(N)/T1 TSQ1P1
ELSE TSQ1P1
DELTA = T1**2 - TWO*F(N)*T2 TSQ1P1
IF(DELTA.GE.ZERO) THEN TSQ1P1
ROOT = (-T1+SIGN(ONE,T1) * SQRT(DELTA))/T2 TSQ1P1
ELSE TSQ1P1
ROOT = -T1/T2 TSQ1P1
ENDIF TSQ1P1
ENDIF TSQ1P1
TSQ1P1
c compute tensor residual TSQ1P1
TSQ1P1
RESTNS = ABS(F(N)+AJA(N,N)*ROOT+HALF*ANLS(N,1)*(S(N+2,1)**2)* TSQ1P1
+ (ROOT**2)) TSQ1P1
RETURN TSQ1P1
END TSQ1P1
SUBROUTINE TSQFCN(P,X,SUM,AJA,ANLS,S,F,WRK1,WRK2,WRK3, TSQFCN
+ WRK4,WRK5,MAXM,MAXN,M,N,Q) TSQFCN
TSQFCN
INTEGER MAXM,MAXN,M,N,P,Q TSQFCN
DOUBLE PRECISION X(P),AJA(MAXM,N),ANLS(MAXM,P),S(MAXN,P) TSQFCN
DOUBLE PRECISION F(M),WRK1(M),WRK2(P),WRK3(P),WRK4(M),WRK5(M) TSQFCN
TSQFCN
C********************************************************************* TSQFCN
C THIS ROUTINE IS USED TO EVALUATE THE RESIDUAL OF THE LAST M-N+P TSQFCN
C QUADRATIC EQUATIONS IN P UNKNOWNS OF THE TENSOR MODEL. NOTE THAT TSQFCN
C THIS ROUTINE IS CALLED BY UNCMIN TO SOLVE THE NONLINEAR LEAST SQUARES TSQFCN
C PART OF THE TENSOR MODEL. TSQFCN
C********************************************************************* TSQFCN
C TSQFCN
C INPUT PARAMETERS : TSQFCN
C ----------------- TSQFCN
C TSQFCN
C P : DIMENSION OF THE PROBLEM SOLVED BY UNCMIN TSQFCN
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSQFCN
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSQFCN
C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSQFCN
C F : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY AN TSQFCN
C ORTHOGONAL MATRIX TSQFCN
C WRK1,WRK2,WRK3,WRK4,WRK5 : WORKING VECTORS TSQFCN
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSQFCN
C MAXN : LEADING DIMENSION OF S TSQFCN
C M,N : DIMENSION OF PROBLEM TSQFCN
C Q : NUMERICAL RANK OF JACOBIAN : TSQFCN
C Q > P : JACOBIAN IS SINGULAR TSQFCN
C Q = P : OTHERWISE TSQFCN
C TSQFCN
C INPUT-OUTPUT PARAMETERS : TSQFCN
C ----------------------- TSQFCN
C TSQFCN
C X : NULL VECTOR ON ENTRY AND APPROXIMATION OF THE SOLUTION TSQFCN
C TO THE SYSTEM OF M-N+Q QUADRATIC EQUATIONS IN P UNKNOWNS TSQFCN
C OF THE TENSOR MODEL ON EXIT TSQFCN
C TSQFCN
C OUTPUT PARAMETERS : TSQFCN
C ----------------- TSQFCN
C TSQFCN
C SUM : RESIDUAL OF THE LAST M-N+P QUADRATIC EQUATIONS IN P TSQFCN
C UNKNOWNS OF THE TENSOR MODEL TSQFCN
C TSQFCN
C SUBPROGRAMS CALLED: TSQFCN
C TSQFCN
C LEVEL 1 BLAS ... DNRM2 TSQFCN
C LEVEL 2 BLAS ... DGEMV TSQFCN
C TENSOLVE ... TSSTMX TSQFCN
C TSQFCN
C********************************************************************* TSQFCN
TSQFCN
INTEGER I TSQFCN
DOUBLE PRECISION SUM,ZERO,HALF,ONE TSQFCN
DOUBLE PRECISION DNRM2 TSQFCN
DATA ZERO,HALF,ONE/0.0D0,0.50D0,1.0D0/ TSQFCN
TSQFCN
c compute the lower right (m-n+q) x p submatrix of AJA times X TSQFCN
TSQFCN
CALL DGEMV('N',M-N+Q,P,ONE,AJA(N-Q+1,N-P+1),MAXM, TSQFCN
+ X,1,ZERO,WRK1,1) TSQFCN
TSQFCN
c compute S-trans times X TSQFCN
TSQFCN
CALL TSSTMX(S,X,MAXN,N,P,WRK2,WRK3) TSQFCN
TSQFCN
c compute 0.5 * (S-trans times X)**2 TSQFCN
TSQFCN
DO 10 I = 1, P TSQFCN
WRK2(I) = HALF * WRK3(I)**2 TSQFCN
10 CONTINUE TSQFCN
TSQFCN
c compute 0.5 * (down (m-n+q) x p submatrix of ANLS) * TSQFCN
c (S-trans times X)**2 TSQFCN
TSQFCN
CALL DGEMV('N',M-N+Q,P,ONE,ANLS(N-Q+1,1),MAXM, TSQFCN
+ WRK2,1,ZERO,WRK4,1) TSQFCN
TSQFCN
DO 20 I = 1,M-N+Q TSQFCN
WRK5(I) = WRK4(I)+F(N-Q+I)+WRK1(I) TSQFCN
20 CONTINUE TSQFCN
TSQFCN
SUM = HALF*DNRM2(M-N+Q,WRK5,1)**2 TSQFCN
TSQFCN
RETURN TSQFCN
END TSQFCN
SUBROUTINE TSQLFC(QL,NR,M,N,IERR) TSQLFC
TSQLFC
INTEGER NR,M,N,IERR TSQLFC
DOUBLE PRECISION QL(NR,N) TSQLFC
TSQLFC
C********************************************************************** TSQLFC
C THIS ROUTINE PERFORMS A QL DECOMPOSITION. TSQLFC
C********************************************************************** TSQLFC
C TSQLFC
C INPUT PARAMETERS : TSQLFC
C ---------------- TSQLFC
C TSQLFC
C NR : LEADING DIMENSION OF MATRIX QL TSQLFC
C M : ROW DIMENSION OF QL TSQLFC
C N : COLUMN DIMENSION OF QL TSQLFC
C TSQLFC
C INPUT-OUTPUT PARAMETERS : TSQLFC
C ----------------------- TSQLFC
C TSQLFC
C QL : INPUT MATRIX ON ENTRY AND FACTORED MATRIX ON EXIT TSQLFC
C TSQLFC
C OUTPUT PARAMETERS : TSQLFC
C ------------------ TSQLFC
C TSQLFC
C IERR : RETURN CODE WITH THE FOLLOWING MEANINGS : TSQLFC
C IERR = 1 : SINGULARITY DETECTED TSQLFC
C IERR = 0 : OTHERWISE TSQLFC
C TSQLFC
C SUBPROGRAMS CALLED: TSQLFC
C TSQLFC
C LEVEL 1 BLAS ... DAXPY,DDOT,DSCAL TSQLFC
C TSQLFC
C********************************************************************** TSQLFC
TSQLFC
INTEGER I,J,K TSQLFC
DOUBLE PRECISION NU,SIGMA,TAU,RNU,ZERO,ONE TSQLFC
INTRINSIC ABS,MAX,SQRT TSQLFC
DOUBLE PRECISION DDOT,DNRM2 TSQLFC
DATA ZERO,ONE/0.0D0,1.0D0/ TSQLFC
TSQLFC
c zero out rows m+1 and m+2 of matrix QL TSQLFC
TSQLFC
DO 10 J = 1,N TSQLFC
QL(M+1,J) = ZERO TSQLFC
QL(M+2,J) = ZERO TSQLFC
10 CONTINUE TSQLFC
TSQLFC
IERR = 0 TSQLFC
TSQLFC
K = 1 TSQLFC
TSQLFC
20 CONTINUE TSQLFC
IF((K.LT.M).AND.(K.LE.N)) THEN TSQLFC
TSQLFC
c find NU = max element of col K on or above l-diagonal TSQLFC
TSQLFC
NU = ZERO TSQLFC
DO 40 I = 1,M+1-K TSQLFC
NU = MAX(NU,ABS(QL(I,K))) TSQLFC
40 CONTINUE TSQLFC
TSQLFC
IF(NU.NE.ZERO) THEN TSQLFC
TSQLFC
c normalize col K on or above l-diagonal TSQLFC
TSQLFC
RNU = ONE/NU TSQLFC
CALL DSCAL(M-K+1,RNU,QL(1,K),1) TSQLFC
TSQLFC
c code to find SIGMA = SGN(QL(M+1-K,K))*l2-norm of col K on or TSQLFC
c above l-diagonal TSQLFC
TSQLFC
SIGMA = DNRM2(M-K+1,QL(1,K),1) TSQLFC
SIGMA = SIGN(SIGMA,QL(M+1-K,K)) TSQLFC
TSQLFC
c store last element(1st in normal QR) of U-vector in QL(M+1-K,K) TSQLFC
TSQLFC
QL(M+1-K,K) = QL(M+1-K,K)+SIGMA TSQLFC
TSQLFC
c store value of <U,U>/2 in QL(M+1,K) TSQLFC
TSQLFC
QL(M+1,K) = SIGMA*QL(M+1-K,K) TSQLFC
IF(QL(M+1,K).EQ.ZERO) THEN TSQLFC
IERR = 1 TSQLFC
RETURN TSQLFC
ENDIF TSQLFC
TSQLFC
c store L(M+1-K,K) in QL(M+2,K) TSQLFC
TSQLFC
QL(M+2,K) = -NU*SIGMA TSQLFC
TSQLFC
c code to get (I-2U*UT/<U,U>)*QL for kth iteration TSQLFC
TSQLFC
IF(K.LT.N) THEN TSQLFC
DO 50 J = K+1,N TSQLFC
TSQLFC
c loop to get TAU = <U,J-TH COL OF QL> TSQLFC
TSQLFC
TAU = DDOT(M-K+1,QL(1,K),1,QL(1,J),1) TSQLFC
TAU = -TAU/QL(M+1,K) TSQLFC
TSQLFC
c loop to get (I-2U*UT/<U,U>)*j-th col of QL TSQLFC
TSQLFC
CALL DAXPY(M-K+1,TAU,QL(1,K),1,QL(1,J),1) TSQLFC
TSQLFC
50 CONTINUE TSQLFC
ENDIF TSQLFC
K = K+1 TSQLFC
ELSE TSQLFC
IERR = 1 TSQLFC
RETURN TSQLFC
ENDIF TSQLFC
TSQLFC
GOTO 20 TSQLFC
TSQLFC
ENDIF TSQLFC
TSQLFC
IF(M.EQ.N) QL(M+2,N) = QL(1,N) TSQLFC
TSQLFC
IF(QL(M+2,N).EQ.ZERO) THEN TSQLFC
IERR = 1 TSQLFC
ENDIF TSQLFC
TSQLFC
RETURN TSQLFC
END TSQLFC
SUBROUTINE TSQMLV(NR,N,P,Q,V,W,TRANS) TSQMLV
TSQMLV
INTEGER NR,N,P TSQMLV
DOUBLE PRECISION Q(NR,P),V(N),W(N) TSQMLV
TSQMLV
C********************************************************************** TSQMLV
C THIS ROUTINE MULTIPLYS AN ORTHOGONAL MATRTIX Q OR ITS TRANSPOSE BY TSQMLV
C A VECTOR. TSQMLV
C********************************************************************** TSQMLV
C TSQMLV
C INPUT PARAMETERS : TSQMLV
C ---------------- TSQMLV
C TSQMLV
C NR : LEADING DIMENSION OF MATRIX Q TSQMLV
C N : DIMENSION OF VECTORS V AND W TSQMLV
C P : COLUMN DIMENSION OF MATRIX Q TSQMLV
C Q : ORTHOGONAL MATRIX (OBTAINED FROM TSQLFC SUBROUTINE) TSQMLV
C V : VECTOR TO BE MULTIPLIED BY Q TSQMLV
C TRANS : BOOLEAN PARAMETER: TSQMLV
C = TRUE : PERFORM Q-TRANS*V TSQMLV
C = FALSE : PERFORM Q*V TSQMLV
C TSQMLV
C OUTPUT PARAMETERS : TSQMLV
C ----------------- TSQMLV
C TSQMLV
C W : VECTOR Q*V OR Q-TRANS*V ON EXIT TSQMLV
C TSQMLV
C SUBPROGRAMS CALLED: TSQMLV
C TSQMLV
C LEVEL 1 BLAS ... DAXPY,DCOPY,DDOT TSQMLV
C TSQMLV
C********************************************************************** TSQMLV
TSQMLV
INTEGER J,K TSQMLV
DOUBLE PRECISION TAU,CONST TSQMLV
LOGICAL TRANS TSQMLV
DOUBLE PRECISION DDOT TSQMLV
TSQMLV
CALL DCOPY(N,V,1,W,1) TSQMLV
TSQMLV
DO 40 J = 1,P TSQMLV
IF(TRANS) THEN TSQMLV
K = P+1-J TSQMLV
ELSE TSQMLV
K = J TSQMLV
ENDIF TSQMLV
TAU = DDOT(N-K+1,Q(1,K),1,W,1) TSQMLV
CONST = -TAU/Q(N+1,K) TSQMLV
CALL DAXPY(N-K+1,CONST,Q(1,K),1,W,1) TSQMLV
40 CONTINUE TSQMLV
TSQMLV
RETURN TSQMLV
END TSQMLV
SUBROUTINE TSQMTS(ANLS,QHAT,NR,MJ,N,M,P,LB,WRK1,ZERO1) TSQMTS
TSQMTS
INTEGER NR,M,N,P,MJ,LB,ZERO1 TSQMTS
DOUBLE PRECISION ANLS(NR,P),QHAT(NR,N),WRK1(M) TSQMTS
TSQMTS
C********************************************************************** TSQMTS
C THIS ROUTINE MULTIPLIES AN ORTHOGONAL MATRIX QHAT BY THE TENSOR TSQMTS
C MATRIX ANLS. TSQMTS
C********************************************************************** TSQMTS
C TSQMTS
C INPUT PARAMETERS : TSQMTS
C ---------------- TSQMTS
C TSQMTS
C QHAT : ORTHOGONAL MATRIX (OBTAINED FROM TSQRFC SUBROUTINE) TSQMTS
C NR : LEADIND DIMENSION OF MATRIX QHAT TSQMTS
C MJ : ROW DIMENSION OF QHAT TSQMTS
C N : COLUMN DIMENSION OF MATRIX QHAT TSQMTS
C M : ROW DIMENSION OF MATRIX ANLS TSQMTS
C P : COLUMN DIMENSION OF MATRIX ANLS TSQMTS
C LB : STARTING COLUMN FROM WHICH QR DECOMPOSITION WAS PERFORMED TSQMTS
C WRK1 : WORKING VECTOR TSQMTS
C ZERO1: FIRST ZERO COLUMN OF MATRIX QHAT IN CASE OF SINGULARITY TSQMTS
C TSQMTS
C INPUT-OUTPUT PARAMETERS : TSQMTS
C ------------------------ TSQMTS
C TSQMTS
C ANLS : MATRIX TO BE MULTIPLIED BY AN ORTHOGONAL MATRIX TSQMTS
C ON ENTRY AND THE MATRIX QHAT*ANLS ON EXIT TSQMTS
C TSQMTS
C SUBPROGRAMS CALLED: TSQMTS
C TSQMTS
C LEVEL 1 BLAS ... DCOPY TSQMTS
C TENSOLVE ... TSQMUV TSQMTS
C TSQMTS
C********************************************************************** TSQMTS
TSQMTS
INTEGER J TSQMTS
TSQMTS
DO 40 J = 1,P TSQMTS
CALL TSQMUV(QHAT,ANLS(1,J),WRK1,NR,MJ,LB,ZERO1,.FALSE.) TSQMTS
CALL DCOPY(M,WRK1,1,ANLS(1,J),1) TSQMTS
40 CONTINUE TSQMTS
TSQMTS
RETURN TSQMTS
END TSQMTS
SUBROUTINE TSQMUV(Q,V,W,NR,M,LB,ZERO1,TRANSP) TSQMUV
TSQMUV
INTEGER NR,M,LB,ZERO1 TSQMUV
DOUBLE PRECISION Q(NR,*),V(M),W(M) TSQMUV
LOGICAL TRANSP TSQMUV
TSQMUV
C********************************************************************** TSQMUV
C THIS SUBROUTINE MULTIPLIES AN ORTHOGONAL MATRIX BY A VECTOR. TSQMUV
C********************************************************************** TSQMUV
C TSQMUV
C INPUT PARAMETERS : TSQMUV
C ----------------- TSQMUV
C TSQMUV
C Q : FACTORED MATRIX (OBTAINED FROM SUBROUTINE TSQRFC) TSQMUV
C V : VECTOR TO BE MULTIPLIED BY THE ORTHOGONAL MATRIX Q TSQMUV
C NR : LEADING DIMENSION OF MATRIX Q TSQMUV
C M : ROW DIMENSION OF MATRIX Q TSQMUV
C LB : STARTING COLUMN FROM WHICH QR DECOMPOSITION WAS PERFORMED TSQMUV
C ZERO1: FIRST ZERO COLUMN OF THE MATRIX Q TSQMUV
C TRANSP : BOOLEAN PARAMETER : TSQMUV
C = TRUE : PERFORM Q-TRANS*V TSQMUV
C = FALSE : PERFORM Q*V TSQMUV
C TSQMUV
C OUTPUT PARAMETERS : TSQMUV
C ----------------- TSQMUV
C TSQMUV
C W : Q*V OR Q-TRANS*V ON EXIT TSQMUV
C TSQMUV
C SUBPROGRAMS CALLED: TSQMUV
C TSQMUV
C LEVEL 1 BLAS ... DAXPY,DCOPY,DDOT TSQMUV
C TSQMUV
C******************************************************************** TSQMUV
C TSQMUV
INTEGER UB,A,B,C,K TSQMUV
DOUBLE PRECISION TAU,CONST TSQMUV
DOUBLE PRECISION DDOT TSQMUV
TSQMUV
c copy the vector V to W TSQMUV
TSQMUV
CALL DCOPY(M,V,1,W,1) TSQMUV
TSQMUV
UB = ZERO1-1 TSQMUV
IF(TRANSP) THEN TSQMUV
A = UB TSQMUV
B = LB TSQMUV
C = -1 TSQMUV
ELSE TSQMUV
A = LB TSQMUV
B = UB TSQMUV
C = 1 TSQMUV
ENDIF TSQMUV
TSQMUV
DO 50 K = A,B,C TSQMUV
TAU = DDOT(M-K+1,Q(K,K),1,W(K),1) TSQMUV
CONST = -TAU/Q(M+1,K) TSQMUV
CALL DAXPY(M-K+1,CONST,Q(K,K),1,W(K),1) TSQMUV
50 CONTINUE TSQMUV
TSQMUV
RETURN TSQMUV
END TSQMUV
SUBROUTINE TSQRFC(QR,NR,N,M,LB,UB,IERR,EPSM,AL2NRM,CURPOS,ZERO1) TSQRFC
TSQRFC
INTEGER NR,N,M,LB,UB,IERR,ZERO1 TSQRFC
INTEGER CURPOS(N) TSQRFC
DOUBLE PRECISION QR(NR,N),AL2NRM(M),EPSM TSQRFC
TSQRFC
C********************************************************************** TSQRFC
C THIS ROUTINE PERFORMS COLUMN-PIVOTED QR DECOMPOSITION ON AN M*N TSQRFC
C MATRIX. THE DECOMPOSITION IS DONE BETWEEN COLS LB AND UB. TSQRFC
C********************************************************************** TSQRFC
C TSQRFC
C INPUT PARAMETERS : TSQRFC
C ----------------- TSQRFC
C TSQRFC
C NR : LEADING DIMENSION OF MATRIX QR TSQRFC
C N : COLUMN DIMENSION OF MATRIX QR TSQRFC
C M : ROW DIMENSION OF MATRIX QR TSQRFC
C LB,UB : SUBSPACE OF QR DECOMPOSITION TSQRFC
C EPSM : MACHINE PRECISION TSQRFC
C AL2NRM: WORKING VECTOR TSQRFC
C TSQRFC
C INPUT-OUTPUT PARAMETERS : TSQRFC
C ------------------------ TSQRFC
C QR : INPUT MATRIX ON ENTRY AND FACTORED MATRIX ON EXIT TSQRFC
C TSQRFC
C OUTPUT PARAMETERS : TSQRFC
C ------------------ TSQRFC
C TSQRFC
C IERR : RETURN CODE WITH TH FOLLOWING MEANINGS: TSQRFC
C IERR = 1 : SINGULARITY DETECTED TSQRFC
C IERR = 0 : OTHERWISE TSQRFC
C CURPOS : PIVOT VECTOR TSQRFC
C ZERO1 : FIRST ZERO COLUMN OF MATRIX QR IN CASE OF SINGULARITY TSQRFC
C TSQRFC
C SUBPROGRAMS CALLED: TSQRFC
C TSQRFC
C LEVEL 1 BLAS ... DAXPY,DDOT,DNRM2,DSCAL,DSWAP,IDAMAX TSQRFC
C TSQRFC
C **********************************************************************TSQRFC
TSQRFC
INTEGER COLPIV,I,J,K,L TSQRFC
DOUBLE PRECISION COLMAX,SIGMA,TAU,AMAX TSQRFC
DOUBLE PRECISION NU,RNU,ZERO,ONE TSQRFC
DOUBLE PRECISION DDOT,DNRM2 TSQRFC
INTEGER IDAMAX TSQRFC
INTRINSIC ABS,MAX,SIGN,SQRT TSQRFC
DATA ZERO,ONE/0.0D0,1.0D0/ TSQRFC
TSQRFC
c zero rows m+1 and m+2 of QR matrix TSQRFC
TSQRFC
DO 10 J = 1,N TSQRFC
CURPOS(J) = J TSQRFC
10 CONTINUE TSQRFC
TSQRFC
DO 20 J = LB,UB TSQRFC
QR(M+1,J) = ZERO TSQRFC
QR(M+2,J) = ZERO TSQRFC
20 CONTINUE TSQRFC
TSQRFC
IERR = 0 TSQRFC
ZERO1 = UB+1 TSQRFC
K = LB TSQRFC
TSQRFC
c get L2NORM**2 of columns (LB to UB) TSQRFC
TSQRFC
DO 30 J = K,UB TSQRFC
AL2NRM(J) = DNRM2(M-K+1,QR(K,J),1)**2 TSQRFC
30 CONTINUE TSQRFC
TSQRFC
40 CONTINUE TSQRFC
IF((K.LT.M).AND.(K.LE.UB)) THEN TSQRFC
TSQRFC
AMAX = ZERO TSQRFC
DO 60 J = K,UB TSQRFC
IF(AL2NRM(J).GE.AMAX) THEN TSQRFC
AMAX = AL2NRM(J) TSQRFC
COLPIV = J TSQRFC
ENDIF TSQRFC
60 CONTINUE TSQRFC
TSQRFC
IF(AMAX.EQ.ZERO) THEN TSQRFC
IERR = 1 TSQRFC
ZERO1 = K TSQRFC
RETURN TSQRFC
ENDIF TSQRFC
TSQRFC
IF(K.EQ.LB) THEN TSQRFC
COLMAX = AMAX TSQRFC
ENDIF TSQRFC
TSQRFC
IF(AL2NRM(COLPIV).LE.EPSM*COLMAX) THEN TSQRFC
IERR = 1 TSQRFC
ZERO1 = K TSQRFC
RETURN TSQRFC
ENDIF TSQRFC
TSQRFC
IF(COLPIV .NE. K) THEN TSQRFC
CALL DSWAP(M+2,QR(1,COLPIV),1,QR(1,K),1) TSQRFC
L = CURPOS(K) TSQRFC
CURPOS(K) = CURPOS(COLPIV) TSQRFC
CURPOS(COLPIV) = L TSQRFC
CALL DSWAP(1,AL2NRM(COLPIV),1,AL2NRM(K),1) TSQRFC
ENDIF TSQRFC
TSQRFC
c find NU = max element of col K on or below diagonal TSQRFC
TSQRFC
L = IDAMAX(M-K+1,QR(K,K),1) + K - 1 TSQRFC
NU = ABS(QR(L,K)) TSQRFC
TSQRFC
IF(NU.EQ.ZERO) THEN TSQRFC
IERR = 1 TSQRFC
ZERO1 = K TSQRFC
RETURN TSQRFC
ENDIF TSQRFC
TSQRFC
c normalize col K on or below diagonal TSQRFC
TSQRFC
RNU = ONE/NU TSQRFC
CALL DSCAL(M-K+1,RNU,QR(K,K),1) TSQRFC
TSQRFC
c code to find SIGMA = SGN(QR(K,K))*l2-norm of col K on or TSQRFC
c below diagonal TSQRFC
TSQRFC
SIGMA = DNRM2(M-K+1,QR(K,K),1) TSQRFC
SIGMA = SIGN(SIGMA,QR(K,K)) TSQRFC
TSQRFC
c store 1st element of U-vector in QR(K,K) TSQRFC
TSQRFC
QR(K,K) = QR(K,K)+SIGMA TSQRFC
TSQRFC
c store value of <U,U>/2 in QR(M+1,K) TSQRFC
TSQRFC
QR(M+1,K) = SIGMA*QR(K,K) TSQRFC
IF(QR(M+1,K).EQ.ZERO) THEN TSQRFC
IERR = 1 TSQRFC
ZERO1 = K TSQRFC
RETURN TSQRFC
ENDIF TSQRFC
TSQRFC
c store R(K,K) in QR(M+2,K) TSQRFC
TSQRFC
QR(M+2,K) = -NU*SIGMA TSQRFC
IF(QR(M+2,K).EQ.ZERO) THEN TSQRFC
IERR = 1 TSQRFC
ZERO1 = K TSQRFC
RETURN TSQRFC
ENDIF TSQRFC
TSQRFC
c code to get (I-2U*UT/<U,U>)*QR for kth iteration TSQRFC
TSQRFC
IF(K.LT.N) THEN TSQRFC
DO 130 J = K+1,N TSQRFC
TSQRFC
c loop to get UT*J-TH col of QR TSQRFC
TSQRFC
TAU = DDOT(M-K+1,QR(K,K),1,QR(K,J),1) TSQRFC
TAU = -TAU/QR(M+1,K) TSQRFC
TSQRFC
c loop to get (I-2U*UT/<U,U>)*j-th col of QR TSQRFC
TSQRFC
CALL DAXPY(M-K+1,TAU,QR(K,K),1,QR(K,J),1) TSQRFC
TSQRFC
130 CONTINUE TSQRFC
ENDIF TSQRFC
TSQRFC
c update l2norm**2 (K+1 to UB) TSQRFC
TSQRFC
DO 140 I = K+1,UB TSQRFC
AL2NRM(I) = AL2NRM(I)-QR(K,I)**2 TSQRFC
140 CONTINUE TSQRFC
TSQRFC
K = K+1 TSQRFC
GOTO 40 TSQRFC
TSQRFC
ELSE TSQRFC
TSQRFC
IF(LB.EQ.UB) COLMAX = AL2NRM(LB) TSQRFC
TSQRFC
ENDIF TSQRFC
TSQRFC
IF(M.EQ.UB) QR(M+2,UB) = QR(M,M) TSQRFC
IF(ABS(QR(M+2,UB)).LE.EPSM*COLMAX) THEN TSQRFC
IERR = 1 TSQRFC
ZERO1 = UB TSQRFC
ENDIF TSQRFC
TSQRFC
RETURN TSQRFC
END TSQRFC
SUBROUTINE TSRSID(ITN,METHOD,FQ,D,CURPOS,PIVOT,PBAR,AJA,ANLS, TSRSID
+ SHAT,FLAG,NWTAKE,IERR,MAXM,MAXN,M,N,P,WRK1, TSRSID
+ VN,VNP,VNS,SCRES) TSRSID
TSRSID
INTEGER MAXM,MAXN,M,N,P,IERR,FLAG,ITN,METHOD TSRSID
INTEGER CURPOS(N),PIVOT(N),PBAR(N) TSRSID
DOUBLE PRECISION SCRES,D(N),VN(M),VNP(M),VNS(M),AJA(MAXM,N) TSRSID
DOUBLE PRECISION WRK1(M),SHAT(MAXN,P),FQ(M) TSRSID
DOUBLE PRECISION ANLS(MAXM,P) TSRSID
LOGICAL NWTAKE TSRSID
TSRSID
C********************************************************************** TSRSID
C THIS ROUTINE COMPUTES || F + J*D + 0.5*A*D**2 ||**2 IN THE L2 TSRSID
C NORM SENS AT A GIVEN STEP D. TSRSID
C********************************************************************** TSRSID
C TSRSID
C INPUT PARAMETERS TSRSID
C ---------------- TSRSID
C TSRSID
C ITN : CURRENT ITERATION NUMBER TSRSID
C METHOD: METHOD TO BE USED TSRSID
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSRSID
C ORTHOGONAL MATRICES TSRSID
C D : STEP AT WHICH TO EVALUATE THE TENSOR MODEL TSRSID
C CURPOS: PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA TSRSID
C FROM COLUMN 1 TO N-P) TSRSID
C PIVOT : PIVOT VECTOR ( USED DURING THE FACTORIZATION OF AJA TSRSID
C FROM COLUMN N-P+1 TO N) TSRSID
C PBAR : PIVOT VECTOR (USED DURING THE FACTORIZATION OF AJA TSRSID
C IF IT IS SINGULAR TSRSID
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSRSID
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSRSID
C SHAT : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS AFTER TSRSID
C A QL FACTORIZATION TSRSID
C FLAG : RETURN CODE WITH THE FOLLOWING MEANINGS: TSRSID
C FLAG = 0 : NO SINGULARITY DETECTED DURING FACTORIZATIONTSRSID
C OF THE JACOBIAN FROM COLUMN 1 TO N TSRSID
C FLAG = 1 : SINGULARITY DETECTED DURING FACTORIZATION TSRSID
C OF THE JACOBIAN FROM COLUMN 1 TO N-P TSRSID
C FLAG = 2 : SINGULARITY DETECTED DURING FACTORIZATION TSRSID
C OF THE JACOBIAN FROM COLUMN N-P+1 TO N TSRSID
C NWTAKE: LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS: TSRSID
C NWTAKE = .TRUE. : NEWTON STEP TAKEN TSRSID
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSRSID
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE: TSRSID
C IERR = 0 : NO SINGULARITY DETECTED TSRSID
C IERR = 1 : SINGULARITY DETECTED TSRSID
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSRSID
C MAXN : LEADING DIMENSION OF SHAT TSRSID
C M,N : DIMENSIONS OF THE PROBLEM TSRSID
C P : COLUMN DIMENSION OF THE MATRICES SHAT AND ANLS TSRSID
C WRK1,VN,VNP,VNS : WORKSPACE VECTORS TSRSID
C TSRSID
C OUTPUT PARAMETERS TSRSID
C ----------------- TSRSID
C TSRSID
C SCRES : || F + J*D + 0.5*A*D**2 ||**2 TSRSID
C TSRSID
C SUBPROGRAMS CALLED: TSRSID
C TSRSID
C LEVEL 1 BLAS ... DLOAD,DNRM2 TSRSID
C LEVEL 2 BLAS ... DGEMV TSRSID
C TENSOLVE ... TSJMUV,TSUDQV TSRSID
C TSRSID
C **********************************************************************TSRSID
TSRSID
INTEGER I,J,LEN TSRSID
DOUBLE PRECISION ZERO,HALF,ONE TSRSID
DOUBLE PRECISION DNRM2 TSRSID
DATA ZERO,HALF,ONE/0.0D0,0.50D0,1.0D0/ TSRSID
TSRSID
CALL TSJMUV(ITN,METHOD,D,CURPOS,PIVOT,PBAR,AJA,SHAT,FLAG, TSRSID
+ IERR,MAXM,MAXN,M,N,P,VN,VNP,VNS,WRK1) TSRSID
TSRSID
CALL DLOAD (M,ZERO,WRK1(N+1),1) TSRSID
TSRSID
LEN = M TSRSID
IF(IERR .GT. 0) LEN = M + N TSRSID
TSRSID
DO 20 I = 1, LEN TSRSID
VN(I) = WRK1(I) + FQ(I) TSRSID
20 CONTINUE TSRSID
TSRSID
IF( .NOT. NWTAKE) THEN TSRSID
CALL TSUDQV(SHAT,VNS,MAXN,N,P,VNP) TSRSID
DO 30 J = 1, P TSRSID
VNP(J) = VNP(J) ** 2 TSRSID
30 CONTINUE TSRSID
CALL DGEMV('N',LEN,P,HALF,ANLS,MAXM,VNP,1,ONE,VN,1) TSRSID
ENDIF TSRSID
TSRSID
SCRES = DNRM2(LEN,VN,1) TSRSID
TSRSID
RETURN TSRSID
END TSRSID
SUBROUTINE TSSCLF(F,DF,M) TSSCLF
TSSCLF
INTEGER M TSSCLF
DOUBLE PRECISION F(M),DF(M) TSSCLF
TSSCLF
C******************************************************************* TSSCLF
C THIS ROUTINE SCALES A FUNCTION VALUE F. TSSCLF
C******************************************************************* TSSCLF
C TSSCLF
C INPUT PARAMETERS : TSSCLF
C ------------------ TSSCLF
C TSSCLF
C DF : DIAGONAL SCALING MATRIX FOR F TSSCLF
C M : DIMENSION OF F TSSCLF
C TSSCLF
C INPUT-OUTPUT PARAMETERS : TSSCLF
C ------------------ TSSCLF
C TSSCLF
C F : UNSCALED FUNCTION VALUE ON ENTRY AND SCALED FUNCTION TSSCLF
C VALUE ON EXIT TSSCLF
C TSSCLF
C********************************************************************* TSSCLF
TSSCLF
INTEGER I TSSCLF
TSSCLF
DO 10 I = 1,M TSSCLF
F(I) = DF(I)*F(I) TSSCLF
10 CONTINUE TSSCLF
TSSCLF
RETURN TSSCLF
END TSSCLF
SUBROUTINE TSSCLJ(X,DX,TYPX,F,DF,FHAT,NR,M,N,EPSM,JACFLG, TSSCLJ
+ TSFVEC,TSANJA,AJA) TSSCLJ
TSSCLJ
INTEGER NR,M,N,JACFLG TSSCLJ
DOUBLE PRECISION EPSM TSSCLJ
DOUBLE PRECISION X(N),DX(N),TYPX(N),F(M),DF(M) TSSCLJ
DOUBLE PRECISION AJA(NR,N),FHAT(M) TSSCLJ
EXTERNAL TSFVEC,TSANJA TSSCLJ
TSSCLJ
C********************************************************************** TSSCLJ
C THIS ROUTINE COMPUTES THE JACOBIAN MATRIX AT THE CURRENT ITERATE TSSCLJ
C X AND SCALES ITS VALUE. TSSCLJ
C********************************************************************** TSSCLJ
C TSSCLJ
C INPUT PARAMETERS : TSSCLJ
C ----------------- TSSCLJ
C TSSCLJ
C X : SCALED CURRENT ITERATE TSSCLJ
C DX : DIAGONAL SCALING MATRIX FOR X TSSCLJ
C F : SCALED FUNCTION VALUE AT X TSSCLJ
C DF : DIAGONAL SCALING MATRIX FOR F TSSCLJ
C FHAT : WORKSPACE ARRAY TSSCLJ
C NR : LEADING DIMENSION OF AJA TSSCLJ
C M,N : DIMENSIONS OF PROBLEM TSSCLJ
C EPSM : MACHINE PRECISION TSSCLJ
C JACFLG : JACOBIAN FLAG TSSCLJ
C TSFVEC : SUBROUTINE TO EVALUATE FUNCTION TSSCLJ
C TSANJA : SUBROUTINE TO EVALUATE ANALYTIC JACOBIAN TSSCLJ
C TSSCLJ
C TSSCLJ
C INPUT-OUTPUT PARAMETERS : TSSCLJ
C ------------------------ TSSCLJ
C TSSCLJ
C AJA : SCALED JACOBIAN AT CURRENT ITERATE TSSCLJ
C TSSCLJ
C SUBPROGRAMS CALLED: TSSCLJ
C TSSCLJ
C TENSOLVE ... TSUNSX,TSUNSF,TSFDFJ,TSSCLF,TSSCLX TSSCLJ
C USER ... TSANJA TSSCLJ
C TSSCLJ
C******************************************************************** TSSCLJ
TSSCLJ
INTEGER I,J TSSCLJ
TSSCLJ
c unscale X AND F TSSCLJ
TSSCLJ
CALL TSUNSX(X,DX,N) TSSCLJ
CALL TSUNSF(F,DF,M) TSSCLJ
TSSCLJ
c compute the finite difference or analytic Jacobian at X TSSCLJ
TSSCLJ
IF(JACFLG.EQ.0) THEN TSSCLJ
CALL TSFDFJ(X,F,NR,M,N,EPSM,TSFVEC,FHAT,AJA) TSSCLJ
ELSE TSSCLJ
CALL TSANJA(AJA,X,NR,M,N) TSSCLJ
ENDIF TSSCLJ
TSSCLJ
c scale the Jacobian matrix TSSCLJ
TSSCLJ
DO 20 J = 1,N TSSCLJ
DO 10 I = 1,M TSSCLJ
AJA(I,J) = AJA(I,J)*DF(I)*TYPX(J) TSSCLJ
10 CONTINUE TSSCLJ
20 CONTINUE TSSCLJ
TSSCLJ
c scale back X AND F TSSCLJ
TSSCLJ
CALL TSSCLF(F,DF,M) TSSCLJ
CALL TSSCLX(X,DX,N) TSSCLJ
TSSCLJ
RETURN TSSCLJ
END TSSCLJ
SUBROUTINE TSSCLX(X,DX,N) TSSCLX
TSSCLX
INTEGER N TSSCLX
DOUBLE PRECISION X(N),DX(N) TSSCLX
TSSCLX
C********************************************************************** TSSCLX
C THIS ROUTINE SCALES A VECTOR X. TSSCLX
C********************************************************************** TSSCLX
C TSSCLX
C INPUT PARAMETERS : TSSCLX
C ------------------ TSSCLX
C TSSCLX
C DX : DIAGONAL SCALING MATRIX FOR X TSSCLX
C N : DIMENSION OF X TSSCLX
C TSSCLX
C OUTPUT PARAMETERS : TSSCLX
C ------------------ TSSCLX
C TSSCLX
C X : SCALED VECTOR X TSSCLX
C TSSCLX
C********************************************************************** TSSCLX
TSSCLX
INTEGER I TSSCLX
TSSCLX
DO 10 I = 1,N TSSCLX
X(I) = DX(I)*X(I) TSSCLX
10 CONTINUE TSSCLX
TSSCLX
RETURN TSSCLX
END TSSCLX
SUBROUTINE TSSLCT(RESIDT,RESIDN,ITRMCD,FC,M,N,DN,DT,G,DF,NWTAKE) TSSLCT
TSSLCT
INTEGER M,N,ITRMCD TSSLCT
DOUBLE PRECISION RESIDT,RESIDN TSSLCT
DOUBLE PRECISION FC(M),DF(N),DN(N),DT(N),G(N) TSSLCT
LOGICAL NWTAKE TSSLCT
TSSLCT
C********************************************************************* TSSLCT
C THIS ROUTINE DECIDES WHICH DIRECTION TO CHOOSE: THE TENSOR OR THE TSSLCT
C STANDARD DIRECTION. THE STANDARD DIRECTION IS CHOSEN WHENEVER TSSLCT
C THE TENSOR DIRECTION IS NOT DESCENT OR THE TENSOR DIRECTION IS TO A TSSLCT
C MINIMIZER OF THE TENSOR MODEL AND DOESN'T PROVIDE ENOUGH DECREASE TSSLCT
C IN THE TENSOR MODEL, OR THE QUADRATIC SYSTEM OF Q EQUATIONS IN P TSSLCT
C UNKNOWNS CANNOT BE SOLVED BY UNCMIN WITHIN 150 ITERATIONS. TSSLCT
C********************************************************************* TSSLCT
C TSSLCT
C INPUT PARAMETERS : TSSLCT
C ------------------ TSSLCT
C TSSLCT
C RESIDT : TENSOR RESIDUAL TSSLCT
C RESIDN : NEWTON RESIDUAL TSSLCT
C ITRMCD : UNCMIN TERMINATION CODE TSSLCT
C FC : FUNCTION VALUE AT XC TSSLCT
C M,N: DIMENSIONS OF PROBLEM TSSLCT
C DN : STANDARD STEP TSSLCT
C DT : TENSOR STEP TSSLCT
C G : GRADIENT VALUE AT XC TSSLCT
C TSSLCT
C TSSLCT
C OUTPUT PARAMETERS : TSSLCT
C ----------------- TSSLCT
C TSSLCT
C DF : EITHER THE STANDARD OR TENSOR STEP ON EXIT TSSLCT
C NWTAKE : BOOLEAN VALUE WITH THE FOLLOWING MEANINGS: TSSLCT
C =.TRUE. : STANDARD STEP IS TAKEN TSSLCT
C =.FALSE. : TENSOR STEP IS TAKEN TSSLCT
C TSSLCT
C SUBPROGRAMS CALLED: TSSLCT
C TSSLCT
C LEVEL 1 BLAS .... DCOPY,DDOT,DNRM2 TSSLCT
C TSSLCT
C********************************************************************* TSSLCT
TSSLCT
DOUBLE PRECISION ANRMFC,DTNORM,GNORM TSSLCT
DOUBLE PRECISION TEMP,TEMP1,BETA,GAMA TSSLCT
DOUBLE PRECISION TENTH,ONETT,HALF TSSLCT
DOUBLE PRECISION DNRM2,DDOT TSSLCT
DATA ONETT,TENTH,HALF/1.0D-4,1.0D-1,0.50D0/ TSSLCT
TSSLCT
NWTAKE = .FALSE. TSSLCT
ANRMFC = DNRM2(M,FC,1) TSSLCT
DTNORM = DNRM2(N,DT,1) TSSLCT
GNORM = DNRM2(N,G,1) TSSLCT
TEMP = DDOT(N,DT,1,G,1) TSSLCT
TSSLCT
GAMA = HALF TSSLCT
IF(M.GT.N) THEN TSSLCT
BETA = TENTH TSSLCT
ELSE TSSLCT
BETA = ONETT TSSLCT
ENDIF TSSLCT
TSSLCT
TEMP1 = -BETA*DTNORM*GNORM TSSLCT
TSSLCT
IF(RESIDT.GE.GAMA*(ANRMFC+RESIDN).OR.(TEMP.GT.TEMP1).OR. TSSLCT
+ ITRMCD.EQ.4) THEN TSSLCT
CALL DCOPY(N,DN,1,DF,1) TSSLCT
NWTAKE = .TRUE. TSSLCT
ELSE TSSLCT
CALL DCOPY(N,DT,1,DF,1) TSSLCT
ENDIF TSSLCT
TSSLCT
RETURN TSSLCT
END TSSLCT
SUBROUTINE TSSMIN(ANLS,FQ,ADT,AG,CONST1,CONST2,DLT,NR,M,N, TSSMIN
+ P,NWTAKE,IERR,EPSM,VN,VNP,VNS,SOL) TSSMIN
TSSMIN
DOUBLE PRECISION DLT,EPSM TSSMIN
INTEGER NR,M,N,P,IERR TSSMIN
DOUBLE PRECISION ADT(N),AG(N),VN(M),VNP(M) TSSMIN
DOUBLE PRECISION VNS(M),ANLS(NR,P),FQ(M) TSSMIN
DOUBLE PRECISION CONST1(P),CONST2(P) TSSMIN
LOGICAL NWTAKE TSSMIN
TSSMIN
C***********************************************************************TSSMIN
C THIS ROUTINE MINIMIZES THE TENSOR MODEL OVER THE SUBSPACE SPANNED BY TSSMIN
C THE TENSOR STEP AND THE STEEPEST DECENT DIRECTION. IF THE NEWTON STEP TSSMIN
C WERE CHOSEN, IT WILL MINIMIZE THE NEWTON MODEL OVER THE SUBSPACE TSSMIN
C SPANNED BY THE NEWTON STEP AND THE STEEPEST DESCENT DIRECTION. TSSMIN
C***********************************************************************TSSMIN
C TSSMIN
C INPUT PARAMETERS : TSSMIN
C ----------------- TSSMIN
C TSSMIN
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSSMIN
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSSMIN
C ORTHOGONAL MATRICES TSSMIN
C ADT : JACOBIAN MATRIX TIMES DT (SEE SUBROUTINE TS2DTR) TSSMIN
C AG : JACOBIAN MATRIX TIMES GBAR (SEE SUBROUTINE TS2DTR) TSSMIN
C CONST1: SHAT-TRANS * DT (SEE SUBROUTINE TS2DTR) TSSMIN
C CONST2: SHAT-TRANS * GBAR (SEE SUBROUTINE TS2DTR) TSSMIN
C ALPHA : POINT AT WHICH DERIVATIVE IS EVALUATED TSSMIN
C DLT : CURRENT TRUST RADIUS TSSMIN
C NR : LEADING DIMENSION OF ANLS TSSMIN
C M,N: DIMENSIONS OF THE PROBLEM TSSMIN
C P: COLUMN DIMENSION OF MATRIX ANLS TSSMIN
C NWTAKE : LOGICAL VARIABLE WITH THE FOLLOWING MEANINGS : TSSMIN
C NWTAKE = .TRUE. : STANDARD STEP TAKEN TSSMIN
C NWTAKE = .FALSE. : TENSOR STEP TAKEN TSSMIN
C IERR : RETURN CODE FROM QRP FACTORIZATION ROUTINE TSSMIN
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSSMIN
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED TSSMIN
C EPSM : MACHINE PRECISION TSSMIN
C VN,VNP,VNS : WORKING VECTORS TSSMIN
C TSSMIN
C TSSMIN
C OUTPUT PARAMETERS : TSSMIN
C ----------------- TSSMIN
C TSSMIN
C SOL : GLOBAL MINIMIZER OF THE ONE VARIABLE FUNCTION IN ALPHA TSSMIN
C ||F + J*(ALPHA*DT + BETA*GBAR) + 0.5*A*(ALPHA*DT + TSSMIN
C BETA*GBAR)**2||**2 WHERE BETA = SQRT(DLT**2 - ALPHA**2) TSSMIN
C TSSMIN
C SUBPROGRAMS CALLED: TSSMIN
C TSSMIN
C TENSOLVE ... TSFAFA,TSMFDA,TSLMIN TSSMIN
C TSSMIN
C***********************************************************************TSSMIN
TSSMIN
INTEGER INT TSSMIN
DOUBLE PRECISION SOL,TOL,DL,S,SP,C,TSFAFA,A,TSMFDA TSSMIN
DOUBLE PRECISION D,S1,B,Q,BC,OPTIM,AC,GLOPT,QINCR TSSMIN
DOUBLE PRECISION ZERO,OHUND,TENTH,TWO,THREE,TEN TSSMIN
LOGICAL FIRST TSSMIN
DATA ZERO,OHUND,TENTH,TWO,THREE,TEN/0.0D0,1.0D-2,1.0D-1,2.0D0, TSSMIN
+ 3.0D0,10.0D0/ TSSMIN
TSSMIN
FIRST = .TRUE. TSSMIN
TOL = EPSM**(TWO/THREE) TSSMIN
INT = 40 TSSMIN
DL = TOL TSSMIN
IF(DLT.LT.TOL) THEN TSSMIN
DL = TOL*TENTH TSSMIN
ELSEIF(DLT.GT.TOL*TEN) THEN TSSMIN
DL = TOL*TEN TSSMIN
ENDIF TSSMIN
IF(DLT.LE.OHUND) THEN TSSMIN
INT = 10 TSSMIN
ENDIF TSSMIN
TSSMIN
c find global minimizer of FALPHA TSSMIN
TSSMIN
QINCR = TWO*(DLT-DL)/INT TSSMIN
DO 10 S = -DLT+DL,DLT*(INT-2)/INT,QINCR TSSMIN
SP = S TSSMIN
S1 = S+QINCR TSSMIN
TSSMIN
c evaluate FALPHA(SP) and the derivative of FALPHA at SP TSSMIN
TSSMIN
IF(FIRST) THEN TSSMIN
C = TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,SP,DLT, TSSMIN
+ NR,M,N,P,NWTAKE,IERR,VN) TSSMIN
A = TSMFDA(ANLS,ADT,AG,CONST1,CONST2,SP,DLT, TSSMIN
+ NR,M,N,P,NWTAKE,IERR,VN,VNP) TSSMIN
ELSE TSSMIN
C = D TSSMIN
A = B TSSMIN
ENDIF TSSMIN
TSSMIN
c evaluate FALPHA(S1) and the derivative of FALPHA at S1 TSSMIN
TSSMIN
D = TSFAFA(ANLS,FQ,ADT,AG,CONST1,CONST2,S1,DLT, TSSMIN
+ NR,M,N,P,NWTAKE,IERR,VN) TSSMIN
B = TSMFDA(ANLS,ADT,AG,CONST1,CONST2,S1,DLT, TSSMIN
+ NR,M,N,P,NWTAKE,IERR,VN,VNP) TSSMIN
TSSMIN
c minimize FALPHA in the subinterval [SP,S1] TSSMIN
TSSMIN
IF((A.LE.ZERO).AND.(B.GE.ZERO)) THEN TSSMIN
IF(C.GT.D) THEN TSSMIN
Q = D TSSMIN
BC = B TSSMIN
CALL TSLMIN(S1,SP,BC,Q,ANLS,FQ,ADT,AG,CONST1,CONST2, TSSMIN
+ DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP, TSSMIN
+ VNS,OPTIM) TSSMIN
ELSE TSSMIN
Q = C TSSMIN
AC = A TSSMIN
CALL TSLMIN(SP,S1,AC,Q,ANLS,FQ,ADT,AG,CONST1,CONST2, TSSMIN
+ DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP, TSSMIN
+ VNS,OPTIM) TSSMIN
ENDIF TSSMIN
ELSEIF((A.LE.ZERO).AND.(B.LE.ZERO)) THEN TSSMIN
IF(C.LE.D) THEN TSSMIN
Q = C TSSMIN
AC = A TSSMIN
CALL TSLMIN(SP,S1,AC,Q,ANLS,FQ,ADT,AG,CONST1,CONST2, TSSMIN
+ DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP, TSSMIN
+ VNS,OPTIM) TSSMIN
ELSE TSSMIN
OPTIM = S1 TSSMIN
Q = D TSSMIN
ENDIF TSSMIN
ELSEIF((A.GE.ZERO).AND.(B.GE.ZERO)) THEN TSSMIN
IF(C.GE.D) THEN TSSMIN
Q = D TSSMIN
BC = B TSSMIN
CALL TSLMIN(S1,SP,BC,Q,ANLS,FQ,ADT,AG,CONST1,CONST2, TSSMIN
+ DLT,NR,M,N,P,NWTAKE,IERR,TOL,VN,VNP, TSSMIN
+ VNS,OPTIM) TSSMIN
ELSE TSSMIN
OPTIM = SP TSSMIN
Q = C TSSMIN
ENDIF TSSMIN
ENDIF TSSMIN
TSSMIN
c update the global minimizer of FALPHA so far TSSMIN
TSSMIN
IF(FIRST) THEN TSSMIN
IF(A.GT.ZERO .AND. B.LT.ZERO) THEN TSSMIN
GLOPT = C TSSMIN
SOL = SP TSSMIN
IF(C.GT.D) THEN TSSMIN
GLOPT = D TSSMIN
SOL = S1 TSSMIN
ENDIF TSSMIN
ELSE TSSMIN
GLOPT = Q TSSMIN
SOL = OPTIM TSSMIN
ENDIF TSSMIN
FIRST = .FALSE. TSSMIN
ELSEIF(GLOPT.GE.Q) THEN TSSMIN
GLOPT = Q TSSMIN
SOL = OPTIM TSSMIN
ENDIF TSSMIN
10 CONTINUE TSSMIN
TSSMIN
RETURN TSSMIN
END TSSMIN
SUBROUTINE TSSMRD(VECT,RESNEW,X,MU,IERR,M,N) TSSMRD
TSSMRD
INTEGER M,N,IERR TSSMRD
DOUBLE PRECISION RESNEW,MU TSSMRD
DOUBLE PRECISION VECT(M),X(N) TSSMRD
TSSMRD
C********************************************************************** TSSMRD
C THIS ROUTINE COMPUTES THE RESIDUAL OF THE STANDARD MODEL. TSSMRD
C********************************************************************** TSSMRD
C TSSMRD
C INPUT PARAMETERS : TSSMRD
C ----------------- TSSMRD
C TSSMRD
C VECT : RIGHT HAND SIDE VECTOR OF THE NEWTON/GAUSS-NEWTON TSSMRD
C EQUATIONS AFTER BEING MULTIPLIED BY ORTHOGONAL MATRICES TSSMRD
C (SEE SUBROUTINE TSCPSS) TSSMRD
C X : STANDARD STEP COMPUTED BY THE SUBROUTINE TSCPSS TSSMRD
C MU : A SMALL PERTURBATION USED IN COMPUTING THE STANDARD TSSMRD
C STEP WHEN THE JACOBIAN IS SINGULAR TSSMRD
C IERR : RETURN CODE WITH THE FOLLOWING MEANINGS : TSSMRD
C IERR = 0 : NO SINGULARITY OF JACOBIAN DETECTED TSSMRD
C IERR = 1 : OTHERWISE TSSMRD
C M,N : DIMENSION OF PROBLEM TSSMRD
C TSSMRD
C OUTPUT PARAMETERS : TSSMRD
C ------------------ TSSMRD
C TSSMRD
C RESNEW : RESIDUAL OF THE STANDARD MODEL TSSMRD
C TSSMRD
C SUBPROGRAMS CALLED: TSSMRD
C TSSMRD
C LEVEL 1 BLAS ... DNRM2 TSSMRD
C TSSMRD
C********************************************************************** TSSMRD
TSSMRD
DOUBLE PRECISION TEMP,PROD TSSMRD
DOUBLE PRECISION DNRM2 TSSMRD
INTRINSIC SQRT TSSMRD
TSSMRD
IF(IERR .EQ .0) THEN TSSMRD
RESNEW = DNRM2(M-N,VECT(N+1),1) TSSMRD
ELSE TSSMRD
TEMP = DNRM2(M,VECT(N+1),1)**2 TSSMRD
PROD = MU * DNRM2(N,X,1)**2 TSSMRD
RESNEW = SQRT(TEMP-PROD) TSSMRD
ENDIF TSSMRD
TSSMRD
RETURN TSSMRD
END TSSMRD
SUBROUTINE TSSQP1(AJA,ANLS,S,F,MAXM,MAXN,M,N,Q,ROOT,RESTNS) TSSQP1
TSSQP1
INTEGER MAXM,MAXN,M,N,Q TSSQP1
DOUBLE PRECISION ROOT,RESTNS TSSQP1
DOUBLE PRECISION AJA(MAXM,N),S(MAXN,*),ANLS(MAXM,*),F(M) TSSQP1
TSSQP1
C********************************************************************** TSSQP1
C THIS ROUTINE SOLVES THE LOWER M-N+Q QUADRATIC EQUATIONS IN P UNKNOWNS TSSQP1
C OF THE TENSOR MODEL WHEN Q > 1 AND P = 1. TSSQP1
C********************************************************************** TSSQP1
C TSSQP1
C INPUT PARAMETERS : TSSQP1
C ----------------- TSSQP1
C TSSQP1
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSSQP1
C ANLS : TENSOR TERM MATRIX AT CURRENT ITERATE TSSQP1
C S : MATRIX OF PAST LINEARLY INDEPENDENT DIRECTIONS TSSQP1
C F : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY AN TSSQP1
C ORTHOGONAL MATRIX TSSQP1
C MAXM : LEADING DIMENSION OF AJA AND ANLS TSSQP1
C MAXN : LEADING DIMENSION OF S TSSQP1
C M,N : ROW AND COLUMN DIMENSIONS OF AJA TSSQP1
C Q : NUMERICAL RANK OF JACOBIAN : TSSQP1
C Q > P : JACOBIAN IS SINGULAR TSSQP1
C Q = P : OTHERWISE TSSQP1
C TSSQP1
C TSSQP1
C OUTPUT PARAMETERS : TSSQP1
C ----------------- TSSQP1
C TSSQP1
C ROOT : SOLUTION TO THE SYSTEM TSSQP1
C RESTNS : TENSOR RESIDUAL TSSQP1
C TSSQP1
C SUBPROGRAMS CALLED: TSSQP1
C TSSQP1
C LEVEL 1 BLAS ... DDOT,DNRM2 TSSQP1
C TSSQP1
C **********************************************************************TSSQP1
TSSQP1
INTEGER I TSSQP1
DOUBLE PRECISION TEMP,A,B,C,D,RES1,RES2,RES3,RES,S1,S2,S3 TSSQP1
DOUBLE PRECISION T,T0,T1,T2,T3,PI,A1,A2,A3,THETA TSSQP1
DOUBLE PRECISION ZERO,QUART,HALF,ONE,TWO,THREE,FOUR,NINE TSSQP1
DOUBLE PRECISION TSEVEN,FFOUR,ONETRD TSSQP1
DOUBLE PRECISION DDOT,DNRM2 TSSQP1
INTRINSIC ABS,SQRT TSSQP1
PARAMETER (PI = 3.141592653589793D0) TSSQP1
DATA ZERO,QUART,HALF,ONE,TWO,THREE,FOUR,NINE,TSEVEN,FFOUR/0.0D0, TSSQP1
+ 0.250D0,0.50D0,1.0D0,2.0D0,3.0D0,4.0D0,9.0D0,27.0D0,54.0D0/ TSSQP1
TSSQP1
c compute the coefficients of a third degree polynomial TSSQP1
TSSQP1
ONETRD = ONE/THREE TSSQP1
A = ZERO TSSQP1
B = ZERO TSSQP1
C = ZERO TSSQP1
TSSQP1
TEMP = DNRM2(M-N+Q,F(N-Q+1),1)**2 TSSQP1
D = TWO*DDOT(M-N+Q,AJA(N-Q+1,N),1,F(N-Q+1),1) TSSQP1
T0 = S(N+2,1)**2 TSSQP1
T1 = T0**2 TSSQP1
DO 10 I = N-Q+1,M TSSQP1
T2 = AJA(I,N) TSSQP1
T3 = ANLS(I,1) * T0 TSSQP1
C = C + TWO * (T2**2 + F(I) * T3) TSSQP1
B = B + THREE * T2 * T3 TSSQP1
A = A + ANLS(I,1)**2 * T1 TSSQP1
10 CONTINUE TSSQP1
TSSQP1
c compute the roots of the third degree polynomial TSSQP1
TSSQP1
IF(A.EQ.ZERO) THEN TSSQP1
IF(B.NE.ZERO) THEN TSSQP1
T0 = SQRT(C**2-FOUR*B*D) TSSQP1
T1 = TWO*B TSSQP1
S1 = (-C+T0)/T1 TSSQP1
S2 = (-C-T0)/T1 TSSQP1
RES1 = ABS(TEMP+D*S1+HALF*C*(S1**2)+ONETRD*B*(S1**3)) TSSQP1
RES2 = ABS(TEMP+D*S2+HALF*C*(S2**2)+ONETRD*B*(S2**3)) TSSQP1
IF(RES1 .GT. RES2) THEN TSSQP1
ROOT = S2 TSSQP1
RES = RES2 TSSQP1
ELSE TSSQP1
ROOT = S1 TSSQP1
RES = RES1 TSSQP1
ENDIF TSSQP1
RESTNS = SQRT(RES) TSSQP1
RETURN TSSQP1
ELSEIF(C.NE.ZERO) THEN TSSQP1
ROOT = -D/C TSSQP1
RES = ABS(TEMP+D*ROOT+HALF*C*(ROOT**2)) TSSQP1
RESTNS = SQRT(RES) TSSQP1
RETURN TSSQP1
ELSE TSSQP1
ROOT = ZERO TSSQP1
RESTNS = SQRT(TEMP) TSSQP1
RETURN TSSQP1
ENDIF TSSQP1
ELSEIF(D.EQ.ZERO) THEN TSSQP1
ROOT = ZERO TSSQP1
RESTNS = SQRT(TEMP) TSSQP1
RETURN TSSQP1
ENDIF TSSQP1
T3 = D TSSQP1
TSSQP1
A1 = B/A TSSQP1
A2 = C/A TSSQP1
A3 = D/A TSSQP1
T0 = (THREE*A2-A1**2)/NINE TSSQP1
T1 = (NINE*A1*A2-TSEVEN*A3-TWO*A1**3)/FFOUR TSSQP1
D = T0**3 + T1**2 TSSQP1
TSSQP1
IF(D.GT.0) THEN TSSQP1
T2 = T1+SQRT(D) TSSQP1
T = T1-SQRT(D) TSSQP1
IF(T.LT.0) THEN TSSQP1
T = -(-T)**ONETRD TSSQP1
ELSE TSSQP1
T = T**ONETRD TSSQP1
ENDIF TSSQP1
IF(T2.LT.0)THEN TSSQP1
T2 = -(-T2)**ONETRD TSSQP1
ELSE TSSQP1
T2 = T2**ONETRD TSSQP1
ENDIF TSSQP1
S1 = T2+T-A1/THREE TSSQP1
S3 = S1 TSSQP1
S2 = S1 TSSQP1
ELSE TSSQP1
T = T1/SQRT(-T0**3) TSSQP1
THETA = ACOS(T) TSSQP1
THETA = THETA/THREE TSSQP1
T = TWO*SQRT(-T0) TSSQP1
S1 = T*COS(THETA)-A1/THREE TSSQP1
S2 = T*COS(THETA+PI*TWO/THREE)-A1/THREE TSSQP1
S3 = T*COS(THETA+PI*FOUR/THREE)-A1/THREE TSSQP1
ENDIF TSSQP1
TSSQP1
c compute the tensor residual for each root TSSQP1
TSSQP1
RES1 = ABS(TEMP+T3*S1+HALF*C*(S1**2)+ONETRD*B*(S1**3)+ TSSQP1
+ QUART*A*(S1**4)) TSSQP1
RES2 = ABS(TEMP+T3*S2+HALF*C*(S2**2)+ONETRD*B*(S2**3)+ TSSQP1
+ QUART*A*(S2**4)) TSSQP1
RES3 = ABS(TEMP+T3*S3+HALF*C*(S3**2)+ONETRD*B*(S3**3)+ TSSQP1
+ QUART*A*(S3**4)) TSSQP1
TSSQP1
c select the root that produces the smallest tensor residual TSSQP1
TSSQP1
RES = RES1 TSSQP1
ROOT = S1 TSSQP1
IF(RES .GT. RES2) THEN TSSQP1
RES = RES2 TSSQP1
ROOT = S2 TSSQP1
ENDIF TSSQP1
IF(RES .GT. RES3) THEN TSSQP1
RES = RES3 TSSQP1
ROOT = S3 TSSQP1
ENDIF TSSQP1
RESTNS = SQRT(RES) TSSQP1
TSSQP1
RETURN TSSQP1
END TSSQP1
SUBROUTINE TSSSTP(AJA,FN,M,N,NR,EPSM,IGLOBL,AL2NRM,Y,W, TSSSTP
+ DN,FQ,PIVOT,PBAR,IERR) TSSSTP
TSSSTP
INTEGER NR,M,N,IERR,IGLOBL TSSSTP
INTEGER PIVOT(N),PBAR(N) TSSSTP
DOUBLE PRECISION EPSM,AJA(NR,N),AL2NRM(M),FN(M) TSSSTP
DOUBLE PRECISION DN(N),Y(N),W(M),FQ(M) TSSSTP
TSSSTP
C********************************************************************** TSSSTP
C THIS ROUTINE FINDS THE STANDARD STEP WHEN THE ITERATION NUMBER IS TSSSTP
C EQUAL TO 1 OR THE INPUT PARAMETER "METHOD" IS SET TO 0. TSSSTP
C********************************************************************** TSSSTP
C TSSSTP
C INPUT PARAMETERS : TSSSTP
C ----------------- TSSSTP
C TSSSTP
C AJA : JACOBIAN MATRIX AT CURRENT ITERATE TSSSTP
C FN : FUNCTION VALUE AT CURRENT ITERATE TSSSTP
C M,N : DIMENSIONS OF PROBLEM TSSSTP
C NR : LEADING DIMENSION OF AJA TSSSTP
C EPSM : MACHINE EPSILON TSSSTP
C IGLOBL: GLOBAL STRATEGY USED : TSSSTP
C = 0 : LINE SEARCH USED TSSSTP
C = 1 : 2-DIMENSIONAL TRUST REGION USED TSSSTP
C AL2NRM,Y,W : WORKING VECTORS TSSSTP
C TSSSTP
C TSSSTP
C OUTPUT PARAMETERS : TSSSTP
C ------------------ TSSSTP
C TSSSTP
C DN : STANDARD STEP TSSSTP
C FQ : FUNCTION VALUE AT CURRENT ITERATE MULTIPLIED BY TSSSTP
C ORTHOGONAL MATRICES (THIS IS USED IN THE CASE WHERE TSSSTP
C THE GLOBAL STRATEGY IS THE 2-DIMENSIONAL TSSSTP
C TRUST REGION) TSSSTP
C PIVOT,PBAR : PIVOT VECTORS TSSSTP
C IERR : RETURNED CODE WITH THE FOLLOWING MEANING : TSSSTP
C IERR = 1 : SINGULARITY OF JACOBIAN DETECTED (ZERO1 TSSSTP
C IS USED TO KEEP TRACK OF THE FIRST COLUMN TSSSTP
C WHERE SINGULARITY IS DETECTED) TSSSTP
C IERR = 0 : OTHERWISE TSSSTP
C TSSSTP
C SUBPROGRAMS CALLED: TSSSTP
C TSSSTP
C LEVEL 1 BLAS ... DCOPY,DLOAD,DSCAL TSSSTP
C TENSOLVE ... TSQRFC,TSQMUV,TSBSLV,TSPRMV,TSCPMU TSSSTP
C TSSSTP
C********************************************************************** TSSSTP
TSSSTP
INTEGER ZERO1,ZEROTM,I,J TSSSTP
DOUBLE PRECISION MU,ZERO,ONE TSSSTP
DATA ZERO,ONE/0.0D0,1.0D0/ TSSSTP
TSSSTP
c perform a QR factorization of AJA TSSSTP
TSSSTP
CALL TSQRFC(AJA,NR,N,M,1,N,IERR,EPSM,AL2NRM,PIVOT,ZERO1) TSSSTP
TSSSTP
CALL DSCAL(M,-ONE,FN,1) TSSSTP
TSSSTP
IF(IERR.EQ.0) THEN TSSSTP
IF(M.EQ.N) THEN TSSSTP
ZEROTM = ZERO1-1 TSSSTP
ELSE TSSSTP
ZEROTM = ZERO1 TSSSTP
ENDIF TSSSTP
TSSSTP
c multiply (-FN) by the orthogonal matrix resulting from the QR TSSSTP
c decomposition of AJA TSSSTP
TSSSTP
CALL TSQMUV(AJA,FN,W,NR,M,1,ZEROTM,.FALSE.) TSSSTP
TSSSTP
c solve AJA*DN = W TSSSTP
TSSSTP
CALL TSBSLV(AJA,NR,M,N,W,Y) TSSSTP
CALL TSPRMV(DN,Y,PIVOT,N,0) TSSSTP
TSSSTP
IF(IGLOBL.EQ.1) THEN TSSSTP
IERR = 0 TSSSTP
CALL DCOPY(M,W,1,FQ,1) TSSSTP
CALL DSCAL(M,-ONE,FQ,1) TSSSTP
ENDIF TSSSTP
RETURN TSSSTP
ELSE TSSSTP
TSSSTP
c AJA is singular TSSSTP
TSSSTP
CALL TSQMUV(AJA,FN,W,NR,M,1,ZERO1,.FALSE.) TSSSTP
TSSSTP
c solve ( AJA-trans AJA + MU I ) DN = -AJA-trans FN TSSSTP
TSSSTP
c put the diagonal elements stored in row m+2 of AJA into their TSSSTP
c propre positions and zero out the unwanted portions of AJA TSSSTP
TSSSTP
DO 10 J = 1, ZERO1-1 TSSSTP
AJA(J,J) = AJA(M+2,J) TSSSTP
CALL DLOAD (M+N-J,ZERO,AJA(J+1,J),1) TSSSTP
10 CONTINUE TSSSTP
TSSSTP
DO 20 J = ZERO1, N TSSSTP
CALL DLOAD (M+N-ZERO1+1,ZERO,AJA(ZERO1,J),1) TSSSTP
20 CONTINUE TSSSTP
TSSSTP
c compute a perturbation MU TSSSTP
TSSSTP
CALL TSCPMU(AJA,NR,N,EPSM,MU) TSSSTP
TSSSTP
c form the augmented Jacobian matrix by adding an nxn diag(mu) in TSSSTP
c the bottom of AJA TSSSTP
TSSSTP
DO 70 I = M+1,M+N TSSSTP
AJA(I,I-M) = MU TSSSTP
70 CONTINUE TSSSTP
TSSSTP
c factorize the transformed matrix AJA from 1 to n and compute TSSSTP
c the standard step DN TSSSTP
TSSSTP
CALL TSQRFC(AJA,NR,N,M+N,1,N,IERR,EPSM,AL2NRM,PBAR,ZERO1) TSSSTP
CALL TSQMUV(AJA,W,FQ,NR,M+N,1,N+1,.FALSE.) TSSSTP
CALL TSBSLV(AJA,NR,M+N,N,FQ,DN) TSSSTP
CALL TSPRMV(Y,DN,PBAR,N,0) TSSSTP
CALL TSPRMV(DN,Y,PIVOT,N,0) TSSSTP
ENDIF TSSSTP
TSSSTP
IF(IGLOBL.EQ.1) THEN TSSSTP
IERR = 1 TSSSTP
CALL DSCAL(M+N,-ONE,FQ,1) TSSSTP
ENDIF TSSSTP
TSSSTP
END TSSSTP
SUBROUTINE TSSTMX(S,X,NR,N,P,WRK1,WRK2) TSSTMX
INTEGER NR,N,P TSSTMX
DOUBLE PRECISION X(P),S(NR,P),WRK1(P),WRK2(P) TSSTMX
TSSTMX
C********************************************************************* TSSTMX
C THIS ROUTINE COMPUTES SHAT-TRANS * X, WHERE SHAT IS AN UPSIDE DOWN TSSTMX
C TRIANGULAR MATRIX RESULTING FROM A QL FACTORIZATION OF A MATRIX TSSTMX
C A AND X IS A VECTOR. TSSTMX
C********************************************************************* TSSTMX
C TSSTMX
C INPUT PARAMETERS : TSSTMX
C ----------------- TSSTMX
C TSSTMX
C SHAT : UPSIDE DOWN TRIANGULAR MATRIX RESULTING FROM A QL TSSTMX
C FACTORIZATION TSSTMX
C X : VECTOR TO BE MULTIPLIED BY SHAT TSSTMX
C NR : LEADING DIMENSION OF SHAT TSSTMX
C N : ROW DIMENSION OF THE MATRIX A TSSTMX
C P : COLUMN DIMENSION OF SHAT TSSTMX
C WRK : WORKSPACE TSSTMX
C TSSTMX
C TSSTMX
C OUTPUT PARAMETERS : TSSTMX
C ----------------- TSSTMX
C TSSTMX
C WRK2 : SHAT-TRANS * X TSSTMX
C TSSTMX
C SUBPROGRAMS CALLED: TSSTMX
C TSSTMX
C LEVEL 1 BLAS ... DCOPY,DDOT,DLOAD TSSTMX
C TSSTMX
C********************************************************************* TSSTMX
TSSTMX
INTEGER COL TSSTMX
DOUBLE PRECISION ZERO TSSTMX
DOUBLE PRECISION DDOT TSSTMX
DATA ZERO/0.0D0/ TSSTMX
TSSTMX
CALL DLOAD(P,ZERO,WRK1,1) TSSTMX
TSSTMX
WRK2(1) = S(N+2,1) * X(P) TSSTMX
IF(P .GT. 1) THEN TSSTMX
WRK1(P) = S(N,2) TSSTMX
WRK1(P-1) = S(N+2,2) TSSTMX
WRK2(2) = DDOT(P,WRK1,1,X,1) TSSTMX
DO 10 COL = 3, P TSSTMX
CALL DCOPY(COL-1,S(N-COL+2,COL),1,WRK1(P-COL+2),1) TSSTMX
WRK1(P-COL+1) = S(N+2,COL) TSSTMX
WRK2(COL) = DDOT(P,WRK1,1,X,1) TSSTMX
10 CONTINUE TSSTMX
ENDIF TSSTMX
TSSTMX
RETURN TSSTMX
END TSSTMX
SUBROUTINE TSTRUD(M,N,X,F,G,SC,NWTAKE,STEPMX,STEPTL,DLT,MXTAKE, TSTRUD
+ DXN,DFN,TSFVEC,SCRES,IRETCD,XPLSP,FPLSP,FPREV, TSTRUD
+ XPLS,FP,FPLS) TSTRUD
TSTRUD
INTEGER M,N,IRETCD TSTRUD
DOUBLE PRECISION F,STEPMX,STEPTL,DLT,SCRES,FPLSP,FPLS TSTRUD
DOUBLE PRECISION X(N),XPLS(N),G(N) TSTRUD
DOUBLE PRECISION SC(N),XPLSP(N),FPREV(M),FP(M) TSTRUD
DOUBLE PRECISION DXN(N),DFN(M) TSTRUD
LOGICAL NWTAKE,MXTAKE TSTRUD
EXTERNAL TSFVEC TSTRUD
TSTRUD
C***********************************************************************TSTRUD
C THIS ROUTINE DECIDES WHETHER TO ACCEPT XPLS=X+SC AS THE NEXT ITERATE TSTRUD
C AND UPDATES THE TRUST REGION DLT. TSTRUD
C***********************************************************************TSTRUD
C TSTRUD
C TSTRUD
C TSTRUD
C PARAMETERS TSTRUD
C ---------- TSTRUD
C M,N --> DIMENSIONS OF PROBLEM TSTRUD
C X(N) --> OLD ITERATE X[K-1] TSTRUD
C F --> 0.50D0 * || FC ||**2 TSTRUD
C G(N) --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE TSTRUD
C SC(N) --> CURRENT STEP TSTRUD
C NWTAKE --> BOOLEAN, =.TRUE. IF INPUT STEP TAKEN TSTRUD
C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE TSTRUD
C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES TSTRUD
C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM TSTRUD
C DLT <--> TRUST REGION RADIUS TSTRUD
C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED TSTRUD
C DXN --->DIAGONAL SCALING MATRIX FOR X TSTRUD
C DFN --->DIAGONAL SCALING MATRIX FOR F TSTRUD
C TSFVEC --->SUBROUTINE TO EVALUATE FUNCTION TSTRUD
C TSTRUD
C IRETCD <--> RETURN CODE TSTRUD
C =0 XPLS ACCEPTED AS NEXT ITERATE; TSTRUD
C DLT TRUST REGION FOR NEXT ITERATION. TSTRUD
C =1 XPLS UNSATISFACTORY BUT ACCEPTED AS NEXT ITERATETSTRUD
C BECAUSE XPLS-X .LT. SMALLEST ALLOWABLE TSTRUD
C STEP LENGTH. TSTRUD
C =2 F(XPLS) TOO LARGE. CONTINUE CURRENT ITERATION TSTRUD
C WITH NEW REDUCED DLT. TSTRUD
C =3 F(XPLS) SUFFICIENTLY SMALL, BUT QUADRATIC MODEL TSTRUD
C PREDICTS F(XPLS) SUFFICIENTLY WELL TO CONTINUE TSTRUD
C CURRENT ITERATION WITH NEW DOUBLED DLT. TSTRUD
C XPLSP(N) <--> WORKSPACE [VALUE NEEDS TO BE RETAINED BETWEEN TSTRUD
C SUCCESSIVE CALLS OF K-TH GLOBAL STEP] TSTRUD
C FPLSP <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] TSTRUD
C FPREV ---> WORKING VECTOR TSTRUD
C XPLS(N) <-- NEW ITERATE X[K] TSTRUD
C FP <-- FUNCTION VALUE AT NEXT ITERATE TSTRUD
C FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS) TSTRUD
C TSTRUD
C SUBPROGRAMS CALLED: TSTRUD
C TSTRUD
C LEVEL 1 BLAS ... DCOPY,DDOT,DNRM2 TSTRUD
C TENSOLVE ... TSFSCL TSTRUD
C TSTRUD
C********************************************************************** TSTRUD
TSTRUD
INTEGER I TSTRUD
DOUBLE PRECISION STEPLN,DLTFP,SLOPE,DLTF,SLP,PQ,RLN,DLTMP TSTRUD
DOUBLE PRECISION DNRM2,DDOT TSTRUD
INTRINSIC ABS,MAX TSTRUD
DOUBLE PRECISION ZERO,TENTH,HALF,ZNN,ONE,TWO TSTRUD
DATA ZERO,TENTH,HALF,ZNN,ONE,TWO/0.0D0,0.10D0,0.50D0, TSTRUD
+ 0.99D0,1.0D0,2.0D0/ TSTRUD
TSTRUD
MXTAKE = .FALSE. TSTRUD
DO 90 I = 1,N TSTRUD
XPLS(I) = X(I)+SC(I) TSTRUD
90 CONTINUE TSTRUD
STEPLN = DNRM2(N,SC,1) TSTRUD
CALL TSFSCL(XPLS,DXN,DFN,M,N,TSFVEC,FP) TSTRUD
FPLS = HALF*DNRM2(M,FP,1)**2 TSTRUD
DLTF = FPLS-F TSTRUD
SLOPE = DDOT(N,G,1,SC,1) TSTRUD
SLP = HALF*SCRES**2-F TSTRUD
TSTRUD
c next statement added for case of compilers which do not optimize TSTRUD
c evaluation of next "IF" statement (in which case fplsp could be TSTRUD
c undefined) TSTRUD
c TSTRUD
IF(IRETCD.EQ.4) FPLSP=0.0D0 TSTRUD
C$ WRITE(*,961) IRETCD,FPLS,FPLSP,DLTF,SLP TSTRUD
IF(IRETCD.NE.3 .OR. (FPLS.LT.FPLSP .AND. DLTF.LE. 1.D-4*SLP)) TSTRUD
+ GO TO 130 TSTRUD
C IF(IRETCD.EQ.3 .AND. (FPLS.GE.FPLSP .OR. DLTF.GT. 1.D-4*SLP)) TSTRUD
C THEN TSTRUD
C TSTRUD
C reset XPLS to XPLSP and terminate global step TSTRUD
C TSTRUD
IRETCD = 0 TSTRUD
CALL DCOPY(N,XPLSP,1,XPLS,1) TSTRUD
FPLS = FPLSP TSTRUD
DLT = HALF*DLT TSTRUD
CALL DCOPY(M,FPREV,1,FP,1) TSTRUD
C$ WRITE(*,951) TSTRUD
GO TO 230 TSTRUD
C ELSE TSTRUD
C TSTRUD
C FPLS too large TSTRUD
C TSTRUD
130 IF(DLTF.LE. 1.D-4*SLP) GO TO 170 TSTRUD
C IF(DLTF.GT. 1.D-4*SLP) TSTRUD
C THEN TSTRUD
C$ WRITE(*,952) TSTRUD
PQ = ONE TSTRUD
RLN = ZERO TSTRUD
DO 140 I = 1,N TSTRUD
RLN = MAX(RLN,ABS(SC(I))/MAX(ABS(XPLS(I)),ONE/PQ)) TSTRUD
140 CONTINUE TSTRUD
C$ WRITE(*,962) RLN TSTRUD
IF(RLN.GE.STEPTL) GO TO 150 TSTRUD
C IF(RLN.LT.STEPTL) TSTRUD
C THEN TSTRUD
C TSTRUD
C cannot find satisfactory XPLS sufficiently distinct from X TSTRUD
C TSTRUD
IRETCD = 1 TSTRUD
C$ WRITE(*,954) TSTRUD
GO TO 230 TSTRUD
C ELSE TSTRUD
C TSTRUD
C reduce trust region and continue global step TSTRUD
C TSTRUD
150 IRETCD = 2 TSTRUD
DLTMP = -SLOPE*STEPLN/(TWO*(DLTF-SLOPE)) TSTRUD
C$ WRITE(*,963) DLTMP TSTRUD
IF(DLTMP.GE. TENTH*DLT) GO TO 155 TSTRUD
C IF(DLTMP.LT. TENTH*DLT) TSTRUD
C THEN TSTRUD
DLT = TENTH*DLT TSTRUD
GO TO 160 TSTRUD
C ELSE TSTRUD
155 IF(DLTMP.GT.HALF*DLT) THEN TSTRUD
DLT = HALF*DLT TSTRUD
ELSE TSTRUD
DLT = DLTMP TSTRUD
ENDIF TSTRUD
C ENDIF TSTRUD
160 CONTINUE TSTRUD
C$ WRITE(*,955) TSTRUD
GO TO 230 TSTRUD
C ENDIF TSTRUD
C ELSE TSTRUD
C TSTRUD
C FPLS sufficiently small TSTRUD
C TSTRUD
170 CONTINUE TSTRUD
C$ WRITE(*,958) TSTRUD
DLTFP = HALF*SCRES**2-F TSTRUD
C$ WRITE(*,964) DLTFP,NWTAKE TSTRUD
IF(IRETCD.EQ.2 .OR. ((ABS(DLTFP-DLTF).GT. TENTH*ABS(DLTF)) TSTRUD
+ .AND. (DLTFP.GT.SLOPE)).OR.NWTAKE TSTRUD
+ .OR. (DLT.GT. ZNN*STEPMX)) GO TO 210 TSTRUD
C IF(IRETCD.NE.2 .AND. (ABS(DLTFP-DLTF) .LE. .1*ABS(DLTF)) TSTRUD
C + .AND. (.NOT.NWTAKE) .AND. (DLT.LE. .99*STEPMX)) TSTRUD
C THEN TSTRUD
C TSTRUD
C double trust region and continue global step TSTRUD
C TSTRUD
IRETCD = 3 TSTRUD
CALL DCOPY(N,XPLS,1,XPLSP,1) TSTRUD
FPLSP = FPLS TSTRUD
DLT = MIN(TWO*DLT,STEPMX) TSTRUD
CALL DCOPY(M,FP,1,FPREV,1) TSTRUD
C$ WRITE(*,959) TSTRUD
GO TO 230 TSTRUD
C ELSE TSTRUD
C TSTRUD
C accept XPLS as next iterate. Choose new trust region. TSTRUD
C TSTRUD
210 CONTINUE TSTRUD
C$ WRITE(*,960) TSTRUD
IRETCD = 0 TSTRUD
IF(DLT.GT. ZNN*STEPMX) MXTAKE = .TRUE. TSTRUD
IF(DLTF.LT. TENTH*DLTFP) GO TO 220 TSTRUD
C IF(DLTF.GE. TENTH*DLTFP) TSTRUD
C THEN TSTRUD
C TSTRUD
C Decrease trust region for next iteration TSTRUD
C TSTRUD
DLT = HALF*DLT TSTRUD
GO TO 230 TSTRUD
C ELSE TSTRUD
C Check whether to increase trust region for next iteration TSTRUD
C TSTRUD
220 IF(DLTF.LE. .75D0*DLTFP) DLT = MIN(TWO*DLT,STEPMX) TSTRUD
C ENDIF TSTRUD
C ENDIF TSTRUD
C ENDIF TSTRUD
C ENDIF TSTRUD
230 CONTINUE TSTRUD
C$ WRITE(*,953) TSTRUD
C$ WRITE(*,956) IRETCD,MXTAKE,DLT,FPLS TSTRUD
C$ WRITE(*,957) TSTRUD
C$ WRITE(*,965) (XPLS(I),I = 1,N) TSTRUD
RETURN TSTRUD
C TSTRUD
951 FORMAT(' TSTRUD RESET XPLS TO XPLSP. TERMINATION GLOBAL STEP') TSTRUD
952 FORMAT(' TSTRUD FPLS TOO LARGE') TSTRUD
953 FORMAT(' TSTRUD VALUES AFTER CALL TO TREGUP') TSTRUD
954 FORMAT(' TSTRUD CANNOT FIND SATISFACTORY XPLS DISTINCT FROM', TSTRUD
+ ' X. TERMINATE GLOBAL STEP') TSTRUD
955 FORMAT(' TSTRUD REDUCE TRUST REGION. CONTINUE GLOBAL STEP') TSTRUD
956 FORMAT(' TSTRUD IRETCD=',I3/ TSTRUD
+ ' TSTRUD MXTAKE=',L1/ TSTRUD
+ ' TSTRUD DLT =',E20.13/ TSTRUD
+ ' TSTRUD FPLS =',E20.13) TSTRUD
957 FORMAT(' TSTRUD NEW ITERATE (XPLS)') TSTRUD
958 FORMAT(' TSTRUD FPLS SUFFICIENTLY SMALL') TSTRUD
959 FORMAT(' TSTRUD DOUBLE TRUST REGION. CONTINUE GLOBAL STEP') TSTRUD
960 FORMAT(' TSTRUD ACCEPT XPLS AS NEW ITERATE. CHOOSE NEW', TSTRUD
+ ' TRUST REGION. TERMINATE GLOBAL STEP') TSTRUD
961 FORMAT(' TSTRUD IRETCD=',I5/ TSTRUD
+ ' TSTRUD FPLS =',E20.13/ TSTRUD
+ ' TSTRUD FPLSP =',E20.13/ TSTRUD
+ ' TSTRUD DLTF =',E20.13/ TSTRUD
+ ' TSTRUD SLP =',E20.13) TSTRUD
962 FORMAT(' TSTRUD RLN =',E20.13) TSTRUD
963 FORMAT(' TSTRUD DLTMP =',E20.13) TSTRUD
964 FORMAT(' TSTRUD DLTFP =',E20.13/ TSTRUD
+ ' TSTRUD NWTAKE=',L1) TSTRUD
965 FORMAT(' TSTRUD ',5(E20.13,3X)) TSTRUD
END TSTRUD
SUBROUTINE TSUDQV(SHAT,WRK1,NR,N,P,CONST1) TSUDQV
TSUDQV
INTEGER NR,N,P TSUDQV
DOUBLE PRECISION CONST1(P),SHAT(NR,P),WRK1(N) TSUDQV
TSUDQV
C********************************************************************** TSUDQV
C THIS ROUTINE COMPUTES SHAT-TRANS * WRK1, WHERE SHAT IS AN UPSIDE TSUDQV
C DOWN TRIANGULAR MATRIX RESULTING FROM A QL FACTORIZATION OF A TSUDQV
C MATRIX A AND WRK1 IS A VECTOR OF LENGTH N. TSUDQV
C********************************************************************** TSUDQV
C TSUDQV
C INPUT PARAMETERS TSUDQV
C ---------------- TSUDQV
C TSUDQV
C SHAT : UPSIDE DOWN TRIANGULAR MATRIX RESULTING FROM A QL TSUDQV
C FACTORIZATION TSUDQV
C WRK1 : VECTOR TO BE MULTIPLIED BY SHAT TSUDQV
C NR : LEADING DIMENSION OF SHAT TSUDQV
C N : DIMENSION OF MATRIX A TSUDQV
C P : COLUMN DIMENSION OF SHAT TSUDQV
C TSUDQV
C OUTPUT PARAMETERS TSUDQV
C ----------------- TSUDQV
C TSUDQV
C CONST1 : SHAT * WRK1 TSUDQV
C TSUDQV
C SUBPROGRAMS CALLED: TSUDQV
C TSUDQV
C LEVEL 1 BLAS ... DDOT TSUDQV
C TSUDQV
C **********************************************************************TSUDQV
TSUDQV
INTEGER J TSUDQV
DOUBLE PRECISION DDOT TSUDQV
TSUDQV
CONST1(1) = SHAT(N+2,1) * WRK1(N) TSUDQV
IF(P .GT. 1) THEN TSUDQV
CONST1(2) = SHAT(N,2) * WRK1(N) + SHAT(N+2,2) * WRK1(N-1) TSUDQV
ENDIF TSUDQV
DO 20 J = 3,P TSUDQV
CONST1(J) = DDOT(J-1,SHAT(N-J+2,J),1,WRK1(N-J+2),1) TSUDQV
+ + SHAT(N+2,J) * WRK1(N-J+1) TSUDQV
20 CONTINUE TSUDQV
TSUDQV
RETURN TSUDQV
END TSUDQV
SUBROUTINE TSUNSF(F,DF,M) TSUNSF
TSUNSF
INTEGER M TSUNSF
DOUBLE PRECISION F(M),DF(M) TSUNSF
TSUNSF
C********************************************************************* TSUNSF
C THIS ROUTINE UNSCALES A FUNCTION VALUE F. TSUNSF
C********************************************************************* TSUNSF
C TSUNSF
C INPUT PARAMETERS : TSUNSF
C ------------------ TSUNSF
C TSUNSF
C DF : DIAGONAL SCALING MATRIX FOR F TSUNSF
C M : DIMENSION OF F TSUNSF
C TSUNSF
C INPUT-OUTPUT PARAMETERS : TSUNSF
C ------------------ TSUNSF
C TSUNSF
C F : SCALED FUNCTION VALUE ON ENTRY AND UNSCALED FUNCTION TSUNSF
C VALUE ON EXIT TSUNSF
C TSUNSF
C********************************************************************** TSUNSF
TSUNSF
INTEGER I TSUNSF
TSUNSF
DO 10 I = 1,M TSUNSF
F(I) = F(I)/DF(I) TSUNSF
10 CONTINUE TSUNSF
TSUNSF
RETURN TSUNSF
END TSUNSF
SUBROUTINE TSUNSX(X,DX,N) TSUNSX
TSUNSX
INTEGER N TSUNSX
DOUBLE PRECISION X(N),DX(N) TSUNSX
TSUNSX
C********************************************************************** TSUNSX
C THIS ROUTINE UNSCALES A VECTOR X. TSUNSX
C********************************************************************** TSUNSX
C TSUNSX
C INPUT PARAMETERS : TSUNSX
C ------------------ TSUNSX
C TSUNSX
C DX : DIAGONAL SCALING MATRIX FOR X TSUNSX
C N : DIMENSION OF X TSUNSX
C TSUNSX
C OUTPUT PARAMETERS : TSUNSX
C ------------------ TSUNSX
C TSUNSX
C X : UNSCALED VECTOR X TSUNSX
C TSUNSX
C********************************************************************** TSUNSX
TSUNSX
INTEGER I TSUNSX
TSUNSX
DO 10 I = 1,N TSUNSX
X(I) = X(I)/DX(I) TSUNSX
10 CONTINUE TSUNSX
TSUNSX
RETURN TSUNSX
END TSUNSX
SUBROUTINE TSUPSF(FC,XC,XP,SQRN,ITN,MAXM,MAXN,M,N,STEP,S,FV) TSUPSF
TSUPSF
INTEGER MAXM,MAXN,M,N,ITN,SQRN TSUPSF
DOUBLE PRECISION S(MAXN,*),FV(MAXM,*) TSUPSF
DOUBLE PRECISION FC(M),XC(N),XP(N),STEP(N) TSUPSF
TSUPSF
C********************************************************************** TSUPSF
C THIS ROUTINE UPDATE THE MATRIX S OF PAST DIRECTIONS AND THE MATRIX TSUPSF
C FV OF FUNCTION VALUES. TSUPSF
C********************************************************************** TSUPSF
C TSUPSF
C INPUT PARAMETERS : TSUPSF
C ---------------- TSUPSF
C TSUPSF
C FC : FUNCTION VALUE AT CURRENT ITERATE TSUPSF
C XC : CURRENT ITERATE X[K-1] TSUPSF
C XP : NEW ITERATE X[K] TSUPSF
C SQRN: MAXIMUM COLUMN DIMENSION OF S AND FV TSUPSF
C ITN : ITERATION NUMBER TSUPSF
C MAXM: LEADING DIMENSION OF FV TSUPSF
C MAXN: LEADING DIMENSION OF S TSUPSF
C M : ROW DIMENSION OF MATRIX FV TSUPSF
C N : ROW DIMENSION OF MATRIX S TSUPSF
C STEP: WORKING VECTOR TSUPSF
C TSUPSF
C TSUPSF
C INPUT-OUTPUT PARAMETERS : TSUPSF
C ----------------------- TSUPSF
C TSUPSF
C S : MATRIX OF PAST DIRECTIONS (XK - XC) TSUPSF
C FV : MATRIX OF PAST FUNCTIONS VALUES TSUPSF
C TSUPSF
C SUBPROGRAMS CALLED: TSUPSF
C TSUPSF
C LEVEL 1 BLAS ... DCOPY TSUPSF
C TSUPSF
C********************************************************************** TSUPSF
TSUPSF
INTEGER I,J,L TSUPSF
TSUPSF
c update FV TSUPSF
TSUPSF
IF(SQRN.LT.ITN) THEN TSUPSF
L = SQRN TSUPSF
ELSE TSUPSF
L = ITN TSUPSF
ENDIF TSUPSF
DO 10 J = L-1,1,-1 TSUPSF
CALL DCOPY(M,FV(1,J),1,FV(1,J+1),1) TSUPSF
10 CONTINUE TSUPSF
TSUPSF
CALL DCOPY(M,FC,1,FV(1,1),1) TSUPSF
TSUPSF
c update S TSUPSF
TSUPSF
DO 30 I = 1,N TSUPSF
STEP(I) = XC(I)-XP(I) TSUPSF
30 CONTINUE TSUPSF
TSUPSF
DO 50 J = L-1,1,-1 TSUPSF
DO 40 I = 1,N TSUPSF
S(I,J+1) = S(I,J) + STEP(I) TSUPSF
40 CONTINUE TSUPSF
50 CONTINUE TSUPSF
CALL DCOPY(N,STEP,1,S(1,1),1) TSUPSF
TSUPSF
RETURN TSUPSF
END TSUPSF
SUBROUTINE TSUTMD(AJA,D,NR,M,N,V) TSUTMD
TSUTMD
INTEGER NR,M,N TSUTMD
DOUBLE PRECISION AJA(NR,N),D(N),V(N) TSUTMD
TSUTMD
C********************************************************************** TSUTMD
C THIS ROUTINE MULTIPLIES AN UPPER TRIANGULAR MATRIX (AS STORED IN TSUTMD
C STEWART) BY A VECTOR D. TSUTMD
C********************************************************************** TSUTMD
C TSUTMD
C INPUT PARAMETERS : TSUTMD
C ----------------- TSUTMD
C TSUTMD
C AJA : JACOBIAN AT CURRENT ITERATE TSUTMD
C D : VECTOR TO BE MULTIPLIED BY AJA TSUTMD
C NR : LEADING DIMENSION OF AJA TSUTMD
C M,N : DIMENSIONS OF PROBLEM TSUTMD
C TSUTMD
C OUTPUT PARAMETERS : TSUTMD
C ----------------- TSUTMD
C TSUTMD
C V : VECTOR AJA * D ON EXIT TSUTMD
C TSUTMD
C SUBPROGRAMS CALLED: TSUTMD
C TSUTMD
C LEVEL 1 BLAS ... DAXPY TSUTMD
C TSUTMD
C********************************************************************** TSUTMD
TSUTMD
INTEGER J TSUTMD
TSUTMD
V(1) = AJA(M+2,1) * D(1) + AJA(1,2) * D(2) TSUTMD
V(2) = AJA(M+2,2) * D(2) TSUTMD
DO 20 J = 3, N TSUTMD
V(J) = AJA(M+2,J) * D(J) TSUTMD
CALL DAXPY(J-1,D(J),AJA(1,J),1,V,1) TSUTMD
20 CONTINUE TSUTMD
TSUTMD
RETURN TSUTMD
END TSUTMD
|
|
IP卡
狗仔卡
发表于 2003-6-10 05:51:41
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
