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=