C234567==1=========2=========3=========4=========5=========6=========7=
C
C                       MANPAK
C                     Version 2                 June 10, 1996
C
C   This is a package of utility programs for computations with
C   sub-manifolds of R^NVAR that are implicitly defined by a system
C   of nonlinear equations
C
C               F(x) = 0.
C
C   Here F is a mapping from R^NVAR to R^NALG, NALG < NVAR, which
C   is sufficiently smooth on an open set E in R^NVAR and satisfies
C
C      rank DF(x) = NALG for all x in M = { u in E;  F(u) = 0 }
C
C   Then M is a submanifold of R^NVAR of dimension MDIM = NVAR-NALG.
C   This setting is assumed for all routines. The first and second
C   derivative of F are denoted by DF and D2F, respectively.
C
C   The routines in the package establish several types of local
C   parametrizations (coordinate systems) and -- once available --
C   compute points on the manifold with given local coordinates.
C   In addition, there are routines for computing the first and
C   second derivatives of these local parametrizations and some
C   other quantities, such as sensitivity measures and the second
C   fundamental tensor.
C
C   For a discussion of the algorithms see
C
C      W. C. Rheinboldt, MANPAK: A Set of Algorithms for
C      Computations on Implicitly Defined Manifolds
C      Inst. for Comp. Math. and Appl., Univ. of Pittsburgh,
C      Tech. Report. TR-ICMA-96-198, July 1996
C      J. Comp. and Math. Applic. submitted
C
C   The routines use reverse communication for calls to user functions;
C   that is, when needed they request from the user the evaluation of
C   the function F or its derivatives at a given point. For details see
C   the individual routines.
C
C   TWO different data structures are used for the local
C   parametrizations called TANGENTIAL and GENERAL, respectively.
C
C   In the TANGENTAL case the data structure consists of a block
C   of five arrays:
C
C      X(NVAR),UBX(NVAR,MDIM),DFX(NALG,NVAR),TAUX(MDIM),JPX(MDIM),
C
C   Here X is the center point of the coordinate system, DFX contains
C   the LQ-factorization (with row pivoting) of the Jacobian DF(X)
C   and the columns of UBX form an orthonormal basis of the tangent
C   space at X. The array TAUX contains the diagonal terms of the
C   Householder transformations used in the LQ-factorization,
C   and JPX provides the corresponding pivot sequence.
C
C   In the GENERAL case the data structure consists of the block
C   of four arrays:
C
C     X(NVAR), UBX(NVAR,MDIM), AUGMT(NVAR,NVAR), JPAUG(MDIM)
C
C   Here X is the center of the coordinate system, the columns of
C   the array UBX form an orthonormal basis of the local coordinate
C   space, and the array AUGMT contains the LU-factorization of the
C   NVAR x NVAR augmented matrix
C
C                 ( DF(X) )
C         AUGMT = (       )
C                 ( UBX^T )
C
C   where DF(X) is the Jacobian of F at X. Finally, JPAUG is the
C   pivot array of the LU factorization.
C
C   The routines use the support routines in the package MANAUX.
C   All error information is handled by the routine MSGPRT included
C   in MANUAX
C
C   Routines in the package
C   -----------------------
C   SUBROUTINE AUGM  ...routine for generating the augmented matrix
C                       AUGMT and computing its LU factorization.
C   SUBROUTINE COBAS ...general-purpose routine for computing an
C                       orthonormal basis matrix UBX as the basis
C                       of the nullspace (kernel) of a given
C                       NALG x NVAR matrix with rank = NALG
C   SUBROUTINE CURVT ...for computing approximations of the curvature
C                       and principal normal of a path on the
C                       manifold passing through three consecutive
C                       points. This also provides an approximatation
C                       of the diagonal terms of the second fundamental
C                       tensor.
C   SUBROUTINE DGPHI ...for the first derivative of the local
C                       parametrization mapping of a general local
C                       coordinate system.
C   SUBROUTINE D2GPHI ..for the second derivative of the local
C                       parametrization mapping of a general local
C                       coordinate system.
C   SUBROUTINE DTPHI ...for differentiating the local parametrization
C                       mapping of a tangential local coordinate system
C   SUBROUTINE GCBAS ...for computing at X an orthonormal basis
C                       matrix UBX of a general local coordinate
C                       space that consists of MDIM canonical basis
C                       vectors of R^NVAR.
C   SUBROUTINE GLOB  ...for globalizing a vector Y of local
C                       coordinates; that is, for computing the
C                       point X = XC + UBXC*Y in the affine space
C                       XC + span(UBXC)
C   SUBROUTINE GNBAS ...for computing at X an orthonormal basis
C                       matrix UBX of a general local coordinate
C                       space that contains the NVAR-th canonical
C                       basis vector
C   SUBROUTINE GPHI  ...for computing a point on M with specified
C                       local coordinates in a general local
C                       coordinate system.
C   SUBROUTINE MOVFR ...uses a moving frame algorithm for ordering
C                       a local coordinate basis UBX2 such that it
C                       has the same orientation as another such
C                       basis UBX1
C   SUBROUTINE ORIENT...uses a combinatorial algorithm for ordering
C                       a local coordinate basis UBX2 such that it
C                       has the same orientation as another such
C                       basis UBX1,
C   SUBROUTINE PROJ  ...for computing the point Y = UBXC^T(X - XC)
C                       for a given global vector X; that is, the
C                       orthogonal projection of X - XC onto
C                       XC + span(UBXC)
C   SUBROUTINE SENMAP...for computing the sensitivity map at a
C                       given point with respect to specified natural
C                       parameters
C   SUBROUTINE SENSNR...for computing the Euclidean norm of the
C                       sensitivity map at a given point with
C                       respect to specified natural parameters
C   SUBROUTINE TPHI  ...for computing a point on M with specified
C                       local coordinates in a tangential local
C                       coordinate system
C   SUBROUTINE TSFT  ...for computing a component of the second
C                       fundamental tensor in a tangential local
C                       coordinate system.
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE AUGM( NVAR,MDIM,AUGMT,LDA,UBXC,LDU,JPAUG,IER )
C
      INTEGER NVAR,MDIM,LDA,LDU,IER,JPAUG(*)
      DOUBLE PRECISION AUGMT(LDA,*),UBXC(LDU,*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Call the routine with the Jacobian DF(X) stored in the first
C     NALG = NVAR - MDIM rows of the array AUGMT and an orthonormal
C     basis matrix of a local parametrization at a point XC stored
C     in the array UBXC.
C
C     The routine sets up the augmented matrix
C
C                      ( DF(X)  )
C           AUGMT  =   (        )
C            `         ( UBXC^T )
C
C     and computes its LU factorization.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   Dimension of the ambient space
C     MDIM   I   IN   Dimension of the manifold
C     AUGMT  D   IN   Array of dimension LDA x NVAR containing in its
C                     first NALG rows the Jacobian at a point X.
C                OUT  The triangular factors of the augmented matrix.
C     LDA    I   IN   Leading dimension of AUGMT, LDA >= NVAR
C     UBXC   D   IN   Array of dimension NVAR x MDIM, the basis  of
C                     a local coordinate system at XC.
C     LDU    I   IN   Leading dimension dimension of UBXC
C     JPAUG  I   OUT  Array of dimension NVAR, the pivot array of the
C                     factorization of AUGMT
C     IER    I   IN   error indicator
C                     IER = 0 : no error
C                     IER = -1: data error reported by LUF
C                     IER = 1 : numerically singular matrix
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....External subroutines called
C
      EXTERNAL LUF,MSGPRT
C
C.....Local variables
C
      INTEGER I,J,K,NALG
C
C.......................Executable statements..........................
C
C.....Store the coordinate basis UBXC in the last MDIM rows of AUGMT
C
      NALG = NVAR - MDIM
      DO 20 J = 1,MDIM
         K = NALG + J
         DO 10 I = 1,NVAR
            AUGMT(K,I) = UBXC(I,J)
   10    CONTINUE
   20 CONTINUE
C
C.....Compute the LU factorization
C
      CALL LUF( NVAR,AUGMT,LDA,JPAUG,IER )
      IF( IER .NE. 0 ) THEN
         IF( IER .GT. 0 ) THEN
            CALL MSGPRT( 'AUGM','Data error in LU factorization' )
            IER = -1
         ELSE
            CALL MSGPRT( 'AUGM',
     &        'Numerically singular matrix in LU factorization' )
            IER = 1
         ENDIF
         RETURN
      ENDIF
C
      IER = 0
      RETURN
C
C.....End of AUGM
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE COBAS( NVAR,NALG,W,LDW,TAU,IPV,UBAS,LDU,
     &                  WRK,SAFMIN,IER )
C
      INTEGER NVAR,NALG,LDW,IPV(*),LDU,IER
      DOUBLE PRECISION W(LDW,*),TAU(*),UBAS(LDU,*),WRK(*),SAFMIN
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     The routine computes an orthonormal basis matrix of the given
C     matrix NALG x NVAR matrix W with rank W = NALG, using the
C     LQ-factorization with pivting of W. The MDIM = NVAR - NALG
C     basis vectors are returned as the columns of the NVAR x MDIM
C     array UBAS. In addition, the arrays W, TAU, and IPV contain
C     the LQ-factorization of W as returned by LQF.
C
C     Geometrically speaking, if V denotes the NALG-dimensional
C     subspace of R^NVAR spanned by the rows of W, then UBAS
C     provides a basis of the orthogonal complement of V.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   Dimension of the ambient space
C     NALG   I   IN   Dimension of the space V
C     W      D   IN   Array of dimension LDW x NVAR, containing
C                     the given matrix
C                OUT  The LQ-factorization of the matrix W
C     LDW    I   IN   Leading  dimension of W, LDW >= NALG
C     TAU    D   OUT  Array of dimension NALG for the factors
C                     of the Householder transformations
C     IPV    I   OUT  Integer array of dimension NALG for the
C                     pivot sequence
C     UBAS   D   OUT  Array of dimension NVAR x MDIM for the basis
C                     of ker W
C     LDU    I   IN   Leading dimension of UBAS, LDU >= NVAR
C     WRK    D   WK   Work array of dimension NALG
C     SAFMIN D   IN   Safe minimum such that
C                     1.0D0/SAFMIN does not overflow
C     IER    I   OUT  Error indicator
C                     IER = 0 : no error
C                     IER = -1: input error reported by a subroutine
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL LQF,MSGPRT
C
C.....Local variables
C
      INTEGER I,J,MDIM,NVP1
C
C.....Parameters
C
      DOUBLE PRECISION ZER, ONE
      CHARACTER*(*)LNAME
      PARAMETER( ZER=0.0D0, ONE=1.0D0, LNAME='COBAS' )
C
C.......................Executable statements..........................
C
C.....Use LQF to compute the QR-factorization of W
C
      CALL LQF( NALG,NVAR,W,LDW,IPV,TAU,WRK,SAFMIN,IER )
      IF( IER .NE. 0 )THEN
         CALL MSGPRT(LNAME,'Error in the QR-factorization ')
         IER = -1
         RETURN
      ENDIF
C
C.....Compute the basis of ker(W) and store in UBAS
C
      MDIM = NVAR - NALG
      DO 20 J = 1,MDIM
         DO 10 I = 1, NVAR
            UBAS(I,J) = ZER
   10    CONTINUE
         UBAS(NALG+J,J) = ONE
   20 CONTINUE
C
      NVP1 = NVAR + 1
      DO 30 J = NALG,1,-1
         CALL HOUSL( NVP1-J,MDIM,W(J,J),LDW,TAU(J),UBAS(J,1),LDU,IER )
         IF( IER .NE. 0) THEN
            CALL MSGPRT(LNAME,
     &           'Error in the Householder transformation')
            IER = -1
            RETURN
         ENDIF
   30 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of COBAS
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE CURVT( NVAR,XC,XL,XR,XCTAN,CURV,PNORML,
     &                  PNMAX,EPMACH,IER )
C
      INTEGER NVAR,IER
      DOUBLE PRECISION CURV,EPMACH,PNMAX
      DOUBLE PRECISION XC(NVAR),XL(NVAR),XR(NVAR),XCTAN(NVAR)
      DOUBLE PRECISION PNORML(NVAR)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Let XL, XC, XR be a sequence of three points located
C     closely together along a path on the manifold and XCTAN the
C     tangent vector of the path at the middle point XC. In
C     general, these points will define a circle in a
C     two-dimensional subspace of the ambient space R^NVAR. Let
C
C         a = ||XL - XC||, b = ||XR - XC||, c = ||XL - XR||
C         s = 0.5*(a + b + c)
C
C     where the Euclidean norm is used. Then, Heron's formula
C
C         CURV = (4.0/a*b*c)*sqrt{ s*(s-a)(s-b)*(s-c) }
C
C     is used to approximate the curvature of the circle through
C     the points. The principal normal vector of the curve is
C     approximated by the vector orthogonal to XCTAN in the plane
C     spanned by the three points. Then the vector
C
C           CURV * PNORML/||PNORML|| , ||.|| Euclidean norm
C
C     represents an approximation of the diagonal component
C     V_xc(XCTAN,XCTAN) of the second fundamental tensor of M at XC.
C
C     Call the routine with a sequence of three points XL, XC, XR and
C     the tangent vector XCTAN. These data are best constructed as
C     follows. Let XC be the center point of a tangential local
C     coordinate with the orthonormal basis matrix UBXC. Choose a
C     tangent vector XCTAN = UBXC*YC of Euclidean norm one and
C     compute XL and XR by applying TPHI with H*YC and -H*YC,
C     respectively, with a suitably small step H > 0.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I  IN   dimension of the ambient space
C     MDIM   I  IN   dimension of the manifold
C     XC     D  IN   array of dimension NVAR, the center point of
C                    the local basis
C     XL     D  IN   array of dimension NVAR, the "left" neighbor
C     XR     D  IN   array of dimension NVAR, the "right" neighbor
C     XCTAN  D  IN   array of dimension NVAR, the tangent vector
C                    of the path at XC
C               OUT  the same vector normalized to Euclidean
C                    norm one.
C     CURV   D  OUT  the computed curvature
C     PNORML D  OUT  array of dimension NVAR, the computed normal
C                    vector of the path at XC
C     PNMAX  D  OUT  the maximum norm of PNORML
C     EPMACH D  IN   Relative machine accuracy
C     IER    I  OUT  error indicator
C                    IER = 0   computation successful
C                    IER = 1   the normal vector is zero.
C                              CURV + PNMAC = 0.0 are returned
C                    IER = -1  one or several of the sides of
C                              the triangle generated by the
C                              three points are zero.
C                    IER = -2  the tangent vector is zero
C                    IER = -3  the triangle is degenerate
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL MSGPRT
C
C.....Functions called
C
      DOUBLE PRECISION ABS,DDIST2,SQRT
C
C.....Parameters
C
      DOUBLE PRECISION ONE,ZER,DMAX
      CHARACTER*(*)LNAME
      PARAMETER( ONE=1.0D0, ZER=0.0D0, DMAX=0.99D0 ,LNAME='UCURV')
C
C.....Local variables
C
      INTEGER I
      DOUBLE PRECISION A,B,C,AS,BS,ALPH,DIFF,FACT1,FACT2,QUOT
      DOUBLE PRECISION SCAL,SUM,SUMINV,TFACT,TMP,XCTNRM
C
C.......................Executable statements..........................
C
C.....Compute the lengths of the sides of the triangle
C
      A = DDIST2(NVAR,XL,1,XC,1,2)
      IF (A .EQ. ZER) THEN
         CALL MSGPRT(LNAME,' The side XL--XC of'//
     &                ' the triangle has zero length ')
         IER = -1
         RETURN
      ENDIF
      B = DDIST2(NVAR,XR,1,XC,1,2)
      IF (B .EQ. ZER) THEN
         CALL MSGPRT(LNAME,' The side XR--XC of'//
     &                ' the triangle has zero length ')
         IER = -1
         RETURN
      ENDIF
      C = DDIST2(NVAR,XL,1,XR,1,2)
      IF (C .EQ. ZER) THEN
         CALL MSGPRT(LNAME,' The hypothenuse XL--XR of'//
     &                ' the triangle has zero length')
         IER = -1
         RETURN
      ENDIF
C
C.....Normalize the tangent vector
C
      XCTNRM = DDIST2(NVAR,XCTAN,1,XCTAN,1,1)
      IF (XCTNRM .EQ. ZER) THEN
         CALL MSGPRT(LNAME,' The tangent vector'//
     &                ' XCTAN has zero length ')
         IER = -2
         RETURN
      ENDIF
      TMP = ONE/XCTNRM
      DO 10 I = 1,NVAR
         XCTAN(I) = TMP*XCTAN(I)
   10 CONTINUE
      TFACT = XCTNRM*XCTNRM
C
C.....Calculate the raw principal normal and its inner product
C..... with the unit tangent vector
C
      SUM = ZER
      DO 20 I = 1,NVAR
         TMP = (XL(I)-XC(I))/A + (XR(I)-XC(I))/B
         PNORML(I) = TMP
         SUM = SUM + TMP*XCTAN(I)
   20 CONTINUE
C
C.....Compute the approximate normal vector by orthogonalizing
C.....XCTAN and the raw principal normal
C
      DO 30 I = 1,NVAR
         PNORML(I) = PNORML(I) - SUM*XCTAN(I)
   30 CONTINUE
      TMP = DDIST2(NVAR,PNORML,1,PNORML,1,1)
      IF (TMP .EQ. ZER) THEN
         CALL MSGPRT(LNAME,' Zero normal vector'//
     &               ' PNORML and CURV are set to zero')
         CURV = ZER
         PNMAX = 0
         IER = 1
      ENDIF
      TFACT = TFACT/TMP
C
C.....Scale the sides of the triangle
C
      IF (A .GT. B) THEN
         TMP = A
         A   = B
         B   = TMP
      ENDIF
      TMP = C
      IF (TMP .LT. A) THEN
         C = B
         B = A
         A = TMP
      ELSEIF (TMP .LT. B) THEN
         C = B
         B = TMP
      ENDIF
C
      SCAL = ONE/C
      AS = SCAL*A
      BS = SCAL*B
C
      SUM = AS + BS
      TMP = AS*BS
      IF (ABS(SUM).LT.ONE .OR. TMP.EQ.ZER) THEN
         CALL MSGPRT(LNAME,' The triangle is degenerate')
         IER = -3
         RETURN
      ENDIF
      QUOT = SUM/TMP
C
C.....Compute the curvature
C
      DIFF = BS - AS
      IF (DIFF .GE. SQRT(EPMACH)) THEN
         TMP = ONE - DIFF*DIFF
         FACT1 = QUOT*SQRT(TMP)
      ELSE
         FACT1 = QUOT
      ENDIF
C
      SUMINV = ONE/SUM
      TMP = ONE - SUMINV
      IF (TMP .LT. EPMACH) THEN
         CALL MSGPRT(LNAME,' The triangle is degenerate')
         IER = -3
         RETURN
      ENDIF
      IF (SUMINV .LT. DMAX) THEN
         TMP = ONE - SUMINV*SUMINV
         FACT2 = SQRT(TMP)
      ELSE
         ALPH = ACOS(SUMINV)
         FACT2 = SIN(ALPH)
      ENDIF
      CURV = FACT1*FACT2*SCAL
C
C.....Scale the computed principal normal with CURV*TFACT
C
      SCAL = CURV*TFACT
      PNMAX = ZER
      DO 40 I = 1,NVAR
         TMP = SCAL*PNORML(I)
         PNORML(I) = TMP
         IF (ABS(TMP) .GT. PNMAX) PNMAX = ABS(TMP)
   40 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of CURVT
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE DGPHI( TASK,NVAR,MDIM,DPHI,LDP,AUGMT,LDA,
     &                  UBXC,LDU,JPAUG,IER )
C
      CHARACTER*6 TASK
      INTEGER NVAR,MDIM,LDA,LDP,LDU,JPAUG(*),IER
      DOUBLE PRECISION DPHI(LDP,*),AUGMT(LDA,*),UBXC(LDU,*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     A local coordinate system at a point XC of M  induces a
C     parametrization mapping PHI from the space R^MDIM of local
C     coordinates onto M. For a given point X = PHI(Y) near XC the
C     routine computes the NVAR x MDIM derivative DPHI(Y) of PHI at Y.
C
C     Call the routine in either one of the following two different
C     modes under control of the variable TASK:
C
C     TASK = 'FACTOR'
C       Call the routine with the Jacobian DF(X) at the point X
C       stored in the first NALG = NVAR - MDIM rows of the array
C       AUGMT and the orthonormal basis of the local coordinate
C       space at XC stored in the array UBXC. Then AUGM will first
C       be called to form and factor the augmented matrix
C
C                        ( DF(X) )
C       (1)      AUGMT = (       )
C                        ( UBXC' )
C
C     TASK = 'EVAL'
C       Call the routine with the LU-factored augmented matrix (1)
C       in the array AUGMT. In this case the array UBXC is not
C       referenced.
C
C     If X is sufficiently close to XC then the augmented matrix
C     is nonsingular and the desired derivative DPHI is the
C     unique solution of the matrix equation
C
C                 ( DF(X) )         (0)
C                 (       ) DPHI  = ( )
C                 ( UBXC' )         (I)
C
C
C     TASK = 'FACTOR' is the default
C
C  Variables in the calling sequence
C  ---------------------------------
C     TASK   C   IN   Task identifier with the values
C                     TASK = 'FACTOR' (see above)
C                     TASK = 'EVAL'   (see above)
C     NVAR   I   IN   Dimension of the ambient space
C     MDIM   I   IN   Dimension of the manifold
C     DPHI   D   OUT  Array of dimension NVAR x MDIM, the computed
C                     derivative of PHI
C     LDP    I   IN   Leading dimension of DPHI, LDP >= NVAR
C     AUGMT  D        Array of dimension LDA x NVAR
C                IN   For TASK = 'FACTOR' contains in
C                     its NALG = NVAR - MDIM rows the Jacobian DF(X)
C                     For TASK = 'EVAL' the LU factored augmented
C                     matrix (1)
C                OUT  the LU factored augmented matrix (1)
C     LDA    I   IN   Leading dimension of AUGMT, LDB >= NVAR
C     UBXC   D        Array of dimension LDU x MDIM
C                IN   For TASK = 'FACTOR' contains the local basis
C                     matrix at XC
C                     For TASK = 'EVAL' not used
C     LDU    I   IN   Leading dimension of UBXC, LDU >= NVAR
C     JPAUG  D        Array of dimension NVAR
C                IN   For TASK = 'EVAL' the pivot array of the
C                     LU factorization of AUGMT
C                OUT  The pivot array of the LU factorization of AUGMT
C     IER    I   OUT  Error indicator
C                     IER =  0 no error
C                     IER = -1 error return from a subroutine
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL LUF,LUSK,MSGPRT
C
C.....Local variables
C
      INTEGER I,J,NALG
C
C.....Parameters
C
      DOUBLE PRECISION ZER,ONE
      CHARACTER*(*)LNAME
      PARAMETER( ZER = 0.0D0, ONE = 1.0D0, LNAME='DGPHI' )
C
C.......................Executable statements..........................
C
C.....Check TASK
C
      IF( TASK .NE. 'EVAL' )THEN
C
C........Set up and factor the augmented matrix
C
         CALL AUGM(NVAR,MDIM,AUGMT,LDA,UBXC,LDU,JPAUG,IER)
         IF (IER .NE. 0) THEN
            CALL MSGPRT(LNAME,
     &        'Error in forming or factoring the augmented matrix')
            IER = -1
            RETURN
         ENDIF
      ENDIF
C
C.....Solve for the columns of DPHI
C
      NALG = NVAR - MDIM
      DO 20 J = 1,MDIM
         DO 10 I = 1,NVAR
            DPHI(I,J) = ZER
  10     CONTINUE
         DPHI(NALG+J,J) = ONE
  20  CONTINUE
      CALL LUSK(NVAR,AUGMT,LDA,JPAUG,DPHI,NVAR,MDIM,IER)
      IF (IER .NE. 0) THEN
         CALL MSGPRT(LNAME,
     &              'Error in solving the augmented equation')
         IER = -1
         RETURN
      ENDIF
C
      IER = 0
      RETURN
C
C.....End of DGPHI
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE D2GPHI( TASK,NVAR,MDIM,D2PHI,DPNRM,AUGMT,LDA,
     &                   UBXC,LDU,JPAUG,IER )
C
      CHARACTER*6 TASK
      INTEGER NVAR,MDIM,LDA,LDU,JPAUG(*),IER
      DOUBLE PRECISION D2PHI(*),DPNRM,AUGMT(LDA,*),UBXC(LDU,*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     A local coordinate system at a point XC of M  induces a
C     parametrization mapping PHI from the space R^MDIM of local
C     coordinates onto M. For a given point X = PHI(Y) near XC and
C     a local direction vector YP the routine computes the second
C     derivative term D2PHI(Y)(YP,YP) of PHI. The local direction
C     is not explicitly used. Instead, upon input, the array
C     D2PHI is assumed to contain in its first NALG = NVAR - MDIM
C     components the second derivative term D2F(X)(DX,DX) formed
C     with the global direction DX = DPHI(Y)*YP. In addition, the
C     routine returns the maximum norm DPNRM of D2PHI(Y)(YP,YP).
C
C     Call the routine in either one of the following two different
C     modes under control of the variable TASK:
C
C     TASK = 'FACTOR'
C       Call the routine with the Jacobian DF(X) at the point X
C       stored in the first NALG = NVAR - MDIM rows of the array
C       AUGMT and the orthonormal basis of the local coordinate
C       space at XC stored in the array UBXC. Then AUGM will first
C       be called to set up and factor the augmented matrix
C
C                        ( DF(X) )
C       (1)      AUGMT = (       )
C                        ( UBXC' )
C
C     TASK = 'EVAL'
C       Call the routine with the LU-factored augmented matrix (1)
C       in the array AUGMT. In this case the array UBXC is not
C       referenced.
C
C     If X is sufficiently close to XC then the augmented matrix
C     is nonsingular and the desired term D2PHI is the
C     unique solution of the matrix equation.
C
C                ( DF(X) )          ( D2F(X)(DX,DX) )
C                (       ) D2PHI  = (               )
C                ( UBXC' )          (      0        )
C
C     TASK = 'FACTOR' is the default
C
C  Variables in the calling sequence
C  ---------------------------------
C     TASK   C   IN   Task identifier with the values
C                     TASK = 'FACTOR' (see above)
C                     TASK = 'EVAL'   (see above)
C     NVAR   I   IN   Dimension of the ambient space
C     MDIM   I   IN   Dimension of the manifold
C     D2PHI  D        Array of dimension NVAR
C                IN   Contains in its first NALG components the
C                     second derivative term D2F(X)(DX,DX)
C               OUT   The computed derivative term D2PHI(Y)(YP,YP)
C     AUGMT  D        Array of dimension LDA x NVAR
C                IN   For TASK = 'FACTOR' contains in
C                     its NALG = NVAR - MDIM rows the Jacobian DF(X)
C                     For TASK = 'EVAL' the LU factored augmented
C                     matrix (1)
C                OUT  the LU factored augmented matrix (1)
C     LDA    I   IN   Leading dimension of AUGMT, LDB >= NVAR
C     UBXC   D        Array of dimension LDU x MDIM
C                IN   For TASK = 'FACTOR' contains the local basis
C                     matrix at XC
C                     For TASK = 'EVAL' not used
C     LDU    I   IN   Leading dimension of UBXC, LDU >= NVAR
C     JPAUG  D        Array of dimension NVAR
C                IN   For TASK = 'EVAL' the pivot array of the
C                     LU factorization of AUGMT
C                OUT  The pivot array of the LU factorization of AUGMT
C     IER    I   OUT  Error indicator
C                     IER =  0 no error
C                     IER = -1 error return from a subroutine
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL LUF,LUSK,MSGPRT
C
C.....Function called
C
      DOUBLE PRECISION ABS
C
C.....Local variables
C
      INTEGER I,J,NALG
C
C.....Parameters
C
      DOUBLE PRECISION ZER
      CHARACTER*(*)LNAME
      PARAMETER( ZER = 0.0D0, LNAME='D2GPHI' )
C
C.......................Executable statements..........................
C
C.....Check TASK
C
      IF( TASK .NE. 'EVAL' )THEN
C
C........Set up and factor the augmented matrix
C
         CALL AUGM(NVAR,MDIM,AUGMT,LDA,UBXC,LDU,JPAUG,IER)
         IF (IER .NE. 0) THEN
            CALL MSGPRT(LNAME,
     &        'Error in forming or factoring the augmented matrix')
            IER = -1
            RETURN
         ENDIF
      ENDIF
C
C.....Solve for the columns of DPHI
C
      NALG = NVAR - MDIM
      DO 10 I = 1,NALG
         D2PHI(I) = -D2PHI(I)
  10  CONTINUE
      DO 20 J = 1,MDIM
         D2PHI(NALG+J) = ZER
  20  CONTINUE
      CALL LUS1( NVAR,AUGMT,LDA,JPAUG,D2PHI,IER )
      IF (IER .NE. 0) THEN
         CALL MSGPRT(LNAME,
     &             'Error in solving the augmented equation')
         IER = -1
         RETURN
      ENDIF
C
C.....Determine the maximum norm
C
      DPNRM = ZER
      DO 30 I = 1, NVAR
         IF(DPNRM .LT. ABS(D2PHI(I))) DPNRM = ABS(D2PHI(I))
   30 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of D2GPHI
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE DTPHI( NVAR,MDIM,DPHI,LDP,UBXN,LDUN,UBXC,LDUC,
     &                  WRKMAT,LDW,IWRK,WRK,IER )
C
      INTEGER IER,NVAR,MDIM,LDP,LDUN,LDUC,LDW,IWRK(*)
      DOUBLE PRECISION DPHI(LDP,*),UBXN(LDUN,*),UBXC(LDUC,*)
      DOUBLE PRECISION WRKMAT(LDW,*),WRK(*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     A tangential local coordinate system centered at a point
C     XC on the manifold induces a parametrization mapping PHI from
C     the space R^MDIM of local coordinates onto M. For a point
C     XN = PHI(YN) near XC the routine computes the derivative
C     DPHI(YN) of PHI at YN.
C
C     Call the routine with the array UBXC containing an
C     orthonormal tangential basis at XC, and the array UBXN
C     containing an orthonormal tangential basis at XN.
C     The bases may or may not have been computed by TBAS.
C
C     The routine uses the fact that
C
C           DPHI = UBXN * [UBXC' * UBXN]^{-1} ;
C
C     and hence that the solutions of the MDIM x MDIM linear
C     systems [UBXC' * UBXN] *Z = UBXN' are the rows of DPHI.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   Dimension of X
C     MDIM   I   IN   Dimension of the manifold
C     DPHI   D   OUT  Array of dimension NVAR x MDIM, containing
C                     the computed derivative
C     LDP    I   IN   Leading dimension of DPHI, LDP .GE. NVAR
C     UBXN   D   IN   Array of dimension LDUM x MDIM containing
C                     the tangent basis at XN
C     LDUN   I   IN   Leading dimension of UBXN, LDUN >= NVAR
C     UBXC   D   IN   Array of dimension LDUC x MDIM containing
C                     the tangent basis at XC
C     LDUC   I   IN   Leading dimension of UBXC, LDUC .GE. NVAR
C     WRKMAT D   WK   Work array of dimension LDW x MDIM
C     LDW    I   IN   Leading dimension of WRKMAT, LDW >= MDIM.
C     IWRK   I   WK   Integer work array of dimension MDIM
C     WRK    D   WK   Work array of dimension NVAR
C     IER    I   OUT  error indicator
C                     IER =  0 no error
C                     IER = -1 error return from a subroutine
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL LUF,LUS1,MSGPRT
C
C.....Parameter
C
      DOUBLE PRECISION ZER
      CHARACTER*(*)LNAME
      PARAMETER( ZER=0.0D0, LNAME='DTPHI' )
C
C.....Local variables
C
      INTEGER I,J,K
      DOUBLE PRECISION SUM
C
C.......................Executable statements..........................
C
C.....Form the product WRKMAT = UBXN' * UBXC
C
      DO 30 J = 1,MDIM
         DO 20 I = 1,MDIM
            SUM = ZER
            DO 10 K = 1,NVAR
               SUM = SUM + UBXN(K,I)*UBXC(K,J)
   10       CONTINUE
            WRKMAT(I,J) = SUM
   20    CONTINUE
   30 CONTINUE
C
C.....Compute the LU-factorization of WRKMAT
C
      CALL LUF(MDIM,WRKMAT,LDW,IWRK,IER)
      IF(IER .NE. 0) THEN
         CALL MSGPRT(LNAME,'Error in LU factorization')
         IER = -1
         RETURN
      ENDIF
C
C.....Solve WRKMAT * Z = UBXN'
C
      DO 60 K = 1,NVAR
         DO 40 J = 1,MDIM
            WRK(J) = UBXN(K,J)
   40    CONTINUE
         CALL LUS1(MDIM,WRKMAT,MDIM,IWRK,WRK,IER)
         IF(IER .NE. 0) THEN
            CALL MSGPRT(LNAME,'Error in solving the linear system')
            IER = -1
            RETURN
         ENDIF
         DO 50 J = 1,MDIM
            DPHI(K,J) = WRK(J)
   50    CONTINUE
   60 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of DTPHI
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE GCBAS( NVAR,MDIM,UBX,LDU,DFX,LDF,TAUX,JPX,
     &                  WRK,IWRK,SAFMIN,IER )
C
      INTEGER NVAR,MDIM,LDU,LDF,JPX(*),IER,IWRK(*)
      DOUBLE PRECISION UBX(LDU,*),DFX(LDF,*),TAUX(*),WRK(*)
      DOUBLE PRECISION SAFMIN
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Call the routine with the Jacobian DF(X) stored in the first
C     NALG = NVAR - MDIM rows of the array AUGMT. The routine
C     returns in the array UBX an orthonormal basis of a local
C     parametrization consistiing of MDIM canonical basis vectors
C     of the ambient space R^NVAR. (The identifier C in the name
C     of the routine refers to the word "canonical").
C
C     The local coordinate space is defined as the orthogonal
C     complement of a space W obtained from the Jacobian DF(X).
C     More specifically, let
C
C            DF(XC)*P = Q * (R, S)
C
C     be the QR-factorization (with column pivoting) of the Jacobian,
C     where P denotes some NALG (=NVAR-MDIM) dimensional permutation
C     matrix. Then P selects NALG column-vectors of the Jacobian and
C     the NALG x NALG matrix formed of these columns is nonsingular.
C     The selected columns correspond to NALG canonical basis vectors
C     of R^NVAR which span the desired space W. The basis matrix
C     UBXC then consists of the remaining canonical basis vectors
C     not in W.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   Dimension of the ambient space
C     MDIM   I   IN   Dimension of the manifold
C     UBX    D   OUT  Array of dimension LDU x MDIM, the computed
C                     basis of the local coordinate system
C     LDU    I   IN   Row dimension of UBX, LDU >= NVAR
C     DFX    D   IN   Array of dimension LDF x NVAR, containing
C                     the Jacobian DF(XC)
C                OUT  The LQ-factorization of the Jacobian
C     LDB    I   IN   Leading dimension of DFXC, LDB >= NALG
C     TAUX   D   OUT  Array of dimension NALG, the factors of
C                     Householder reflectors
C     JPX    I   OUT  Integer array of dimension NALG for the
C                     pivot sequence
C     IWRK   I   WRK  Work array of dimension NVAR
C     WRK    D   WRK  Work array of dimension NVAR
C     SAFMIN D   IN   Machine constant, safe minimum such that
C                      1.0D0/SAFMIN does not overflow
C     IER    I   IN   error indicator
C                     IER =  0 no error
C                     IER = -1 error return from QRF
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutine called
C
      EXTERNAL QRF,MSGPRT
C
C.....Parameters
C
      DOUBLE PRECISION ZER, ONE
      PARAMETER( ZER=0.0D0, ONE=1.0D0 )
C
C.....Local variables
C
      INTEGER I,J,L,NALG
C
C.......................Executable statements..........................
C
C.....Use QRF to compute the QR-factorization of DFXC with pivoting
C
      NALG = NVAR - MDIM
      CALL QRF(NALG,NVAR,DFX,LDF,JPX,TAUX,WRK,SAFMIN,IER)
      IF(IER .NE. 0)THEN
         CALL MSGPRT('GCBAS',
     &      ' Error in the LQ factorization of the Jacobian')
         IER = -1
         RETURN
      ENDIF
C
C.....Zero the work vetor
C
      DO 10 I = 1, NVAR
         IWRK(I) = 0
   10 CONTINUE
C
C.....Collect the indices of the used canonical basis vectors
C.....of R^NVAR in IWRK and zero the matrix UBX
C
      DO 30 J = 1,MDIM
         L = JPX(J)
         IWRK(L) = 1
         DO 20 I = 1, NVAR
           UBX(I,J) = ZER
   20    CONTINUE
   30 CONTINUE
C
C.....Store the remaning canonical basis vectors in UBX
C
      L = 1
      DO 40 I = 1, NVAR
        IF(IWRK(I) .EQ. 0) THEN
           UBX(I,L) = ONE
           L = L + 1
        ENDIF
   40 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of GCBAS
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE GLOB( NVAR,MDIM,Y,X,XC,UBXC,LDU )
C
      INTEGER NVAR,MDIM,LDU
      DOUBLE PRECISION Y(*),X(NVAR),XC(*),UBXC(LDU,*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Let the NVAR x MDIM array UBXC contain an orthonormal basis
C     of a local parametrization at the point XC. For a given
C     point with local coordinates Y in this parametrization, the
C     routine computes the corresponding point
C
C                 X = XC + UBXC*Y
C
C     in global coordinates belonging to the local coordinate
C     space V = XC + span(UBXC) in R^NVAR.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   dimension of the ambient space
C     MDIM   I   IN   dimension of the space V
C     Y      D   IN   array of dimension MDIM, the given point
C                     of V in local coordinates
C     X      D   OUT  array of dimension NVAR, the computed global
C                     point corresponding to Y
C     XC     D   IN   array of dimension NVAR, the origin of the
C                     basis of V.
C     UBXC   D   IN   array of dimension  LDU x MDIM, containing
C                     the orthonormal local basis at XC
C     LDU    I   IN   leading dimension of UBXC, LDU >= NVAR
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Local variables
C
      INTEGER I,J
      DOUBLE PRECISION SUM
C
C.......................Executable statements..........................
C
      DO 20 I = 1,NVAR
         SUM = XC(I)
         DO 10 J = 1,MDIM
            SUM = SUM + UBXC(I,J)*Y(J)
 10      CONTINUE
         X(I) = SUM
 20   CONTINUE
C
      RETURN
C
C.....End of GLOB
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE GNBAS( NVAR,MDIM,UBX,LDU,DFX,LDF,
     &                  WRK1,WRK2,IWRK,SAFMIN,IER )
C
      INTEGER NVAR,MDIM,LDU,LDF,IER,IWRK(*)
      DOUBLE PRECISION DFX(LDF,*),UBX(LDU,*),WRK1(*),WRK2(*)
      DOUBLE PRECISION SAFMIN
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     The routine establishes an orthonormal basis UBX of a general
C     local coordinate system at a point X of the manifold which
C     contains the NVAR-th canonical basis vector of R^NVAR.
C
C     The local coordinate space is defined as the orthogonal
C     complement of a space W obtained by applying COBAS to the
C     NALG x (NVAR-1) matrix consisting of the first NVAR-1
C     columns of DF(X). Then the NVAR-th canonical basis vector
C     is added to the basis.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR    I   IN   Dimension of the ambient space
C     MDIM    I   IN   Dimension of the manifold
C     UBX     D   OUT  Array of dimension LDU x MDIM, the computed
C                      basis of the local coordinate system
C     LDU     I   IN   Leading dimension of UBX, LDU >= NVAR
C     DFX     D   IN   Array of dimension LDF x NVAR containing
C                      the Jacobian DF(X)
C                 OUT  the Jacobian is destroyed
C     LDF     I   IN   Leading dimension of DFX, LDF >= NALG
C     WRK1    D   WK   Work array of dimension NVAR
C     WRK2    D   WK   Work array of dimension NVAR
C     IWRK    I   WK   Work array of dimension NVAR
C     SAFMIN  D   IN   Machine constant, safe minimum such that
C                      1.0D0/SAFMIN does not overflow
C     IER     I   IN   error indicator
C                      IER =  0 no error
C                      IER = -1 error return from COBAS
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL COBAS,MSGPRT
C
C.....Parameter
C
      DOUBLE PRECISION ZER, ONE
      PARAMETER( ZER = 0.0D0, ONE=1.0D0 )
C
C.....Local variables
C
      INTEGER I,J,NALG,NRED
C
C.......................Executable statements..........................
C
      NALG = NVAR - MDIM
      NRED = NVAR - 1
C
C.....Compute the reduced basis matrix
C
      CALL COBAS( NRED,NALG,DFX,LDF,WRK1,IWRK,UBX,LDU,
     &            WRK2,SAFMIN,IER )
      IF(IER .NE. 0)THEN
         CALL MSGPRT('GNBAS',' Error in computing the basis')
         IER = -1
         RETURN
      ENDIF
C
C.....Augment the basis by the NVAR-th canonical basis vector
C
      DO 10 J = 1,NVAR
         UBX(J,MDIM) = ZER
   10 CONTINUE
      DO 20 I = 1, MDIM
         UBX(NVAR,I) = ZER
   20 CONTINUE
      UBX(NVAR,MDIM) = ONE
C
      IER = 0
      RETURN
C
C.....End of GNBAS
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE GPHI( TASK,NVAR,MDIM,Y,X,FX,XC,UBXC,LDU,
     &                 AUGMT,LDA,JPAUG,EPMACH,ISTEP )
C
      CHARACTER*6 TASK
      INTEGER NVAR,MDIM,LDA,LDU,JPAUG(*),ISTEP
      DOUBLE PRECISION Y(*),X(*),FX(*),XC(*),UBXC(LDU,*)
      DOUBLE PRECISION AUGMT(LDA,*),EPMACH
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Given a general coordinate system of M at XC specified by the
C     arrays [XC, UBXC, AUGMT, JPAUG], the routine computes the
C     point X on M in global coordinates which has the given vector
C     Y as its local representation in the coordinate system.
C
C     The routine uses reverse communication to access the user
C     subroutine for the function F.  This is done under the control
C     of the character variable TASK:
C
C     For the initial call set TASK = 'START'.
C     The routine returns repeatedly with TASK = 'EVAL' signifying
C     a request for the evaluation of F at the global point X.
C     Return to the routine with the computed vector FX = F(X)
C     but no changes of TASK or any other variable.
C
C     By the conventions for general coordinate systems, the array
C     AUGMT contains the LU factorization of the augmented matrix
C
C                    ( DF(XC) )
C            AUGMT = (        )
C                    ( UBXC^T )
C
C     For the given MDIM dimensional local vector Y the routine
C     uses the chord Newton method specified by this augmented
C     matrix to compute the (unique) point X on M with local
C     coordinate Y; that is the solution of the nonlinear system
C
C           ( F(X)           )       (0)
C           (                )   =   ( )
C           ( UBXC^T(X - XC) )       (Y)
C
C     The process converges if ||Y|| is sufficiently small; that is,
C     if X is sufficiently close to the center point XC of the
C     local coordinate system.
C
C     At the final return, TASK may have either one of the following
C     values:
C
C     1. TASK = 'DONE'    process converged
C     2. TASK = 'DIVERG'  divergence detected
C     3. TASK = 'STPCNT'  maximal number of steps NMAX exceeded
C     4. TASK = 'SINGUL'  Singular iteration matrix
C     5. TASK = 'ERROR'   Wrong input
C
C     The failure-exits 2. and 3. signify that ||Y|| was not small
C     enough for the given local coordinate system. This may not be
C     a fatal error but may only require the choice of a smaller Y.
C     On the other hand the failure exits 4, and 5. are fatal.
C
C  Variables in the calling sequence
C  ---------------------------------
C     TASK   C  IN   Task identifier
C                    TASK = 'START'  start the process
C                    TASK = 'EVAL'   the function vector F(X) is
C                                     available in FX
C               OUT  TASK = 'EVAL'   Compute FX = F(X) for given X
C                    TASK = 'DONE'   process converged
C                    TASK = 'DIVERG'  divergence detected
C                    TASK = 'STPCNT' Maximal number of steps exceeded
C                    TASK = 'SINGUL' Singular iteration matrix
C                    TASK = 'ERROR'  Wrong input
C     NVAR   I  IN   Dimension of the ambient space
C     MDIM   I  IN   Dimension of the manifold
C     Y      D  IN   Array of dimension MDIM, the given local point
C     X      D  I/O  Array of dimension NVAR, global point on M
C     FX     D  IN   Array of dimension NVAR which contains in its
C                    first NALG components the computed FX = F(X)
C                    and which will be augmented by the routine
C     XC     D  IN   The center point on M of the local basis
C     UBXC   D  IN   Array of dimension LDU x MDIM, the basis of
C                    the local coordinate system
C     LDU    I  IN   Leading dimension of UBXC, LDU .GE. NVAR
C     AUGMT  D  IN   Array of dimension LDA x NVAR, the LU-factored
C                    augmented matrix
C     LDA    I  IN   Leading dimension of AUGMT, LDA .GE. NVAR
C     JPAUG  I  OUT  Array of dimension NVAR, the pivot array of
C                    the LU-factorization
C     EPMACH D  IN   Machine epsilon
C     ISTEP  I  OUT  Number of iteration steps
C                    ISTEP = -2  Input error, no X returned
C                    ISTEP = -1  Y is small, X = XC is returned
C                    ISTEP = 0   the residual is small, X = UGLOB(Y)
C                                is returned
C                    ISTEP > 0   number of iteration steps taken. When
C                                an error condition occurred then
C                                ISTEP has the value of the current
C                                step
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....External subroutines
C
      EXTERNAL GLOB,LUS1
C
C.....Function calls
C
      DOUBLE PRECISION DDIST2,SQRT
C
C.....Parameters
C
      INTEGER MAXSTP
      DOUBLE PRECISION FAC, FACLOW
      PARAMETER( FAC=2.0D0, FACLOW=1.05D0, MAXSTP=25 )
      DOUBLE PRECISION ZER, SIXTN
      PARAMETER( ZER=0.0D0, SIXTN=1.6D1 )
C
C.....Local variables
C
      INTEGER I,IER
C
C.....Variables saved between calls
C
      INTEGER NALG
      DOUBLE PRECISION FNRM,FNRML,ABSERR,RELERR,RELSIX
      DOUBLE PRECISION STEP,STEPL,TOLX,TOLSIX,XNRM
      SAVE NALG,FNRM,FNRML,ABSERR,RELERR,RELSIX,
     &     STEP,STEPL,TOLX,TOLSIX,XNRM
C
      DATA FNRM/0.0D0/, STEPL/0.0D0/
C
C.......................Executable statements..........................
C
C.....Check TASK
C
      IF (TASK .EQ. 'EVAL') THEN
         GOTO 100
C
      ELSEIF (TASK .EQ. 'START') THEN
C
C........If ||Y|| < EPMACH set X = XC and return
C
         STEP = DDIST2(MDIM,Y,1,Y,1,1)
         IF( STEP .LT. EPMACH ) THEN
            DO 10 I = 1,NVAR
               X(I) = XC(I)
   10       CONTINUE
            ISTEP = -1
            GOTO 200
         ENDIF
C
C........Initialize data
C
         NALG   = NVAR - MDIM
         XNRM   = DDIST2(NVAR,XC,1,XC,1,1)
         RELERR = SQRT(EPMACH)
         ABSERR = EPMACH/SQRT(RELERR)
         TOLX   = RELERR*XNRM + ABSERR
         TOLSIX = SIXTN*TOLX
         RELSIX = SIXTN*RELERR
C
         ISTEP  = 0
C
C........Globalize the point Y
C
         CALL GLOB(NVAR,MDIM,Y,X,XC,UBXC,LDU)
C
C........Ask for evaluation of FX = F(X)
C
         TASK = 'EVAL'
         RETURN
      ELSE
C
C........TASK has unacceptable value
C
         ISTEP = -2
         GOTO 240
      ENDIF
C
C.....Iteration point
C
  100 CONTINUE
C
C.....Compute the function norm
C
      FNRML = FNRM
      FNRM  = DDIST2(NALG,FX,1,FX,1,1)
C
C.....Test for convergence
C
      IF (ISTEP .EQ. 0) THEN
         IF (FNRM .LT. ABSERR) GOTO 200
      ELSE
C
C........Strong acceptance test
C
         IF(FNRM .LE. RELERR .AND. STEP .LE. TOLX) GOTO 200
C
C........Weak acceptance tests
C
         IF(FNRM .LE. ABSERR .OR. STEP .LE. ABSERR) GOTO 200
         IF(ISTEP .GT. 1) THEN
            IF((FNRM+FNRML) .LE. RELERR
     &             .AND. STEP.LE.TOLSIX) GOTO 200
            IF(FNRM .LE. RELSIX
     &             .AND. (STEP+STEPL).LE.TOLX) GOTO 200
         ENDIF
C
C........Divergence test
C
         IF (ISTEP .EQ. 1) THEN
            IF (FNRM .GT. (FAC*FNRML+RELERR)) GOTO 210
         ELSEIF (ISTEP .EQ. 2) THEN
            IF(STEP .GT. (FAC*STEPL+RELERR)) GOTO 210
            IF (FNRM .GT. (FACLOW*FNRML+RELERR)) GOTO 210
         ELSE
            IF(STEP .GT. (FACLOW*STEPL+RELERR)) GOTO 210
            IF (FNRM .GT. (FACLOW*FNRML+RELERR)) GOTO 210
         ENDIF
      ENDIF
C
C.....Check for maximum number of steps
C
      IF (ISTEP .GT. MAXSTP) GOTO 220
C
C.....Augment the function vector
C
      DO 110 I = 1,MDIM
         FX(NALG+I) = ZER
  110 CONTINUE
C
C.....Solve the chord Newton equation
C
      CALL LUS1( NVAR,AUGMT,LDA,JPAUG,FX,IER )
      IF (IER .NE. 0) GOTO 230
C
C.....Form the next iterate
C
      STEPL = STEP
      STEP = DDIST2(NVAR,FX,1,FX,1,1)
      DO 130 I = 1,NVAR
         X(I) = X(I) - FX(I)
  130 CONTINUE
C
      ISTEP = ISTEP + 1
      TASK = 'EVAL'
      RETURN
C
C.....Common returns
C
  200 CONTINUE
      TASK = 'DONE'
      RETURN
C
  210 TASK = 'DIVERG'
      RETURN
C
  220 TASK = 'STPCNT'
      RETURN
C
  230 TASK = 'SINGUL'
      RETURN
C
  240 TASK = 'ERROR'
      RETURN
C
C.....End of GPHI
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE MOVFR( NVAR,MDIM,UBX1,LDU1,UBX2,LDU2,WRKMAT,LDW,
     &                  WRK1,WRK2,WRK3,EPMACH,SAFMIN,IER )
C
      INTEGER NVAR,MDIM,LDU1,LDU2,LDW,IER
      DOUBLE PRECISION UBX1(LDU1,*),UBX2(LDU2,*)
      DOUBLE PRECISION WRKMAT(LDW,*),WRK1(*),WRK2(*),WRK3(*)
      DOUBLE PRECISION EPMACH,SAFMIN
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Given two orthonormal tangent bases UBX1 and UBX2 at some
C     neighboring points of the manifold, the routine uses the moving
C     frame algorithm of W. Rheinboldt, Num. Math. 53, 1988, 165-181,
C     to adjust the orientation of UBX2 to match that of UBX1. The
C     matrix UBX1 is destroyed and UBX2 is overwritten with the
C     reoriented form of UBX2.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   dimension of the ambient space
C     MDIM   I   IN   dimension of the manifold
C     UBX1   D   IN   double precision array of dimension LDU1 x MDIM,
C                     the first local basis
C     LDU1   I   IN   leading dimension of UBX1, LDU1 >= NVAR
C     UBX2   D   IN   double precision array of dimension LDU2 x MDIM,
C                     the second local basis
C     LDU2   I   IN   leading dimension of UBX2, LDU2 >= NVAR
C     WRKMAT D   WK   work array of dimension MDIM x MDIM
C     LDW    I   IN   leading dimension of WRKMAT, LDW .GE. MDIM
C     WRK1   D   WK   Work array of dimension MDIM
C     WRK2   D   WK   Work array of dimension MDIM
C     WRK3   D   WK   Work array of dimension MDIM
C     EPMACH D   IN   Relative machine accuracy
C     SAFMIN D   IN   Safe minimum
C     IER    I   OUT  Error indicator
C                     IER =  0 no error
C                     IER = -1 error returns from subroutines
C                     IER =  1 SVD incomplete
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....External subroutines called
C
      EXTERNAL BIDIA,SVD
C
C.....Parameters
C
      DOUBLE PRECISION ZER
      CHARACTER*(*)LNAME
      PARAMETER( ZER=0.0D0, LNAME='MOVFR' )
C
C.....Local variables
C
      INTEGER I,J,K,JOB
      DOUBLE PRECISION SUM
C
C.......................Executable statements..........................
C
C.....Form the product matrix UBX1^T * UBX2 and store it in WRKMAT
C
      DO 30 J = 1,MDIM
         DO 20 I = 1,MDIM
            SUM = ZER
            DO 10 K = 1,NVAR
               SUM = SUM + UBX2(K,I)*UBX1(K,J)
   10       CONTINUE
            WRKMAT(I,J) = SUM
   20    CONTINUE
   30 CONTINUE
C
C.....Bidiagonalize WRKMAT
C
      JOB = 1
      CALL BIDIA( MDIM,MDIM,WRKMAT,LDW,WRK1,WRK2,UBX1,LDU1,
     &            WRK3,SAFMIN,JOB,IER )
      IF( IER .NE. 0) THEN
         CALL MSGPRT(LNAME,'Error in bidiagonalizing'//
     &               ' the reorientation matrix')
         IER = -1
         RETURN
      ENDIF
C
C.....Compute the SVD
C
      CALL SVD(MDIM,MDIM,WRK1,WRK2,UBX1,LDU1,WRKMAT,LDW,
     &         EPMACH,SAFMIN,1,IER)
      IF (IER .NE. 0) THEN
         IF( IER .GT. 0) THEN
         CALL MSGPRT(LNAME,'Error in the SVD decomposition'//
     &                       ' of the reorientation matrix')
            IER = -1
            RETURN
         ELSE
C
C...........SVD incomplete, no reorientation is done and the
C...........code returns IER = 1 and the unmodified basis UBX2
C
            CALL MSGPRT(LNAME,'SVD incomplete, no reorientation')
            IER = 1
            RETURN
         ENDIF
      ENDIF
C
C.....Rotate the new basis by multiplying it from the right by
C.....A*B^T where A and B denote the matrices of left and right
C.....singular vectors, respectively
C
      DO 80 I = 1,NVAR
         DO 50 J = 1,MDIM
            SUM = ZER
            DO 40 K = 1,MDIM
               SUM = SUM + UBX2(I,K)*UBX1(K,J)
   40       CONTINUE
            WRK1(J) = SUM
   50    CONTINUE
         DO 70 J = 1,MDIM
            SUM = ZER
            DO 60 K = 1,MDIM
               SUM = SUM + WRK1(K)*WRKMAT(J,K)
   60       CONTINUE
            UBX2(I,J) = SUM
   70    CONTINUE
   80 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of MOVFR
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE ORIENT( NVAR,MDIM,UBX1,LDU1,UBX2,LDU2,
     &                   WRKMAT,LDW,IWRK,IER )
C
      INTEGER NVAR,MDIM,LDU1,LDU2,LDW,IWRK(*),IER
      DOUBLE PRECISION UBX1(LDU1,*),UBX2(LDU2,*)
      DOUBLE PRECISION WRKMAT(LDW,*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     For two orthonormal tangent bases UBX1 and UBX2 at some
C     neighboring points of the manifold, the routine uses the
C     following simple combinatorial algorithm to adjust the
C     orientation of UBX2 to match that of UBX1:
C
C     (1) compute the matrix A = UBX1^T * UBX2
C     (2) mark all columns of A as uncovered
C     (3) In each row i of A find the index k(i) of the column which
C         is uncovered and for which the corresponding element of A
C         is largest in modulus.
C     (4) Cover column k.
C     (5) if the element (i,k) of A is negative then change the
C         sign of the k-th vector of UBX2
C     (6) If all columns are covered the map i --> k(i) defines
C         the desired permutation of the columns of UBX2.
C
C     The algorithm may fail even though a more sophisticated method
C     may be able to produce a reorientation. See e.g. the moving
C     frame algorithm implemented in the subroutine MOVFR.
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   dimension of the ambient space
C     MDIM   I   IN   dimension of the manifold
C     UBX1   D   IN   double precision array of dimension LDU1 x MDIM,
C                     the first local basis
C     LDU1   I   IN   row dimension of UBX1, LDU1 >= NVAR
C     UBX2   D   IN   double precision array of dimension LDU2 x MDIM,
C                     the second local basis
C     LDU2   I   IN   row dimension of UBX2, LDU2 >= NVAR
C     WRKMAT D   WK   work array of dimension MDIM x MDIM
C     LDW    I   IN   row dimension of WRKMAT, LDW .GE. MDIM
C     IWRK   I   WK   integer work array of dimension MDIM
C     IER    I   OUT  error indicator
C                     IER =  0 no error
C                     IER = -1 orientation failed
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....External subroutine called
C
      EXTERNAL MSGPRT
C
C.....FORTRAN function called
C
      DOUBLE PRECISION SIGN
C
C.....Parameters
C
      DOUBLE PRECISION ZER,ONE
      PARAMETER( ZER = 0.0D0, ONE=1.0D0 )
C
C.....Local variables
C
      INTEGER J,K,L,MAXL,ITMP
      DOUBLE PRECISION ESGN,TMP,EMAX
C
C.......................Executable statements..........................
C
C.....Form the product matrix A = UBX1^T * UBX2 and store it in
C.....WRKMAT. At the same time set IWRK to zero
C
      IER = 0
      DO 30 L=1,MDIM
         IWRK(L) = 0
         DO 20 J=1,MDIM
            TMP = ZER
            DO 10 K=1,NVAR
               TMP = TMP + UBX1(K,J)*UBX2(K,L)
   10       CONTINUE
            WRKMAT(J,L) = TMP
   20    CONTINUE
   30 CONTINUE
C
C.....In each row determine the largest, uncovered element
C
      DO 60 L = 1,MDIM
         EMAX = ZER
         MAXL = 0
         ESGN = ONE
         DO 40 J = 1,MDIM
            IF( IWRK(J) .EQ. 0) THEN
               TMP = WRKMAT(L,J)
               IF (ABS(TMP) .GT. EMAX) THEN
                  EMAX = ABS(TMP)
                  MAXL = J
                  ESGN = SIGN(ONE,TMP)
               ENDIF
            ENDIF
   40    CONTINUE
         IF (MAXL .NE. 0) THEN
            IWRK(MAXL) = L
         ELSE
            CALL MSGPRT('ORIENT',' Basis reorientation failed')
            IER = -1
            RETURN
         ENDIF
         IF( ESGN .LT. ZER) THEN
            DO 50 K = 1,NVAR
               UBX2(K,MAXL) = -UBX2(K,MAXL)
   50       CONTINUE
         ENDIF
   60 CONTINUE
C
C.....If MDIM = 1 then return
C
      IF (MDIM .EQ. 1) RETURN
C
C.....Reorient the basis
C
      DO 80 J = 1,MDIM
         L = IWRK(J)
         IF (J .NE. L) THEN
            DO 70 K = 1,NVAR
               TMP = UBX2(K,L)
               UBX2(K,L) = UBX2(K,J)
               UBX2(K,J) = TMP
   70       CONTINUE
            ITMP    = IWRK(L)
            IWRK(L) = IWRK(J)
            IWRK(J) = ITMP
         ENDIF
   80 CONTINUE
C
      RETURN
C
C.....End of Orient
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE PROJ( NVAR,MDIM,X,Y,XC,UBXC,LDU )
C
      INTEGER NVAR,MDIM,LDU
      DOUBLE PRECISION X(*),Y(*),XC(*),UBXC(LDU,*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     For an affine sub-space V = XC + span(UBXC) of R^NVAR
C     spanned by the orthonormal columns of the NVAR x MDIM array
C     UBXC and any point X in R^NVAR, the routine computes the
C     orthogonal projection of X - XC onto V in local coordinates
C     with respect to the given basis; that is, the point
C
C              Y = UBXC^T * (X - XC)
C
C  Variables in the calling sequence
C  ---------------------------------
C     NVAR   I   IN   dimension of the ambient space
C     MDIM   I   IN   dimension of the sffime subspace
C     X      D   IN   array of dimension NVAR, the given point
C                     in global coordinates
C     Y      D   OUT  array of dimension MDIM, the computed point
C                     in local coordinates
C     XC     D   IN   array of dimension NVAR, the center point of
C                     the local basis
C     UBXC   D   IN   array of dimension NVAR x MDIM, containing
C                     the orthonormal local basis
C     LDU    I   IN   row dimension of UBXC, LDU >= NVAR
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Parameter
C
      DOUBLE PRECISION ZER
      PARAMETER( ZER=0.0D0 )
C
C.....Local variables
C
      INTEGER I,J
      DOUBLE PRECISION SUM
C
C.......................Executable statements..........................
C
      DO 20 J = 1,MDIM
         SUM = ZER
         DO 10 I = 1,NVAR
            SUM = SUM + UBXC(I,J)*(X(I) - XC(I))
 10      CONTINUE
         Y(J) = SUM
 20   CONTINUE
C
      RETURN
C
C.....End of PROJ
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE SENMAP( NVAR,MDIM,UTAN,LDU,LPARA,SEN,LDS,
     &                   WRKMAT,LDW,WRK,IWRK,IER )
C
      INTEGER NVAR,MDIM,LDU,LDS,LDW,LPARA(*),IWRK(*),IER
      DOUBLE PRECISION UTAN(LDU,*),SEN(LDS,*)
      DOUBLE PRECISION WRKMAT(LDW,*),WRK(*)
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C   Let UTAN be the tangential local coordinate basis at some point
C   x on the manifold M and LPARA an MDIM-dimensional vector
C   specifying the indices of the variables in the natural parameter
C   basis TL of the problem. The routine computes the sensitivity map
C   SEN = UNL*(UL)^-1 at the point x. Here UL is the matrix formed
C   from the rows of UTAN specified by the indices given in LPARA,
C   and UNL the matrix of the remaining rows of UTAN.
C
C  Variables in the calling sequence
C  ---------------------------------
C   NVAR    I   IN   dimension of the ambient space
C   MDIM    I   IN   dimension of the manifold
C   UTAN    D   IN   array of dimension NVAR x MDIM, the tangent
C                    basis at some point of the manifold
C   LDU     I   IN   row dimension of UTAN, LDU >= NVAR
C   LPARA   I   IN   vector of dimension MDIM with the indices of
C                    the natural parameter variables
C   SEN     D   OUT  array of dimension MDIM x MDIM containing the
C                    computed sensitivity map of M at the particular
C                    point
C   LDS     I   IN   row dimension of SENMAP, LDS >= MDIM
C   WRKMAT  D   WRK  work array of dimension MDIM x MDIM
C   LDW     I   IN   row dimension of WRKMAT, LDW >= MDIM
C   WRK     D   WRK  work array of dimension MDIM
C   IWRK    I   WRK  work array of dimension MDIM
C   IER     I   OUT  error indicator
C                    IER =  0 no error
C                    IER =  1 sensivity is undefined
C                    IER = -1 error in the solution of the linear
C                             system with the matrix UL
C                    IER = -2 error in the parameter array LPARA
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL LUF, LUS1, MSGPRT
C
C.....Local variables
C
      INTEGER I,J,K
C
C.....Parameter
C
      CHARACTER*(*)LNAME
      PARAMETER( LNAME='SENMAP')
C
C.......................Executable statements..........................
C
C.....Form the matrix UL in WRKMAT
C
      DO 20 I = 1, MDIM
         K = LPARA(I)
         DO 10 J = 1, MDIM
            WRKMAT(J,I) = UTAN(K,J)
   10    CONTINUE
   20 CONTINUE
C
C.....Compute the LU-factorization
C
      CALL LUF( MDIM,WRKMAT,LDW,IWRK,IER )
      IF( IER .NE. 0 )THEN
         IF( IER .LT. 0 )THEN
            CALL MSGPRT(LNAME,'Input error for LU-factorization')
            IER = -1
         ELSE
            CALL MSGPRT(LNAME,'Numerically singular matrix')
            IER = 1
         ENDIF
         RETURN
      ENDIF
C
C.....Loop over the rows of UTAN not in LPARA
C
      K = 0
      DO 60 J = 1,NVAR
         DO 30 I = 1,MDIM
            IF( LPARA(I) .EQ. J) GOTO 60
   30    CONTINUE
         K = K+1
         DO 40 I = 1,MDIM
            WRK(I) = UTAN(K,I)
   40    CONTINUE
         CALL LUS1( MDIM,WRKMAT,LDW,IWRK,WRK,IER )
         IF( IER .NE. 0 )THEN
           IF( IER .LT. 0 )THEN
              CALL MSGPRT(LNAME,
     &              'Input error in solving the linear system')
              IER = -1
           ELSE
              CALL MSGPRT(LNAME,'Numerically singular matrix')
              IER = 1
           ENDIF
           RETURN
         ENDIF
         DO 50 I = 1,MDIM
            SEN(K,I) = WRK(I)
   50    CONTINUE
   60 CONTINUE
C
      IER = 0
      IF( K + MDIM .NE. NVAR )THEN
         CALL MSGPRT('SENMAP','Error in the parameter array LPARA')
         IER = -1
      ENDIF
C
      RETURN
C
C.....End of SENMAP
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE SENSNR( MDIM,UTAN,LDU,LPARA,SENRM,
     &                   WRKMAT,LDW,WRK1,WRK2,WRK3,
     &                   EPMACH,SAFMIN,IER )
C
      INTEGER MDIM,LDU,LDW,LPARA(*),IER
      DOUBLE PRECISION UTAN(LDU,*),SENRM
      DOUBLE PRECISION WRKMAT(LDW,*),WRK1(*),WRK2(*),WRK3(*)
      DOUBLE PRECISION EPMACH,SAFMIN
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C   Let UTAN be the tangential local coordinate basis at some point
C   x on the manifold M and LPARA an MDIM-dimensional vector
C   specifying the indices of the variables in the natural parameter
C   basis TL of the problem. The routine computes the Euclidean norm
C   of the sensitivity map at the point x; that is,
C
C      SENRM = || S ||  =  dist/ sqrt[ 1 - dist^2]
C
C   Here the distance dist = dist[rge(TL),rge(UTAN)] between the
C   local coordinate spaces is computed from the largest singular
C   value sigma of TL^T * UTAN as
C
C      dist = sin( arcos(sigma) )
C
C   For reference see
C
C      W. C. Rheinboldt, On the sensitivity of solutions of
C      parameterized equations, SIAM J. Num. Anal. 30, 1993,305-320
C
C  Variables in the calling sequence
C  ---------------------------------
C     MDIM    I   IN   dimension of the manifold
C     UTAN    D   IN   array of dimension NVAR x MDIM, the tangent
C                      basis at a point of the manifold
C     LDU     I   IN   row dimension of UTAN, LDU >= NVAR
C     LPARA   I   IN   vector of dimension MDIM with the indices of
C                      the natural parameter variables
C     SENRM   D   OUT  the computed Euclidean norm of the sensitivity
C                      map at the particular point of the manifold
C     WRKMAT  D   WRK  work array of dimension MDIM x MDIM
C     LDW     I   IN   row dimension of WRKMAT, LDW >= MDIM
C     WRK1
C      - WRK3 D   WRK  three work array of dimension MDIM
C     EPMACH  D   IN   machine constant, the smallest number such
C                      that 1.0D0 + EPMACH .NE. 1.0D0
C     IER     I   OUT  error indicator
C                      IER =  0 no error
C                      IER = -1 error returns from subroutine
C                      IER = 1  the value of SENSNR is infinite
C                               in machine arithmetic
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL BIDIA, SVD, MSGPRT
C
C.....Functions called
C
      DOUBLE PRECISION SQRT
C
C.....Local variables
C
      INTEGER I,J,JOB,K
      DOUBLE PRECISION QDUM(1,1),TMP
C
C.....Parameters
C
      INTEGER LDD
      DOUBLE PRECISION ONE, ZER
      CHARACTER*(*)LNAME
      PARAMETER( ONE = 1.0D0, ZER=0.0D0, LDD=1, LNAME='SENSNR' )
C
C.......................Executable statements..........................
C
C.....Generate the matrix TL * UTAN in WRKMAT
C
      DO 20 I = 1, MDIM
         K = LPARA(I)
         DO 10 J = 1, MDIM
            WRKMAT(I,J) = UTAN(K,J)
   10    CONTINUE
   20 CONTINUE
C
C.....Bidiagonalize WRKMAT
C
      JOB = 0
      CALL BIDIA( MDIM,MDIM,WRKMAT,LDW,WRK1,WRK2,QDUM,LDD,
     &            WRK3,SAFMIN,JOB,IER )
      IF( IER .NE. 0 )THEN
         CALL MSGPRT(LNAME,'Error in bidiagonaization')
         IER = -1
         RETURN
      ENDIF
C
C.....Compute the singular values
C
      CALL SVD( MDIM,MDIM,WRK1,WRK2,QDUM,LDD,QDUM,LDD,
     &          EPMACH,SAFMIN,JOB,IER )
      IF( IER .NE. 0 )THEN
         CALL MSGPRT(LNAME,'Error in SVD decomposition')
         IER = -1
         RETURN
      ENDIF
C
C.....Evaluate SENSRM
C
      IER = 0
      TMP = WRK1(MDIM)
      IF( TMP .EQ. ZER )THEN
         IER = 1
      ELSE
         TMP = ONE/(TMP*TMP)
         IF( TMP .LE. EPMACH )THEN
            SENRM = ZER
         ELSE
            SENRM = SQRT(TMP - ONE)
         ENDIF
      ENDIF
C
      RETURN
C
C.....End of SENSNR
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE TPHI( TASK,NVAR,MDIM,Y,X,FX,XC,DFXC,LDF,
     &                 TAUXC,JPXC,UBXC,LDU,WRK,EPMACH,ISTEP )
C
      CHARACTER*6 TASK
      INTEGER NVAR,MDIM,LDF,LDU,ISTEP,JPXC(*)
      DOUBLE PRECISION Y(*),X(*),FX(*)
      DOUBLE PRECISION XC(*),DFXC(LDF,*),TAUXC(*),UBXC(LDU,*)
      DOUBLE PRECISION WRK(*),EPMACH
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     Given a tangential, coordinate system at the point XC defined by
C     the arrays [DFXC,TAUXC,JPXC,UBXC], the point X on M in global
C     coordinates is computed which has the given vector Y as its local
C     representation in the coordinate system.
C
C     COBAS generates the arrays DFXC,TAUC,JPXC containing the
C     LQ-factorization
C                       P*DF(XC) = (L,0)*Q
C     of the Jacobian DF(XC). The chord Newton method
C
C                       X = UGLOB(Y)
C                       Until convergence do
C                            solve L*Z = P * F(X) for Z
C                            X := X - Q*(Z)
C                                       (0)
C     is used to compute the (unique) point X on M with local
C     coordinate Y. The process converges if ||Y|| is sufficiently
C     small; that is, if X is sufficiently close to the center point
C     XC of the local coordinate system.
C
C     The routine uses reverse communication to access the user
C     subroutine for the function F.  This is done under the control
C     of the character variable TASK. For the initial call set
C     TASK = 'START'. The routine returns repeatedly with
C     TASK = 'EVAL' signifying a request for the evaluation of
C     F at the global point X.
C     Return to the routine with the computed vector FX = F(X) but
C     no changes of TASK or of any other variable.
C     At the final return, TASK may have either one of the
C     following values:
C
C     1. TASK = 'DONE'    process converged
C     2. TASK = 'DIVERG'  divergence detected
C     3. TASK = 'STPCNT'  maximal number of steps NMAX exceeded
C     4. TASK = 'SINGUL'  Singular iteration matrix
C     5. TASK = 'ERROR'   Wrong input
C
C     The failure-exits 2. and 3. signify that ||Y|| was not small
C     enough for the given local coordinate system. This may not be
C     a fatal error but may only require the choice of a smaller Y.
C     On the other hand the failure exits 4, and 5. are fatal.
C
C  Variables in the calling sequence
C  ---------------------------------
C     TASK   C   IN   Task identifier
C                     TASK = 'START'  start the process
C                     TASK = 'EVAL'   the function vector F(X) is
C                                     available in FX
C                OUT  TASK = 'EVAL'   compute FX = F(X) for given X
C                     TASK = 'DONE'   process converged
C                     TASK = 'DIVERG' divergence detected
C                     TASK = 'STPCNT' maximal number of steps exceeded
C                     TASK = 'SINGUL' singular iteration matrix
C                     TASK = 'ERROR'  wrong input
C     NVAR   I   IN   dimension of X
C     MDIM   I   IN   dimension of the manifold
C     Y      D   IN   array of dimension MDIM, the point in local
C                     coordinates
C     X      D   I/O  array of dimension NVAR, global point on M
C     FX     D   IN   array of dimension NALG for computed FX = F(X)
C     XC     D   IN   the center point on M of the local coordinate
C                     system
C     DFXC   D   IN   array of dimension LDF x NVAR containing the
C                     LQ factorization of DF(XC)
C     LDF    I   IN   leading dimension of DFXC, LDF .GE. NALG
C     TAUXC  D   IN   array of dimension NALG for the factors of
C                     the Householder reflectors
C     JPXC   D   IN   array of dimension NALG for the pivot sequence
C                     of the LQ-factorization at XC
C     UBXC  D   IN    array of dimension LDU x MDIM containing the
C                     orthogonal basis of the tangent space at XC
C     LDU    I   IN   leading dimension of UBXC, LDU .GE. NVAR
C     WRK    D   WK   work array of dimension NALG
C     EPMACH D   IN   machine constant, the smallest number such
C                     that 1.0D0 + EPMACH .NE. 1.0D0
C     ISTEP  I   OUT  number of iteration steps
C                     ISTEP = -2  Input error, no X returned
C                     ISTEP = -1  Y is small, X = XC is returned
C                     ISTEP = 0   the residual is small, X = UGLOB(Y)
C                                 is returned
C                     ISTEP > 0   number of iteration steps taken. When
C                                 an error condition occurred then
C                                 ISTEP has the value of the current
C                                 step
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....Subroutines called
C
      EXTERNAL GLOB,LQAS
C
C.....Function calls
C
      DOUBLE PRECISION DDIST2,SQRT
C
C.....Parameters
C
      INTEGER MAXSTP
      DOUBLE PRECISION FAC,FACLOW, SIXTN
      PARAMETER( FAC=2.0D0, FACLOW=1.05D0, SIXTN=1.6D1, MAXSTP=25 )
C
C.....Local variables
C
      INTEGER I,IER
C
C.....Variables saved between calls
C
      INTEGER NALG
      DOUBLE PRECISION FNRM,FNRML,ABSERR,RELERR,RELSIX
      DOUBLE PRECISION STEP,STEPL,TOLX,TOLSIX,XNRM
      SAVE NALG,FNRM,FNRML,ABSERR,RELERR,RELSIX,
     &     STEP,STEPL,TOLX,TOLSIX,XNRM
C
      DATA FNRM/0.0D0/, STEPL/0.0D0/
C
C.......................Executable statements..........................
C
C.....Check TASK
C
      IF (TASK .EQ. 'EVAL') THEN
         GOTO 100
C
      ELSEIF (TASK .EQ. 'START') THEN
C
C........If ||Y|| < EPMACH set X = XC and return
C
         STEP = DDIST2(MDIM,Y,1,Y,1,1)
         IF(STEP .LT. EPMACH) THEN
            DO 10 I = 1,NVAR
               X(I) = XC(I)
  10        CONTINUE
            ISTEP = -1
            GOTO 200
         ENDIF
C
C........Initialize data
C
         NALG   = NVAR - MDIM
         XNRM   = DDIST2(NVAR,XC,1,XC,1,1)
         RELERR = SQRT(EPMACH)
         ABSERR = EPMACH/SQRT(RELERR)
         TOLX   = RELERR*XNRM + ABSERR
         TOLSIX = SIXTN*TOLX
         RELSIX = SIXTN*RELERR
C
         ISTEP  = 0
C
C........Globalize the point Y
C
         CALL GLOB(NVAR,MDIM,Y,X,XC,UBXC,LDU)
C
C........Ask for evaluation of FX = F(X)
C
         TASK = 'EVAL'
         RETURN
      ELSE
C
C........TASK has unacceptable value
C
         ISTEP = -2
         GOTO 240
      ENDIF
C
C.....Iteration point
C
  100 CONTINUE
C
C.....Compute the function norm
C
      FNRML = FNRM
      FNRM  = DDIST2(NALG,FX,1,FX,1,1)
C
C.....Test for convergence
C
      IF (ISTEP .EQ. 0) THEN
         IF (FNRM .LT. ABSERR) GOTO 200
      ELSE
C
C........Strong acceptance test
C
         IF(FNRM .LE. RELERR .AND. STEP .LE. TOLX) GOTO 200
C
C........Weak acceptance tests
C
         IF(FNRM .LE. ABSERR .OR. STEP .LE. ABSERR) GOTO 200
         IF(ISTEP .GT. 1) THEN
            IF((FNRM+FNRML) .LE. RELERR
     &             .AND. STEP.LE.TOLSIX) GOTO 200
            IF(FNRM .LE. RELSIX
     &             .AND. (STEP+STEPL).LE.TOLX) GOTO 200
         ENDIF
C
C........Divergence test
C
         IF (ISTEP .EQ. 1) THEN
            IF (FNRM .GT. (FAC*FNRML+RELERR)) GOTO 210
         ELSEIF (ISTEP .EQ. 2) THEN
            IF(STEP .GT. (FAC*STEPL+RELERR)) GOTO 210
            IF (FNRM .GT. (FACLOW*FNRML+RELERR)) GOTO 210
         ELSE
            IF(STEP .GT. (FACLOW*STEPL+RELERR)) GOTO 210
            IF (FNRM .GT. (FACLOW*FNRML+RELERR)) GOTO 210
         ENDIF
      ENDIF
C
C.....Check for maximum number of steps
C
      IF (ISTEP .GT. MAXSTP) GOTO 220
C
C.....compute the solution of DF(XC)*Z = FX which is orthogonal
C.....to  ker DF(XC) and return the result in WRK
C
      CALL LQAS(NALG,NVAR,DFXC,LDF,JPXC,TAUXC,FX,WRK,IER)
      IF(IER .NE. 0) GOTO 230
C
C.....Form the next iterate
C
      STEPL = STEP
      STEP = DDIST2(NVAR,WRK,1,WRK,1,1)
      DO 110 I = 1,NVAR
         X(I) = X(I) - WRK(I)
  110 CONTINUE
C
      ISTEP = ISTEP + 1
      TASK = 'EVAL'
      RETURN
C
C.....Common returns
C
  200 CONTINUE
      TASK = 'DONE'
      RETURN
C
  210 TASK = 'DIVERG'
      RETURN
C
  220 TASK = 'STPCNT'
      RETURN
C
  230 TASK = 'SINGUL'
      RETURN
C
  240 TASK = 'ERROR'
      RETURN
C
C.....End of TPHI
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=
C
      SUBROUTINE TSFT( TASK,NVAR,MDIM,FTEN,FTMX,DFX,LDF,
     &                 TAUX,JPX,WRK,SAFMIN,IER )
C
      CHARACTER*6 TASK
      INTEGER NVAR,MDIM,JPX(*),LDF,IER
      DOUBLE PRECISION FTEN(*),FTMX
      DOUBLE PRECISION DFX(LDF,*),TAUX(*),WRK(*),SAFMIN
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C     For a tangential, coordinate system at the point X of M
C     defined by the arrays  [DFX, TAUX, JPX], and a local tangent
C     vector TY of M at X, the routine overwrites the array FTEN
C     with the component vector V_x(TY,TY) of the second fundamental
C     tensor at X and, in addition, returns its maximum norm in FTMX.
C
C     The tangent direction TY is not explicitly used. Instead,
C     upon input, the array FTEN is assumed to contain in its
C     first NALG = NVAR - MDIM components the second derivative
C     term D2F(X)(DX,DX) formed with the global direction
C     DX = UBX*TY where UBX is the tangent basis at X.
C
C     Call the routine in either one of the following two different
C     modes under control of the variable TASK:
C
C     TASK = 'FACTOR'
C       Call the routine with the Jacobian DF(X) at X stored in the
C       array DFX and the second derivative vector D2F(X)(DX,DX)
C       stored in FTEN. The subroutine LQF will be used to compute
C       the LQ factorization of the Jacobian in the block of
C       arrays [DFX, TAUX, JPX].
C
C     TASK = 'EVAL'
C       Call the routine with the LQ factorization of the Jacobian
C       DF(X) in the block of arrays [DFX, TAUX,JPX] and the
C       second derivative vector D2F(X)(DX,DX) stored in FTEN.
C
C     In either mode the routine uses then the subroutine LQAS to
C     compute the tensor component V_x(TY,TY) in FTEN.
C
C     TASK = 'FACTOR' is the default
C
C  Variables in the calling sequence
C  ---------------------------------
C     TASK   C   IN   Task identifier with the values
C                     TASK = 'FACTOR' (see above)
C                     TASK = 'EVAL'   (see above)
C     NVAR   I   IN   Dimension of the ambient space
C     MDIM   I   IN   Dimension of the manifold
C     FTEN   D   IN   Array of dimension NVAR, containing
C                     D2F(X)(DX,DX) in its first NALG = NVAR-MDIM
C                     components
C                OUT  The computed NVAR-dimensional component
C                     Vx(TY,TY) of the second fundamental tensor
C     FTMX   D   OUT  The maximum norm of FTEN
C     DFX    D        Array of dimension NALG x NVAR
C                IN   If Task = 'FACTOR' the Jacobian DF(X)
C                     IF Task = 'EVAL' the LQ-factored Jacobian
C                OUT  The LQ-factored Jacobian
C     LDF    I   IN   Leading dimension of DFX, LDF >= NVAR
C     TAUX   D        Array of dimension NALG (=NVAR-MDIM)
C                IN   Task = 'EVAL' The diagonal terms of the
C                     Householder transformations used in the
C                     LQ-factorization of DFX
C                OUT  The diagonal terms of the Householder
C                     transformations used in the LQ-factorization
C                     of DFX
C     JPX    D   IN   Array of dimension NALG for the pivot
C                     array of the LQ-factorization
C     WRK    D   WK   Work array of dimension NVAR
C     SAFMIN D   IN   Safe minimum such that
C                     1.0D0/SAFMIN does not overflow
C     IER    I   OUT  error indicator
C                     IER =  0 no error
C                     IER = -1 error return from LQAS
C
C234567--1---------2---------3---------4---------5---------6---------7-
C
C.....External subroutines called
C
      EXTERNAL LQAS,MSGPRT
C
C.....FORTRAN function called
C
      DOUBLE PRECISION ABS
C
C.....Parameter
C
      DOUBLE PRECISION ZER
      PARAMETER( ZER=0.0D0 )
      CHARACTER*(*) LNAME
      PARAMETER( LNAME = 'TSFT')
C
C.....Local variables
C
      INTEGER I,NALG
      DOUBLE PRECISION TMP
C
C.......................Executable statements..........................
C
      NALG = NVAR - MDIM
      IF( TASK .NE. 'EVAL' )THEN
         CALL LQF( NALG,NVAR,DFX,LDF,JPX,TAUX,WRK,SAFMIN,IER )
         IF( IER .NE. 0) THEN
            CALL MSGPRT(LNAME,
     &         'Error in the LQ-factorization of the Jacobian' )
            IER = -1
            RETURN
         ENDIF
      ENDIF
C
C.....Use LQAS to obtain the unique solution of DFX z = FTEN which
C.....is orthogonal to the tangent space at X, and store z in WRK
C
      CALL LQAS( NALG,NVAR,DFX,LDF,JPX,TAUX,FTEN,WRK,IER )
      IF( IER .NE. 0 )THEN
         CALL MSGPRT(LNAME,
     &            'Error in solving for the tensor component')
         IER = -1
         RETURN
      ENDIF
C
C.....Change the sign of WRK and compute its maximal component
C
      FTMX = ZER
      DO 10 I = 1,NVAR
         TMP = WRK(I)
         FTEN(I) = -TMP
         IF( ABS(TMP) .GT. FTMX )FTMX = ABS(TMP)
   10 CONTINUE
C
      IER = 0
      RETURN
C
C.....End of TSFT
C
      END
C
C234567==1=========2=========3=========4=========5=========6=========7=