C ALGORITHM 777, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 23,NO. 4, December, 1997, P. 514--549. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Info # Src # This archive created: Fri Jul 31 08:42:45 1998 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'README.txt' then echo shar: will not over-write existing file "'README.txt'" else cat << \SHAR_EOF > 'README.txt' ! ! HOMPACK90 is a suite of FORTRAN 90 subroutines for solving nonlinear ! systems of equations by homotopy methods. There are subroutines for ! fixed point, zero finding, and general homotopy curve tracking problems, ! utilizing both dense and sparse Jacobian matrices, and implementing ! three different algorithms: ODE-based, normal flow, and augmented ! Jacobian. The (driver) subroutines called by the user are given in the ! table below, and are well documented internally. The user need not ! be concerned with any other subroutines in HOMPACK90. ! ! ! Problem type ! --------|--------|--------|--------|--------|--------| ! x = f(x) | F(x) = 0 |rho(a,lambda,x)=0| ! --------|--------|--------|--------|--------|--------| ! dense | sparse | dense | sparse | dense | sparse | Algorithm ! --------|--------|--------|--------|--------|--------|--------------------- ! FIXPDF | FIXPDS | FIXPDF | FIXPDS | FIXPDF | FIXPDS | ODE based ! --------|--------|--------|--------|--------|--------|--------------------- ! FIXPNF | FIXPNS | FIXPNF | FIXPNS | FIXPNF | FIXPNS | normal flow ! --------|--------|--------|--------|--------|--------|--------------------- ! FIXPQF | FIXPQS | FIXPQF | FIXPQS | FIXPQF | FIXPQS | augmented Jacobian ! --------|--------|--------|--------|--------|--------|--------------------- ! ! ! The sparse subroutines use either the packed skyline storage scheme ! standard in structural mechanics or the compressed sparse row storage ! format, but any sparse storage scheme can be used by replacing some of ! the low-level HOMPACK90 routines with user-written routines. The ! stepping subroutines STEP?? or the reverse call subroutines STEPNX and ! ROOTNX may be of interest to some users with special curve tracking ! needs. ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ORGANIZATIONAL DETAILS. HOMPACK90 is organized in two different ways: by ! algorithm/problem type and by subroutine level. There are three levels ! of subroutines. The top level consists of drivers, one for each problem ! type and algorithm type. Normally these drivers are called by the user, ! and the user need know nothing beyond them. They allocate storage for ! the lower level routines, and all the arrays are variable dimension, so ! there is no limit on problem size. The second subroutine level ! implements the major components of the algorithms such as stepping ! along the homotopy zero curve, computing tangents, and the end game for ! the solution at lambda = 1 . A sophisticated user might call these ! routines directly to have complete control of the algorithm, or for ! some other task such as tracking an arbitrary parametrized curve over ! an arbitrary parameter range. The lowest subroutine level handles the ! numerical linear algebra, and includes some LAPACK and BLAS routines. ! All the linear algebra and associated data structure handling are ! concentrated in these routines, so a user could incorporate his own ! data structures by writing his own versions of these low level routines. ! ! The organization of HOMPACK90 by algorithm/problem type is shown in the ! above table, which lists the driver name for each algorithm and problem ! type. Using brackets to indicate the three subroutine levels described ! above, the natural grouping of the HOMPACK90 routines is: ! ! [FIXPDF] [FODE, ROOT, SINTRP, STEPS] [DGEQPF] ! ! [FIXPDS] [FODEDS, ROOT, SINTRP, STEPDS] [GMFADS, GMRES, ! GMRILUDS, ILUFDS, ILUSOLVDS, MULTDS, MULT2DS, PCGDS, SOLVDS] ! ! [FIXPNF] [ROOTNF, STEPNF, TANGNF] [DGEQPF, DORMQR, ROOT] ! ! [FIXPNS] [ROOTNS, STEPNS, TANGNS] [GMFADS, GMRES, GMRILUDS, ! ILUFDS, ILUSOLVDS, MULTDS, MULT2DS, PCGDS, ROOT, SOLVDS] ! ! [FIXPQF] [ROOTQF, STEPQF, TANGQF] [DGEQRF, DORGQR, UPQRQF] ! ! [FIXPQS] [ROOTNS, STEPQS, TANGNS] [GMFADS, GMRES, GMRILUDS, ! ILUFDS, ILUSOLVDS, MULTDS, MULT2DS, PCGDS, ROOT, SOLVDS] ! ! [POLSYS1H] [FIXPNF, ROOTNF, STEPNF, TANGNF] ! [DGEQPF, DGEQRF, DORMQR, DIVP, FFUNP, GFUNP, HFUNP, HFUN1P, ! INITP, MULP, OTPUTP, POWP, RHO, RHOJAC, ROOT, SCLGNP, STRPTP] ! ! The LAPACK and BLAS subroutines used by HOMPACK90 are ! DCOPY, DDOT, DGEMM, DGEMV, DGEQPF, DGEQR2, DGEQRF, DGER, DLAIC1, ! DLAMCH, DLAPY2, DLARF, DLARFB, DLARFG, DLARFT, DNRM2, DORG2R, DORGQR, ! DORM2R, DORMQR, DSCAL, DSWAP, DTPMV, DTPSV, DTRMM, DTRMV, DTRSV, ! IDAMAX, ILAENV, LSAME, XERBLA. ! ! The user written subroutines, of which exactly two must be supplied ! depending on the driver chosen, are F, FJAC, FJACS, RHO, RHOA, RHOJAC, ! and RHOJS. These external subroutines must conform to the interfaces ! contained in the module HOMOTOPY. The module REAL_PRECISION contains ! machine dependent constants, which must be changed appropriately before ! compilation. The module HOMPACK90_GLOBAL contains global storage, and ! must be used by the user written subroutines. ! ! Testing and installation: HOMPACK90 consists of 4 modules---HOMOTOPY ! (contains interfaces for the user written external subroutines), ! HOMPACK90 (encapsulates all the drivers), HOMPACK90_GLOBAL (global ! dynamic storage), REAL_PRECISION (defines precision of all reals)---and ! external subroutines, all contained in the files hompack90.f, ! blas1.f, blas2.f, blas3.f and lapack.f. ! The file template.f contains templates for the user written ! subroutines. There are three main programs driver[123].f for testing, ! with sample output given in the files RES[123]. driver1.f and ! driver3.f have no input files; driver2.f reads a data file DATA2 and ! writes the solution in a file RES2.OUT (for post-processing), since ! this is normally how the polynomial system driver POLSYS1H would be used. ! ! To test the dense (1), sparse (3), polynomial system (2) algorithms ! respectively in HOMPACK90, compile and link in order the files ! hompack90.f lapack.f blas1.f blas2.f blas3.f driver1.f ! (driver3.f, driver2.f respectively). ! The modules and external subroutines in hompack90.f, blas1.f, blas2.f, ! blas3.f and lapack.f (BLAS and LAPACK routines) can be kept in module ! and object libraries and need not be recompiled. ! ! ! Inquiries should be directed to Layne T. Watson, Departments of Computer ! Science and Mathematics, Virginia Polytechnic Institute & State ! University, Blacksburg, VA 24061-0106; (540) 231-7540; ltw@cs.vt.edu . ! SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'template.f' then echo shar: will not over-write existing file "'template.f'" else cat << \SHAR_EOF > 'template.f' ! Template for user written subroutines. ! ! USE statements for the modules REAL_PRECISION and HOMPACK90_GLOBAL ! should always be present in the user's subroutines. The user ! written subroutines must conform to the interfaces in the module ! HOMOTOPY. ! C C All data and subroutines defining the problem should be coded here. C SUBROUTINE F(X,V) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), DIMENSION(:), INTENT(IN):: X(:) REAL (KIND=R8), DIMENSION(:), INTENT(OUT):: V(:) C C Evaluate F(X) and return in the vector V . C RETURN END SUBROUTINE F SUBROUTINE FJAC(X,V,K) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), DIMENSION(:), INTENT(IN):: X(:) REAL (KIND=R8), DIMENSION(:), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C Return in V the Kth column of the Jacobian matrix of C F(X) evaluated at X . C RETURN END SUBROUTINE FJAC SUBROUTINE RHO(A,LAMBDA,X,V) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: IPAR, PAR ! for POLSYS1H only. REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) C C Evaluate RHO(A,LAMBDA,X) and return in the vector V . C C The following code is specifically for the polynomial system driver C POLSYS1H , and should be used verbatim with POLSYS1H . If the user is C calling FIXP?? or STEP?? directly, he must supply appropriate C replacement code here. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL C Force predicted point to have LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 NPOL=IPAR(1) CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,2*NPOL V(J)=PAR(IPAR(3 + (4-1)) + (J-1)) END DO C RETURN END SUBROUTINE RHO SUBROUTINE RHOA(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(OUT):: A(:) REAL (KIND=R8), INTENT(IN):: LAMBDA,X(:) C C Calculate and return in A the vector Z such that C RHO(Z,LAMBDA,X) = 0 . C RETURN END SUBROUTINE RHOA SUBROUTINE RHOJAC(A,LAMBDA,X,V,K) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: IPAR, PAR ! for POLSYS1H only. REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C Return in the vector V the Kth column of the Jacobian C matrix [D RHO/D LAMBDA, D RHO/DX] evaluated at the point C (A, LAMBDA, X). C C The following code is specifically for the polynomial system driver C POLSYS1H , and should be used verbatim with POLSYS1H . If the user is C calling FIXP?? or STEP?? directly, he must supply appropriate C replacement code here. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL,N2 NPOL=IPAR(1) N2=2*NPOL IF (K .EQ. 1) THEN C Force predicted point to have LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,N2 V(J)=PAR(IPAR(3 + (6-1)) + (J-1)) END DO RETURN ELSE DO J=1,N2 V(J)=PAR(IPAR(3 + (5-1)) + (J-1) + N2*(K-2)) END DO ENDIF C RETURN END SUBROUTINE RHOJAC SUBROUTINE FJACS(X) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL REAL (KIND=R8), DIMENSION(:), INTENT(IN):: X C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C RETURN END SUBROUTINE FJACS SUBROUTINE RHOJS(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL REAL (KIND=R8), INTENT(IN):: A(:),LAMBDA,X(:) C C If MODE = 1, C evaluate the N X N symmetric Jacobian matrix [D RHO/DX] at C (A,X,LAMBDA), and return the result in packed skyline storage format in C QRSPARSE. LENQR is the length of QRSPARSE, and ROWPOS contains the C indices of the diagonal elements of [D RHO/DX] within QRSPARSE. PP C contains -[D RHO/D LAMBDA] evaluated at (A,X,LAMBDA). Note the minus C sign in the definition of PP. The allocatable array COLPOS is not used C in this storage format. C C If MODE = 2, C evaluate the N X (N+1) Jacobian matrix [D RHO/DX, D RHO/DLAMBDA] at C (A,X,LAMBDA), and return the result in sparse row storage format in C QRSPARSE. LENQR is the length of QRSPARSE, ROWPOS contains the indices C of where each row begins within QRSPARSE, and COLPOS (of length LENQR) C contains the column indices of the corresponding elements in QRSPARSE. C Even if zero, the diagonal elements of the Jacobian matrix must be C stored in QRSPARSE. The allocatable array PP is not used in this C storage format. C RETURN END SUBROUTINE RHOJS SHAR_EOF fi # end of overwriting check if test -f 'driver3.f' then echo shar: will not over-write existing file "'driver3.f'" else cat << \SHAR_EOF > 'driver3.f' C MAIN PROGRAM TO TEST FIXPQS, FIXPNS, AND FIXPDS; C C THIS PROGRAM TESTS THE HOMPACK ROUTINES FIXPQS, FIXPNS, AND C FIXPDS. C C THE OUTPUT FROM THIS ROUTINE SHOULD BE AS FOLLOWS, WITH THE C EXECUTION TIMES CORRESPONDING TO A DEC AXP 3000/600. C C TESTING FIXPQS WITH STORAGE MODE = 1 C C LAMBDA = 1.00000000 FLAG = 1 33 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.119 ARC LENGTH = 1.274 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPNS WITH STORAGE MODE = 1 C C LAMBDA = 1.00000000 FLAG = 1 20 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.013 ARC LENGTH = 1.275 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPDS WITH STORAGE MODE = 1 C C LAMBDA = 1.00000000 FLAG = 1 70 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.031 ARC LENGTH = 1.281 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPQS WITH STORAGE MODE = 2 C C LAMBDA = 1.00000000 FLAG = 1 33 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.015 ARC LENGTH = 1.274 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPNS WITH STORAGE MODE = 2 C C LAMBDA = 1.00000000 FLAG = 1 20 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.011 ARC LENGTH = 1.275 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPDS WITH STORAGE MODE = 2 C C LAMBDA = 1.00000000 FLAG = 1 70 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.022 ARC LENGTH = 1.281 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C MODULE SWITCH C C ROWSET IS USED TO INITIALIZE SPARSE MATRIX DATA STRUCTURES ONLY C ONCE, AFTER THEY ARE ALLOCATED. C LOGICAL:: ROWSET END MODULE SWITCH ! PROGRAM TESTS USE SWITCH USE REAL_PRECISION, ONLY : R8 USE HOMPACK90, ONLY : FIXPDS, FIXPNS, FIXPQS IMPLICIT NONE INTEGER, PARAMETER:: N=8, NDIMA=8 REAL (KIND=R8):: A(N),ANSAE,ANSRE,ARCAE,ARCRE, & ARCLEN,DTIME,SSPAR(8),Y(N+1) INTEGER:: IFLAG,II,J,LENQR,MODE,NFE,NP1,TIMENEW(8), & TIMEOLD(8),TRACE CHARACTER (LEN=6):: NAME ! If using a subroutine library of the HOMPACK90 subroutines rather than ! the MODULE HOMPACK90 (as above), then the following INTERFACE ! statements are necessary. ! INTERFACE ! SUBROUTINE FIXPDS(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA, ! & NFE,ARCLEN,MODE,LENQR) ! USE REAL_PRECISION ! INTEGER, INTENT(IN)::LENQR,MODE,N,NDIMA,TRACE ! REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y ! INTEGER, INTENT(IN OUT)::IFLAG ! REAL (KIND=R8), INTENT(IN OUT)::ARCTOL,EPS ! INTEGER, INTENT(OUT)::NFE ! REAL (KIND=R8), INTENT(OUT)::ARCLEN ! END SUBROUTINE FIXPDS C ! SUBROUTINE FIXPNS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, ! & NFE,ARCLEN,MODE,LENQR,SSPAR) ! USE REAL_PRECISION ! INTEGER, INTENT(IN)::LENQR,MODE,N,TRACE ! REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y ! INTEGER, INTENT(IN OUT)::IFLAG ! REAL (KIND=R8), INTENT(IN OUT)::ANSAE,ANSRE,ARCAE,ARCRE, ! & SSPAR(8) ! INTEGER, INTENT(OUT)::NFE ! REAL (KIND=R8), INTENT(OUT)::ARCLEN ! END SUBROUTINE FIXPNS C ! SUBROUTINE FIXPQS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, ! & NFE,ARCLEN,MODE,LENQR,SSPAR) ! USE REAL_PRECISION ! INTEGER, INTENT(IN)::LENQR,MODE,N,TRACE ! REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y ! INTEGER, INTENT(IN OUT)::IFLAG ! REAL (KIND=R8), INTENT(IN OUT)::ANSAE,ANSRE,ARCAE,ARCRE, ! & SSPAR(4) ! INTEGER, INTENT(OUT)::NFE ! REAL (KIND=R8), INTENT(OUT)::ARCLEN ! END SUBROUTINE FIXPQS ! END INTERFACE C C TEST EACH OF THE THREE ALGORITHMS WITH BOTH STORAGE MODES. DO MODE=1,2 SELECT CASE (MODE) CASE(1) LENQR=18 CASE(2) LENQR=36 END SELECT DO II=1,3 C C DEFINE ARGUMENTS FOR CALL TO HOMPACK PROCEDURE. C NP1=N+1 ARCRE=0.5D-4 ARCAE=0.5D-4 ANSRE=1.0D-12 ANSAE=1.0D-12 TRACE=0 SSPAR=0.0 IFLAG=-MODE Y(1:N)=0.5_R8 IF(IFLAG .EQ. -2) A=Y(1:N) ROWSET = .FALSE. C C GET CURRENT DATE AND TIME. C CALL DATE_AND_TIME(VALUES=TIMEOLD) C C CALL TO HOMPACK ROUTINE. C IF (II .EQ. 1) THEN NAME='FIXPQS' CALL FIXPQS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, & NFE,ARCLEN,MODE,LENQR,SSPAR) ELSE IF (II .EQ. 2) THEN NAME='FIXPNS' CALL FIXPNS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, & NFE,ARCLEN,MODE,LENQR,SSPAR) ELSE NAME='FIXPDS' CALL FIXPDS(N,Y,IFLAG,ARCRE,ANSRE,TRACE,A,NDIMA,NFE, & ARCLEN,MODE,LENQR) END IF C C CALCULATE EXECUTION TIME. C CALL DATE_AND_TIME(VALUES=TIMENEW) IF (TIMENEW(8) .LT. TIMEOLD(8)) THEN TIMENEW(8)=TIMENEW(8)+1000 TIMENEW(7)=TIMENEW(7)-1 ENDIF IF (TIMENEW(7) .LT. TIMEOLD(7)) THEN TIMENEW(7)=TIMENEW(7)+60 TIMENEW(6)=TIMENEW(6)-1 ENDIF IF (TIMENEW(6) .LT. TIMEOLD(6)) THEN TIMENEW(6)=TIMENEW(6)+60 TIMENEW(5)=TIMENEW(5)-1 ENDIF IF (TIMENEW(5) .LT. TIMEOLD(5)) TIMENEW(5)=TIMENEW(5)+24 DTIME=DOT_PRODUCT(TIMENEW(5:8)-TIMEOLD(5:8), & (/3600000,60000,1000,1/) )/1000.0 C WRITE (6,45) NAME, MODE 45 FORMAT (//,7X,'TESTING',1X,A6,' WITH STORAGE MODE =',I2) WRITE (6,50) Y(NP1),IFLAG,NFE,DTIME,ARCLEN,(Y(J),J=1,N) 50 FORMAT(/' LAMBDA =',F11.8,' FLAG =',I2,I8,' JACOBIAN ', & 'EVALUATIONS',/,1X,'EXECUTION TIME(SECS) =',F10.3,4X, & 'ARC LENGTH =',F10.3/(1X,4ES16.8)) END DO END DO STOP END PROGRAM TESTS ! ! SAMPLE USER WRITTEN HOMOTOPY SUBROUTINES FOR TESTING FIXP*S. ! SUBROUTINE F(X,V) C*********************************************************************** C C SUBROUTINE F(X,V) COMPUTES F AT THE POINT X, RETURNING THE VALUE IN V. C C*********************************************************************** USE REAL_PRECISION, ONLY : R8 IMPLICIT NONE REAL (KIND=R8), INTENT(IN):: X(:) REAL (KIND=R8), INTENT(OUT):: V(:) V(1)=X(1)**3+6.0*X(2)*X(3)-1+2.0*X(1) V(2)=6.0*X(1)*X(3)+X(2)**4*X(5)-1+3.0*X(2) V(3)=6.0*X(1)*X(2)+X(3)*X(5)-1+4.0*X(3) V(4)=X(4)**3*X(8)-1+2.0*X(4) V(5)=X(2)**5/5.0 + X(3)**2/2.0 + X(8)*X(5)-1+3.0*X(5) V(6)=X(6)*X(8)-1+4.0*X(6) V(7)=X(7)**2*X(8)**3-1+2.0*X(7) V(8)=X(4)**4/4.0 + X(5)**2/2.0 + X(6)**2/2.0 + X(7)**3* & X(8)**2-1+3.0*X(8) RETURN END SUBROUTINE F C SUBROUTINE FJACS(X) C****************************************************************** C C SUBROUTINE FJACS(X) COMPUTES THE JACOBIAN MATRIX OF F AT THE POINT C X, RETURNING THE JACOBIAN MATRIX IN PACKED SKYLINE FORM (MODE=1) C IN THE ARRAYS QRSPARSE AND ROWPOS . C C***************************************************************** USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: QRSPARSE, ROWPOS, COLPOS USE SWITCH REAL (KIND=R8), INTENT(IN):: X(:) C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C INTEGER:: N N=SIZE(X) IF (.NOT. ROWSET) THEN ROWSET=.TRUE. ROWPOS(1:N+1) = (/ 1,2,4,7,8,12,13,14,19 /) END IF QRSPARSE(1)=3.0*X(1)**2+2.0 QRSPARSE(2)=4.0*X(2)**3*X(5)+3.0 QRSPARSE(3)=6.0*X(3) QRSPARSE(4)=X(5)+4.0 QRSPARSE(5)=6.0*X(1) QRSPARSE(6)=6.0*X(2) QRSPARSE(7)=3.0*X(4)**2*X(8)+2.0 QRSPARSE(8)=X(8)+3.0 QRSPARSE(9)=0.0 QRSPARSE(10)=X(3) QRSPARSE(11)=X(2)**4 QRSPARSE(12)=X(8)+4.0 QRSPARSE(13)=2.0*X(7)*X(8)**3+2.0 QRSPARSE(14)=2.0*X(7)**3*X(8)+3.0 QRSPARSE(15)=3.0*X(7)**2*X(8)**2 QRSPARSE(16)=X(6) QRSPARSE(17)=X(5) QRSPARSE(18)=X(4)**3 RETURN END SUBROUTINE FJACS C SUBROUTINE RHO(A,LAMBDA,X,V) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) INTERFACE SUBROUTINE F(X,V) USE REAL_PRECISION REAL (KIND=R8), DIMENSION(:), INTENT(IN):: X REAL (KIND=R8), DIMENSION(:), INTENT(OUT):: V END SUBROUTINE F END INTERFACE C C EVALUATE RHO(A,LAMBDA,X) AND RETURN IN THE VECTOR V . C INTEGER:: N N=SIZE(X) CALL F(X(1:N), V(1:N)) V(1:N) = LAMBDA*V(1:N) + (1.0 - LAMBDA)*(X(1:N) - A(1:N)) C RETURN END SUBROUTINE RHO C SUBROUTINE RHOA(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(OUT):: A(:) REAL (KIND=R8), INTENT(IN):: LAMBDA,X(:) INTERFACE SUBROUTINE F(X,V) USE REAL_PRECISION REAL (KIND=R8), DIMENSION(:), INTENT(IN):: X REAL (KIND=R8), DIMENSION(:), INTENT(OUT):: V END SUBROUTINE F END INTERFACE C C CALCULATE AND RETURN IN A THE VECTOR Z SUCH THAT C RHO(Z,LAMBDA,X) = 0 . C N=NDIMA FOR THIS TEST PROBLEM. C INTEGER:: N N=SIZE(X) CALL F(X(1:N),A(1:N)) A(1:N)=LAMBDA*A(1:N)/(1.0 - LAMBDA) + X(1:N) RETURN END SUBROUTINE RHOA C SUBROUTINE RHOJS(A,LAMBDA,X) C***************************************************************** C C Subroutine RHOJS(A,LAMBDA,X) computes the Jacobian matrix of C rho(a,x,lambda) = lambda*F(x) + (1 - lambda)*(x - a) at the C point (A,X,LAMBDA), returning the Jacobian matrix in sparse row C storage format (MODE = 2) in the arrays QRSPARSE, ROWPOS, and C COLPOS. C C***************************************************************** USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: QRSPARSE, ROWPOS, COLPOS USE SWITCH REAL (KIND=R8), INTENT(IN):: A(:),LAMBDA,X(:) INTERFACE SUBROUTINE F(X,V) USE REAL_PRECISION REAL (KIND=R8), DIMENSION(:), INTENT(IN):: X REAL (KIND=R8), DIMENSION(:), INTENT(OUT):: V END SUBROUTINE F END INTERFACE C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C C--------------------------------------------------------------------- C [QRSPARSE] = [ LAMBDA*DF(X) + (1 - LAMBDA)*I | F(X) - X + A ] . C--------------------------------------------------------------------- C C LOCAL VARIABLES C INTEGER, PARAMETER:: N=8 INTEGER:: J, JPOS, ELEM(N) = (/4,9,14,18,24,27,30,36/) REAL (KIND=R8):: DRHODL(N) C IF (.NOT. ROWSET) THEN ROWSET=.TRUE. ROWPOS(1:N+1) = (/ 1,5,10,15,19,25,28,31,37 /) COLPOS(1:36) = (/1,2,3,9,1,2,3,5,9,1,2,3,5,9,4,5,8,9, & 2,3,4,5,8,9,6,8,9,7,8,9,4,5,6,7,8,9/) END IF C QRSPARSE = 0.0 C ROW 1. QRSPARSE(1)=3.0*X(1)**2+2.0 QRSPARSE(2)=6.0*X(3) QRSPARSE(3)=6.0*X(2) C ROW 2. QRSPARSE(5)=6.0*X(3) QRSPARSE(6)=4.0*X(2)**3*X(5)+3.0 QRSPARSE(7)=6.0*X(1) QRSPARSE(8)=X(2)**4 C ROW 3. QRSPARSE(10)=6.0*X(2) QRSPARSE(11)=6.0*X(1) QRSPARSE(12)=X(5)+4.0 QRSPARSE(13)=X(3) C ROW 4. QRSPARSE(15)=3.0*X(4)**2*X(8)+2.0 QRSPARSE(16)=0.0 QRSPARSE(17)=X(4)**3 C ROW 5. QRSPARSE(19)=X(2)**4 QRSPARSE(20)=X(3) QRSPARSE(21)=0.0 QRSPARSE(22)=X(8)+3.0 QRSPARSE(23)=X(5) C ROW 6. QRSPARSE(25)=X(8)+4.0 QRSPARSE(26)=X(6) COLPOS(25)=6 COLPOS(26)=8 COLPOS(27)=9 C ROW 7. QRSPARSE(28)=2.0*X(7)*X(8)**3+2.0 QRSPARSE(29)=3.0*X(7)**2*X(8)**2 C ROW 8. QRSPARSE(31)=X(4)**3 QRSPARSE(32)=X(5) QRSPARSE(33)=X(6) QRSPARSE(34)=3.0*X(7)**2*X(8)**2 QRSPARSE(35)=2.0*X(7)**3*X(8)+3.0 C QRSPARSE=LAMBDA*QRSPARSE C C FIND INDEX JPOS OF DIAGONAL ELEMENT IN JTH ROW OF QR. C DO J=1,N JPOS=ROWPOS(J) DO IF (COLPOS(JPOS) .EQ. J) EXIT JPOS=JPOS+1 END DO QRSPARSE(JPOS) = QRSPARSE(JPOS) + 1.0 - LAMBDA END DO C C INITIALIZE (N+1)ST COLUMN. C CALL F(X(1:N),DRHODL(1:N)) DRHODL = DRHODL - X(1:N) + A(1:N) QRSPARSE(ELEM) = DRHODL(1:8) C RETURN END SUBROUTINE RHOJS C ********************************************************************** C C THE REST OF THESE SUBROUTINES ARE NOT USED BY PROGRAM TESTS, AND ARE C INCLUDED HERE SIMPLY FOR COMPLETENESS AND AS TEMPLATES FOR THEIR USE. C SUBROUTINE FJAC(X,V,K) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(IN):: X(:) REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C RETURN IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF C F(X) EVALUATED AT X . C V(1)=X(1) ! INTENT(OUT) VARIABLE MUST BE DEFINED. RETURN END SUBROUTINE FJAC C SUBROUTINE RHOJAC(A,LAMBDA,X,V,K) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: IPAR, PAR REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C RETURN IN THE VECTOR V THE KTH COLUMN OF THE JACOBIAN C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT THE POINT C (A, LAMBDA, X). C C THE FOLLOWING CODE IS SPECIFICALLY FOR THE POLYNOMIAL SYSTEM DRIVER C POLSYS1H , AND SHOULD BE USED VERBATIM WITH POLSYS1H . IF THE USER IS C CALLING FIXP?? OR STEP?? DIRECTLY, HE MUST SUPPLY APPROPRIATE C REPLACEMENT CODE HERE. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL,N2 NPOL=IPAR(1) N2=2*NPOL IF (K .EQ. 1) THEN C FORCE PREDICTED POINT TO HAVE LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 ! CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,N2 V(J)=PAR(IPAR(3 + (6-1)) + (J-1)) END DO RETURN ELSE DO J=1,N2 V(J)=PAR(IPAR(3 + (5-1)) + (J-1) + N2*(K-2)) END DO ENDIF C RETURN END SUBROUTINE RHOJAC SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << \SHAR_EOF > 'driver2.f' C MAIN ROUTINE TO TEST POLSYS1H. C C THIS ROUTINE REQUIRES ONE INPUT FILE, DATA2 . C C A SAMPLE INPUT FILE AND ASSOCIATED OUTPUT ARE GIVEN C IN THE COMMENTS THAT FOLLOW. THIS SAMPLE PROBLEM IS C CITED IN THE HOMPACK REPORT. C C***** SAMPLE INPUT DATA (READ FROM THE FILE 'DATA2'): C ' TWO QUADRICS, NO SOLUTIONS AT INFINITY, TWO REAL SOLUTIONS.' C &PROBLEM C IFLGHM = 1 C IFLGSC = 1 C TOTDG = 4 C MAXT = 6 C EPSBIG = 1.D-04 C EPSSML = 1.D-14 C SSPAR(5) = 1.D+00 C NUMRR = 10 C N = 2 / C 00006 NUMTRM(1) C 00002 DEG(1,1,1) C 00000 DEG(1,2,1) C -.00098D+00 C 00000 DEG(1,1,2) C 00002 DEG(1,2,2) C 978000.D+00 C 00001 DEG(1,1,3) C 00001 DEG(1,2,3) C -9.8D+00 C 00001 DEG(1,1,4) C 00000 DEG(1,2,4) C -235.0D+00 C 00000 DEG(1,1,5) C 00001 DEG(1,2,5) C 88900.0D+00 C 00000 DEG(1,1,6) C 00000 DEG(1,2,6) C -1.000D+00 C 00006 NUMTRM(2) C 00002 DEG(2,1,1) C 00000 DEG(2,2,1) C -.0100D+00 C 00000 DEG(2,1,2) C 00002 DEG(2,2,2) C -.9840D+00 C 00001 DEG(2,1,3) C 00001 DEG(2,2,3) C -29.70D+00 C 00001 DEG(2,1,4) C 00000 DEG(2,2,4) C .00987D+00 C 00000 DEG(2,1,5) C 00001 DEG(2,2,5) C -.1240D+00 C 00000 DEG(2,1,6) C 00000 DEG(2,2,6) C -.2500D+00 C***** END OF SAMPLE INPUT DATA. C C***** ASSOCIATED SAMPLE OUTPUT (WRITTEN TO THE FILE 'OUTHP.DAT'): C ! POLSYS1H TEST ROUTINE 7/7/95 ! ! TWO QUADRICS, NO SOLUTIONS AT INFINITY, TWO REAL SOLUTIONS. ! ! IF IFLGHM=1, HOMOGENEOUS; IF IFLGHM=0, INHOMOGENEOUS; IFLGHM= 1 ! ! IF IFLGSC=1, SCLGNP USED; IF IFLGSC=0, NO SCALING; IFLGSC= 1 ! ! TOTDG= 4 MAXT= 6 ! ! EPSBIG, EPSSML = 1.000000000000000E-04 1.000000000000000E-14 ! ! SSPAR(5) = 1.000000000000000E+00 ! ! NUMBER OF EQUATIONS = 2 ! ! NUMBER OF RECALLS WHEN IFLAG=3: 10 ! ! ! ****** COEFFICIENT TABLEAU ****** ! ! NUMT( 1) = 6 ! KDEG( 1, 1, 1) = 2 ! KDEG( 1, 2, 1) = 0 ! COEF( 1, 1) = -9.800000000000000E-04 ! KDEG( 1, 1, 2) = 0 ! KDEG( 1, 2, 2) = 2 ! COEF( 1, 2) = 9.780000000000000E+05 ! KDEG( 1, 1, 3) = 1 ! KDEG( 1, 2, 3) = 1 ! COEF( 1, 3) = -9.800000000000001E+00 ! KDEG( 1, 1, 4) = 1 ! KDEG( 1, 2, 4) = 0 ! COEF( 1, 4) = -2.350000000000000E+02 ! KDEG( 1, 1, 5) = 0 ! KDEG( 1, 2, 5) = 1 ! COEF( 1, 5) = 8.890000000000000E+04 ! KDEG( 1, 1, 6) = 0 ! KDEG( 1, 2, 6) = 0 ! COEF( 1, 6) = -1.000000000000000E+00 ! ! NUMT( 2) = 6 ! KDEG( 2, 1, 1) = 2 ! KDEG( 2, 2, 1) = 0 ! COEF( 2, 1) = -1.000000000000000E-02 ! KDEG( 2, 1, 2) = 0 ! KDEG( 2, 2, 2) = 2 ! COEF( 2, 2) = -9.840000000000000E-01 ! KDEG( 2, 1, 3) = 1 ! KDEG( 2, 2, 3) = 1 ! COEF( 2, 3) = -2.970000000000000E+01 ! KDEG( 2, 1, 4) = 1 ! KDEG( 2, 2, 4) = 0 ! COEF( 2, 4) = 9.870000000000000E-03 ! KDEG( 2, 1, 5) = 0 ! KDEG( 2, 2, 5) = 1 ! COEF( 2, 5) = -1.240000000000000E-01 ! KDEG( 2, 1, 6) = 0 ! KDEG( 2, 2, 6) = 0 ! COEF( 2, 6) = -2.500000000000000E-01 ! ! ! ! IFLG1 = 11 ! ! PATH NUMBER = 1 ! ! FINAL VALUES FOR PATH ! ! ARCLEN = 1.005533190562901E+01 ! NFE = 53 ! IFLG2 = 1 ! REAL, FINITE SOLUTION ! LAMBDA = 1.000000000000003E+00 ! X( 1) = 2.342338519591276E+03 8.841149143431121E-13 ! X( 2) = -7.883448240941412E-01 -9.356862757018485E-16 ! ! X( 3) = -9.493594594086552E-03 -1.064475509002627E-03 ! ! ! PATH NUMBER = 2 ! ! FINAL VALUES FOR PATH ! ! ARCLEN = 1.721129286057142E+00 ! NFE = 37 ! IFLG2 = 1 ! COMPLEX, FINITE SOLUTION ! LAMBDA = 1.000000000000006E+00 ! X( 1) = 1.614785792344189E-02 1.684969554988811E+00 ! X( 2) = 2.679947396144760E-04 4.428029939736605E-03 ! ! X( 3) = -3.819489729424030E-01 3.720689434572830E-01 ! ! ! PATH NUMBER = 3 ! ! FINAL VALUES FOR PATH ! ! ARCLEN = 2.023295279367267E+00 ! NFE = 35 ! IFLG2 = 1 ! COMPLEX, FINITE SOLUTION ! LAMBDA = 1.000000000000000E+00 ! X( 1) = 1.614785792343521E-02 -1.684969554988812E+00 ! X( 2) = 2.679947396144598E-04 -4.428029939736611E-03 ! ! X( 3) = -3.293704938476598E-01 5.566197755230126E-01 ! ! ! PATH NUMBER = 4 ! ! FINAL VALUES FOR PATH ! ! ARCLEN = 4.163266156958467E+00 ! NFE = 46 ! IFLG2 = 1 ! REAL, FINITE SOLUTION ! LAMBDA = 9.999999999999998E-01 ! X( 1) = 9.089212296153869E-02 1.153793567884107E-16 ! X( 2) = -9.114970981974997E-02 1.887399041592030E-17 ! ! X( 3) = -5.736733957279616E-02 1.362436637092185E-01 ! ! ! TOTAL NFE OVER ALL PATHS = 171 C C***** END OF ASSOCIATED SAMPLE OUTPUT. C C ************************************************************* C C PROGRAM DESCRIPTION: 1. READS IN DATA. C 2. GENERATES POLSYS1H INPUT. C 3. CALLS POLSYS1H. C 4. WRITES POLSYS1H OUTPUT. C C DIMENSIONS SHOULD BE SET AS FOLLOWS: C C DIMENSION NUMT(NN),COEF(NN,MMAXT),KDEG(NN,NN+1,MMAXT) C DIMENSION IFLG2(TTOTDG) C DIMENSION LAMBDA(TTOTDG),ROOTS(2,NN+1,TTOTDG),ARCLEN(TTOTDG), C & NFE(TTOTDG) C C WHERE: C N IS THE NUMBER OF EQUATIONS. NN .GE. N. C MAXT IS THE MAXIMUM NUMBER OF TERMS IN ANY ONE EQUATION. C MMAXT .GE. MAXT. C TOTDG IS THE TOTAL DEGREE OF THE SYSTEM. TTOTDG .GE. TOTDG. C C THIS TEST CODE HAS DIMENSIONS SET AS FOLLOWS: C C NN=10, MMAXT=30, TTOTDG=1024 C PROGRAM TESTP USE REAL_PRECISION, ONLY : R8 USE HOMPACK90, ONLY : POLSYS1H INTEGER, PARAMETER:: NN=10,MMAXT=30,TTOTDG=1024 INTEGER:: IFLG1,IFLG2(TTOTDG),IFLGHM,IFLGSC,ITOTIT,J,K, & KDEG(NN,NN+1,MMAXT),L,M,MAXT,N,NFE(TTOTDG),NP1,NT, & NUMRR,NUMT(NN),TOTDG REAL (KIND=R8):: ARCLEN(TTOTDG),COEF(NN,MMAXT),EPSBIG,EPSSML, & LAMBDA(TTOTDG),ROOTS(2,NN+1,TTOTDG),SSPAR(8) CHARACTER (LEN=72):: TITLE ! If using a subroutine library of the HOMPACK90 subroutines rather than ! the MODULE HOMPACK90 (as above), then the following INTERFACE ! statements are necessary. ! INTERFACE ! SUBROUTINE POLSYS1H(N,NUMT,COEF,KDEG,IFLG1,IFLG2,EPSBIG,EPSSML, ! & SSPAR,NUMRR,LAMBDA,ROOTS,ARCLEN,NFE) ! USE HOMOTOPY ! USE REAL_PRECISION ! INTEGER, INTENT(IN):: N,NUMT(:),KDEG(:,:,:),NUMRR ! REAL (KIND=R8), INTENT(IN):: COEF(:,:),EPSBIG,EPSSML ! INTEGER, INTENT(IN OUT):: IFLG1,IFLG2(:) ! REAL (KIND=R8), INTENT(IN OUT):: SSPAR(8) ! REAL (KIND=R8), INTENT(OUT):: LAMBDA(:),ROOTS(:,:,:),ARCLEN(:) ! INTEGER, INTENT(OUT):: NFE(:) ! END SUBROUTINE POLSYS1H ! END INTERFACE C NAMELIST /PROBLEM/ IFLGHM,IFLGSC,TOTDG,MAXT,N, & EPSBIG,EPSSML,SSPAR,NUMRR C OPEN (UNIT=7,FILE='DATA2',ACTION='READ',POSITION='REWIND', & DELIM='APOSTROPHE',STATUS='OLD') OPEN (UNIT=6,FILE='RES2.OUT',ACTION='WRITE', & DELIM='APOSTROPHE',STATUS='REPLACE') C SSPAR(1:8) = 0.0 READ (7,*) TITLE WRITE (6,10) TITLE 10 FORMAT(5X,'POLSYS1H TEST ROUTINE 7/7/95',//,A72) C READ (7, NML=PROBLEM) C WRITE (6,100) IFLGHM 100 FORMAT(/ &' IF IFLGHM=1, HOMOGENEOUS; IF IFLGHM=0, INHOMOGENEOUS; IFLGHM=' & ,I2) WRITE (6,102) IFLGSC 102 FORMAT(/ &' IF IFLGSC=1, SCLGNP USED; IF IFLGSC=0, NO SCALING; IFLGSC=',I5) WRITE (6,104) TOTDG,MAXT 104 FORMAT(/,' TOTDG=',I5,10X,'MAXT=',I5) WRITE (6,106) EPSBIG,EPSSML,SSPAR(5),N,NUMRR 106 FORMAT(/,' EPSBIG, EPSSML =',2ES22.14, & //,' SSPAR(5) =',ES22.14, & //,' NUMBER OF EQUATIONS =',I5, & //,' NUMBER OF RECALLS WHEN IFLAG=3:',I5) C NP1=N+1 C C NOTE THAT THE DEGREES OF VARIABLES IN EACH TERM OF EACH EQUATION C ARE DEFINED BY THE FOLLOWING INDEXING SCHEME: C C KDEG(J, L, K) C C ^ ^ ^ C C E V T C Q A E C U R R C A I M C T A C I B C O L C N E C WRITE(6,200) 200 FORMAT(//,' ****** COEFFICIENT TABLEAU ******') KDEG = 0 !SET UNUSED DEGREES TO ZERO EQN: DO J=1,N READ (7,1000) NUMT(J) WRITE (6,210) J,NUMT(J) 210 FORMAT(/,' NUMT(',I2,') =',I5) NT=NUMT(J) TERMS: DO K=1,NT VARS: DO L=1,N READ (7,1000) KDEG(J,L,K) WRITE (6,220) J,L,K,KDEG(J,L,K) 220 FORMAT(' KDEG(',I2,',',I2,',',I2,') =',I5) END DO VARS READ (7,2000) COEF(J,K) WRITE (6,230) J,K,COEF(J,K) 230 FORMAT(' COEF(',I2,',',I2,') =',ES22.14) END DO TERMS END DO EQN WRITE (6,FMT="(//)") C IFLG1=10*IFLGHM+IFLGSC DO M=1,TOTDG IFLG2(M)=-2 END DO CALL POLSYS1H(N,NUMT(1:N),COEF(1:N,1:MAXT), & KDEG(1:N,1:N+1,1:MAXT),IFLG1,IFLG2(1:TOTDG),EPSBIG,EPSSML, & SSPAR,NUMRR,LAMBDA(1:TOTDG),ROOTS(1:2,1:N+1,1:TOTDG), & ARCLEN(1:TOTDG),NFE(1:TOTDG)) C WRITE (6,240) IFLG1 240 FORMAT(' IFLG1 =',I5,/) ITOTIT = SUM(NFE(1:TOTDG)) DO M=1,TOTDG WRITE (6,260) M 260 FORMAT(' PATH NUMBER =',I5,//' FINAL VALUES FOR PATH'/) WRITE (6,280) ARCLEN(M) 280 FORMAT(' ARCLEN =',ES22.14) WRITE (6,290) NFE(M) 290 FORMAT(' NFE =',I5) WRITE (6,300) IFLG2(M) 300 FORMAT(' IFLG2 =',I3) C C DESIGNATE SOLUTIONS "REAL" OR "COMPLEX" C IF (ANY(ABS(ROOTS(2,1:N,M)) .GE. 1.0E-4)) THEN WRITE (6,779,ADVANCE='NO') 779 FORMAT(' COMPLEX, ') ELSE WRITE (6,780,ADVANCE='NO') 780 FORMAT(' REAL, ') END IF C C DESIGNATE SOLUTION "FINITE" OR "INFINITE" C IF (SUM(ABS(ROOTS(1:2,NP1,M))) .LT. 1.0E-6) THEN WRITE (6,781) 781 FORMAT('INFINITE SOLUTION') ELSE WRITE (6,782) 782 FORMAT('FINITE SOLUTION') END IF C WRITE (6,320) LAMBDA(M),(J,(ROOTS(L,J,M),L=1,2),J=1,N) 320 FORMAT(' LAMBDA =',ES22.14,/,(' X(',I2,') =',2ES22.14)) WRITE (6,330) NP1,ROOTS(1:2,NP1,M) 330 FORMAT(/,' X(',I2,') =',2ES22.14,//) END DO WRITE (6,400) ITOTIT 400 FORMAT(' TOTAL NFE OVER ALL PATHS = ',I10) STOP 1000 FORMAT(I5) 2000 FORMAT(ES22.14) END PROGRAM TESTP ! ! HOMOTOPY subroutines for the polynomial system driver POLSYS1H. ! These subroutines should be used verbatim with POLSYS1H for solving ! polynomial systems of equations. The polynomial coefficients, defined ! as input to POLSYS1H, are accessed by the routines here via the global ! arrays in HOMPACK90_GLOBAL. ! C ################################################################### C ONLY THE SUBROUTINES RHO AND RHOJAC ARE USED BY THE POLYNOMIAL C SYSTEM DRIVER POLSYS1H. ALL THE OTHER ROUTINES HERE ARE PROVIDED C SIMPLY AS TEMPLATES. C ################################################################### ! SUBROUTINE F(X,V) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(IN):: X(:) REAL (KIND=R8), INTENT(OUT):: V(:) C C EVALUATE F(X) AND RETURN IN THE VECTOR V . C V(1)=X(1) ! INTENT(OUT) VARIABLE MUST BE DEFINED. RETURN END SUBROUTINE F SUBROUTINE FJAC(X,V,K) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(IN):: X(:) REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C RETURN IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF C F(X) EVALUATED AT X . C V(1)=X(1) ! INTENT(OUT) VARIABLE MUST BE DEFINED. RETURN END SUBROUTINE FJAC SUBROUTINE RHO(A,LAMBDA,X,V) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: IPAR, PAR REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) C C EVALUATE RHO(A,LAMBDA,X) AND RETURN IN THE VECTOR V . C C THE FOLLOWING CODE IS SPECIFICALLY FOR THE POLYNOMIAL SYSTEM DRIVER C POLSYS1H , AND SHOULD BE USED VERBATIM WITH POLSYS1H . IF THE USER IS C CALLING FIXP?? OR STEP?? DIRECTLY, HE MUST SUPPLY APPROPRIATE C REPLACEMENT CODE HERE. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL C FORCE PREDICTED POINT TO HAVE LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 NPOL=IPAR(1) CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,2*NPOL V(J)=PAR(IPAR(3 + (4-1)) + (J-1)) END DO C RETURN END SUBROUTINE RHO SUBROUTINE RHOA(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(OUT):: A(:) REAL (KIND=R8), INTENT(IN):: LAMBDA,X(:) C C CALCULATE AND RETURN IN A THE VECTOR Z SUCH THAT C RHO(Z,LAMBDA,X) = 0 . C A(1)=LAMBDA ! INTENT(OUT) VARIABLE MUST BE DEFINED. RETURN END SUBROUTINE RHOA SUBROUTINE RHOJAC(A,LAMBDA,X,V,K) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: IPAR, PAR REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C RETURN IN THE VECTOR V THE KTH COLUMN OF THE JACOBIAN C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT THE POINT C (A, LAMBDA, X). C C THE FOLLOWING CODE IS SPECIFICALLY FOR THE POLYNOMIAL SYSTEM DRIVER C POLSYS1H , AND SHOULD BE USED VERBATIM WITH POLSYS1H . IF THE USER IS C CALLING FIXP?? OR STEP?? DIRECTLY, HE MUST SUPPLY APPROPRIATE C REPLACEMENT CODE HERE. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL,N2 NPOL=IPAR(1) N2=2*NPOL IF (K .EQ. 1) THEN C FORCE PREDICTED POINT TO HAVE LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,N2 V(J)=PAR(IPAR(3 + (6-1)) + (J-1)) END DO RETURN ELSE DO J=1,N2 V(J)=PAR(IPAR(3 + (5-1)) + (J-1) + N2*(K-2)) END DO ENDIF C RETURN END SUBROUTINE RHOJAC SUBROUTINE FJACS(X) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: QRSPARSE, ROWPOS, COLPOS REAL (KIND=R8), INTENT(IN):: X(:) C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C RETURN END SUBROUTINE FJACS SUBROUTINE RHOJS(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL, ONLY: QRSPARSE, ROWPOS, COLPOS REAL (KIND=R8), INTENT(IN):: A(:),LAMBDA,X(:) C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C RETURN END SUBROUTINE RHOJS SHAR_EOF fi # end of overwriting check if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << \SHAR_EOF > 'driver1.f' C MAIN PROGRAM TO TEST FIXPQF, FIXPNF, FIXPDF, STEPNX, AND ROOTNX. C BROWN'S FUNCTION, ZERO FINDING. C C THIS PROGRAM TESTS THE HOMPACK ROUTINES FIXPQF, FIXPNF, C FIXPDF, STEPNX, AND ROOTNX. C C THE OUTPUT FROM THIS ROUTINE SHOULD BE AS FOLLOWS, WITH THE C EXECUTION TIMES CORRESPONDING TO A DEC AXP 3000/600. C C TESTING FIXPQF C C LAMBDA = 1.00000000 FLAG = 1 6 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.106 ARCLEN = 2.693 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 C C C TESTING FIXPNF C C LAMBDA = 1.00000000 FLAG = 1 19 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.005 ARCLEN = 2.676 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 C C C TESTING FIXPDF C C LAMBDA = 1.00000000 FLAG = 1 71 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.016 ARCLEN = 2.712 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 C C C TESTING STEPNX AND ROOTNX C C LAMBDA = 1.00000000 FLAG = -1 80 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.020 ARCLEN = 2.711 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 C C PROGRAM TESTF USE REAL_PRECISION, ONLY : R8 USE HOMPACK90, ONLY : FIXPDF, FIXPNF, FIXPQF IMPLICIT NONE INTEGER, PARAMETER:: N=5, NDIMA=5 REAL (KIND=R8):: A(N),ANSAE,ANSRE,ARCAE,ARCRE, & ARCLEN,DTIME,SSPAR(8),Y(N+1) INTEGER:: IFLAG,II,J,NFE,NP1,TIMENEW(8),TIMEOLD(8),TRACE CHARACTER (LEN=6) NAME ! If using a subroutine library of the HOMPACK90 subroutines rather than ! the MODULE HOMPACK90 (as above), then the following INTERFACE ! statements are necessary. ! INTERFACE ! SUBROUTINE FIXPDF(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA, ! & NFE,ARCLEN) ! USE REAL_PRECISION ! INTEGER, INTENT(IN)::N,NDIMA,TRACE ! REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y ! INTEGER, INTENT(IN OUT)::IFLAG ! REAL (KIND=R8), INTENT(IN OUT)::ARCTOL,EPS ! INTEGER, INTENT(OUT)::NFE ! REAL (KIND=R8), INTENT(OUT)::ARCLEN ! END SUBROUTINE FIXPDF C ! SUBROUTINE FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, ! & SSPAR,NFE,ARCLEN,POLY_SWITCH) ! USE REAL_PRECISION ! INTEGER, INTENT(IN)::N,TRACE ! REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y ! INTEGER, INTENT(IN OUT)::IFLAG ! REAL (KIND=R8), INTENT(IN OUT)::ANSAE,ANSRE,ARCAE,ARCRE, ! & SSPAR(8) ! INTEGER, INTENT(OUT)::NFE ! REAL (KIND=R8), INTENT(OUT)::ARCLEN ! LOGICAL, INTENT(IN), OPTIONAL::POLY_SWITCH ! END SUBROUTINE FIXPNF C ! SUBROUTINE FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, ! & SSPAR,NFE,ARCLEN) ! USE REAL_PRECISION ! INTEGER, INTENT(IN)::N,TRACE ! REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y ! INTEGER, INTENT(IN OUT)::IFLAG ! REAL (KIND=R8), INTENT(IN OUT)::ANSAE,ANSRE,ARCAE,ARCRE, ! & SSPAR(4) ! INTEGER, INTENT(OUT)::NFE ! REAL (KIND=R8), INTENT(OUT)::ARCLEN ! END SUBROUTINE FIXPQF ! END INTERFACE C C TEST EACH OF THE THREE ALGORITHMS. C DO II=1,3 C C DEFINE ARGUMENTS FOR CALL TO HOMPACK PROCEDURE. C NP1=N+1 ARCRE=0.5D-4 ARCAE=0.5D-4 ANSRE=1.0D-10 ANSAE=1.0D-10 TRACE=0 SSPAR=0.0 IFLAG=-1 Y(2:NP1)=0.0 C C GET CURRENT DATE AND TIME. C CALL DATE_AND_TIME(VALUES=TIMEOLD) C C CALL TO HOMPACK ROUTINE. C IF (II .EQ. 1) THEN NAME='FIXPQF' CALL FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, & SSPAR,NFE,ARCLEN) ELSE IF (II .EQ. 2) THEN NAME='FIXPNF' CALL FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, & SSPAR,NFE,ARCLEN) ELSE NAME='FIXPDF' CALL FIXPDF(N,Y,IFLAG,ARCRE,ANSRE,TRACE,A,NDIMA,NFE,ARCLEN) END IF C C CALCULATE EXECUTION TIME. C CALL DATE_AND_TIME(VALUES=TIMENEW) IF (TIMENEW(8) .LT. TIMEOLD(8)) THEN TIMENEW(8)=TIMENEW(8)+1000 TIMENEW(7)=TIMENEW(7)-1 ENDIF IF (TIMENEW(7) .LT. TIMEOLD(7)) THEN TIMENEW(7)=TIMENEW(7)+60 TIMENEW(6)=TIMENEW(6)-1 ENDIF IF (TIMENEW(6) .LT. TIMEOLD(6)) THEN TIMENEW(6)=TIMENEW(6)+60 TIMENEW(5)=TIMENEW(5)-1 ENDIF IF (TIMENEW(5) .LT. TIMEOLD(5)) TIMENEW(5)=TIMENEW(5)+24 DTIME=DOT_PRODUCT(TIMENEW(5:8)-TIMEOLD(5:8), & (/3600000,60000,1000,1/) )/1000.0 C WRITE (6,45) NAME 45 FORMAT (//,7X,'TESTING',1X,A6) WRITE (6,50) Y(1),IFLAG,NFE,DTIME,ARCLEN,(Y(J),J=2,NP1) 50 FORMAT(/' LAMBDA =',F11.8,' FLAG =',I2,I8,' JACOBIAN ', & 'EVALUATIONS',/,1X,'EXECUTION TIME(SECS) =',F10.3,4X, & 'ARCLEN =',F10.3/(1X,4ES16.8)) END DO C C TEST REVERSE CALL SUBROUTINES STEPNX AND ROOTNX ON THE SAME C PROBLEM. C CALL MAINX STOP END PROGRAM TESTF ! ! SAMPLE USER WRITTEN HOMOTOPY SUBROUTINES FOR TESTING FIXP*F. ! SUBROUTINE F(X,V) C******************************************************************** C C SUBROUTINE F(X,V) -- EVALUATES BROWN'S FUNCTION AT THE POINT C X, AND RETURNS THE VALUE IN V. C C******************************************************************** USE REAL_PRECISION, ONLY : R8 IMPLICIT NONE REAL (KIND=R8), INTENT(IN):: X(:) REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER:: N N=SIZE(X) V(1)=PRODUCT(X) - 1.0 V(2:N)=SUM(X) - (N+1) + X(2:N) RETURN END SUBROUTINE F SUBROUTINE FJAC(X,V,K) C******************************************************************** C C SUBROUTINE FJAC(X,V,K) -- EVALUATES THE K-TH COLUMN OF C THE JACOBIAN MATRIX FOR BROWN'S FUNCTION EVALUATED AT C THE POINT X, RETURNING THE VALUE IN V. C C******************************************************************** USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(IN):: X(:) REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K INTEGER:: J,N REAL (KIND=R8):: PROD C N=SIZE(X) PROD=1.0 DO J=1,K-1 PROD=PROD*X(J) END DO DO J=K+1,N PROD=PROD*X(J) END DO V(1)=PROD V(2:N)=1.0 IF (K .GT. 1) V(K)=V(K)+1.0 RETURN END SUBROUTINE FJAC C ********************************************************************** C C THE REST OF THESE SUBROUTINES ARE NOT USED BY PROGRAM TESTF, AND ARE C INCLUDED HERE SIMPLY FOR COMPLETENESS AND AS TEMPLATES FOR THEIR USE. C SUBROUTINE RHO(A,LAMBDA,X,V) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) C C EVALUATE RHO(A,LAMBDA,X) AND RETURN IN THE VECTOR V . C C THE FOLLOWING CODE IS SPECIFICALLY FOR THE POLYNOMIAL SYSTEM DRIVER C POLSYS1H , AND SHOULD BE USED VERBATIM WITH POLSYS1H . IF THE USER IS C CALLING FIXP?? OR STEP?? DIRECTLY, HE MUST SUPPLY APPROPRIATE C REPLACEMENT CODE HERE. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL C FORCE PREDICTED POINT TO HAVE LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 NPOL=IPAR(1) ! CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,2*NPOL V(J)=PAR(IPAR(3 + (4-1)) + (J-1)) END DO C RETURN END SUBROUTINE RHO SUBROUTINE RHOA(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 REAL (KIND=R8), INTENT(OUT):: A(:) REAL (KIND=R8), INTENT(IN):: LAMBDA,X(:) C C CALCULATE AND RETURN IN A THE VECTOR Z SUCH THAT C RHO(Z,LAMBDA,X) = 0 . C A(1)=LAMBDA ! INTENT(OUT) VARIABLE MUST BE DEFINED. RETURN END SUBROUTINE RHOA SUBROUTINE RHOJAC(A,LAMBDA,X,V,K) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL REAL (KIND=R8), INTENT(IN):: A(:),X(:) REAL (KIND=R8), INTENT(IN OUT):: LAMBDA REAL (KIND=R8), INTENT(OUT):: V(:) INTEGER, INTENT(IN):: K C C RETURN IN THE VECTOR V THE KTH COLUMN OF THE JACOBIAN C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT THE POINT C (A, LAMBDA, X). C C THE FOLLOWING CODE IS SPECIFICALLY FOR THE POLYNOMIAL SYSTEM DRIVER C POLSYS1H , AND SHOULD BE USED VERBATIM WITH POLSYS1H . IF THE USER IS C CALLING FIXP?? OR STEP?? DIRECTLY, HE MUST SUPPLY APPROPRIATE C REPLACEMENT CODE HERE. INTERFACE SUBROUTINE HFUNP(N,A,LAMBDA,X) USE REAL_PRECISION INTEGER, INTENT(IN):: N REAL (KIND=R8), INTENT(IN):: A(2*N),LAMBDA,X(2*N) END SUBROUTINE HFUNP END INTERFACE INTEGER:: J,NPOL,N2 NPOL=IPAR(1) N2=2*NPOL IF (K .EQ. 1) THEN C FORCE PREDICTED POINT TO HAVE LAMBDA .GE. 0 . IF (LAMBDA .LT. 0.0) LAMBDA=0.0 ! CALL HFUNP(NPOL,A,LAMBDA,X) DO J=1,N2 V(J)=PAR(IPAR(3 + (6-1)) + (J-1)) END DO RETURN ELSE DO J=1,N2 V(J)=PAR(IPAR(3 + (5-1)) + (J-1) + N2*(K-2)) END DO ENDIF C RETURN END SUBROUTINE RHOJAC SUBROUTINE FJACS(X) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL REAL (KIND=R8), INTENT(IN):: X(:) C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C RETURN END SUBROUTINE FJACS SUBROUTINE RHOJS(A,LAMBDA,X) USE REAL_PRECISION, ONLY : R8 USE HOMPACK90_GLOBAL REAL (KIND=R8), INTENT(IN):: A(:),LAMBDA,X(:) C C If MODE = 1, C evaluate the N x N symmetric Jacobian matrix of F(X) at X, and return C the result in packed skyline storage format in QRSPARSE. LENQR is the C length of QRSPARSE, and ROWPOS contains the indices of the diagonal C elements of the Jacobian matrix within QRSPARSE. ROWPOS(N+1) and C ROWPOS(N+2) are set by subroutine FODEDS. The allocatable array COLPOS C is not used by this storage format. C C If MODE = 2, C evaluate the N x N Jacobian matrix of F(X) at X, and return the result C in sparse row storage format in QRSPARSE. LENQR is the length of C QRSPARSE, ROWPOS contains the indices of where each row begins within C QRSPARSE, and COLPOS (of length LENQR) contains the column indices of C the corresponding elements in QRSPARSE. Even if zero, the diagonal C elements of the Jacobian matrix must be stored in QRSPARSE. C RETURN END SUBROUTINE RHOJS C ********************************************************************** C C SUBROUTINE TO TEST THE REVERSE CALL SUBROUTINES STEPNX AND C ROOTNX. THE TEST PROBLEM IS BROWN'S FUNCTION, ZERO FINDING. C THE OUTPUT IS SIMILAR TO THAT FROM THE TEST OF FIXPNF, EXCEPT WITH C MORE JACOBIAN EVALUATIONS SINCE THE UNDEFINED FUNCTION OPTION OF C STEPNX IS USED TO FORCE SMALLER STEPS. C SUBROUTINE MAINX USE REAL_PRECISION, ONLY : R8 USE HOMOTOPY IMPLICIT NONE INTEGER, PARAMETER:: N=5, NDIMA=5 REAL (KIND=R8):: A(NDIMA),ABSERR,ALPHA(3*N+3), & ANSAE,ANSRE,ARCAE,ARCRE,ARCLEN,DTIME,GOFW,H,HOLD, & QR(N,N+2),RELERR,RHOLEN,S,SSPAR(8),TZ(N+1),W(N+1), & WP(N+1),Y(N+1),YOLD(N+1),YOLDS(N+1),YP(N+1),YPOLD(N+1) INTEGER:: IFLAG,ITER=0,J,NFE,NFEC=0,NP1,PIVOT(N+1), & TIMENEW(8),TIMEOLD(8),TRACE LOGICAL:: CRASH, START C INTERFACE SUBROUTINE ROOTNX(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD, & YPOLD,A,GOFW,TZ,W,WP) USE HOMOTOPY USE REAL_PRECISION INTEGER, INTENT(IN):: N INTEGER, INTENT(IN OUT):: NFE,IFLAG REAL (KIND=R8), INTENT(IN):: RELERR,ABSERR REAL (KIND=R8), DIMENSION(:), INTENT(IN):: A REAL (KIND=R8), INTENT(IN OUT):: GOFW REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT):: Y,YP,YOLD,YPOLD, & TZ,W,WP END SUBROUTINE ROOTNX SUBROUTINE STEPNX(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR, & ABSERR,S,Y,YP,YOLD,YPOLD,A,TZ,W,WP,RHOLEN,SSPAR) USE HOMOTOPY USE REAL_PRECISION INTEGER, INTENT(IN):: N INTEGER, INTENT(IN OUT):: NFE,IFLAG LOGICAL, INTENT(IN OUT):: START,CRASH REAL (KIND=R8), INTENT(IN OUT):: HOLD,H,RELERR,ABSERR,S,RHOLEN, & SSPAR(8) REAL (KIND=R8), DIMENSION(:), INTENT(IN):: A REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT):: Y,YP,YOLD,YPOLD, & TZ,W,WP REAL (KIND=R8), DIMENSION(:), ALLOCATABLE, SAVE:: Z0,Z1 END SUBROUTINE STEPNX SUBROUTINE TANGNF(RHOLEN,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT, & NFE,N,IFLAG) USE REAL_PRECISION REAL (KIND=R8):: RHOLEN INTEGER:: IFLAG,N,NFE REAL (KIND=R8):: A(:),Y(:),YP(N+1),YPOLD(N+1) REAL (KIND=R8):: ALPHA(3*N+3),QR(N,N+2),TZ(N+1) INTEGER:: PIVOT(N+1) END SUBROUTINE TANGNF END INTERFACE C C DEFINE ARGUMENTS FOR CALL TO HOMPACK PROCEDURE. C NP1=N+1 NFE=0 ARCRE=0.5D-4 ARCAE=0.5D-4 ANSRE=1.0D-10 ANSAE=1.0D-10 ABSERR=ARCAE; RELERR=ARCRE TRACE=0 SSPAR=0.0 IFLAG=-1 A=0.0 Y(1:NP1)=0.0 YP(1)=1.0; YP(2:NP1)=0.0 YOLD=Y; YPOLD=YP START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 S=0.0 C C GET CURRENT DATE AND TIME. C CALL DATE_AND_TIME(VALUES=TIMEOLD) C C TRACK CURVE TILL LAMBDA > 1.0 . C TRACK: DO WHILE (Y(1) < 1.0_R8) CALL STEPNX (N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR, & ABSERR,S,Y,YP,YOLD,YPOLD,A,TZ,W,WP,RHOLEN,SSPAR) IF (CRASH) CYCLE TRACK SELECT CASE (IFLAG) CASE (-2:0) IF (TRACE .GT. 0) THEN ITER=ITER+1 WRITE (TRACE,11) ITER,NFE,S,Y(1),Y(2:NP1) 11 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =', & F9.4,3X,'LAMBDA =',F7.4,5X,'X VECTOR:'/(1X,6ES12.4)) ENDIF CYCLE TRACK CASE (-12:-10) ! TANGENT VECTOR IF (H > .1_R8) THEN IFLAG = IFLAG - 100 CYCLE TRACK END IF RHOLEN=0.0 CALL TANGNF(RHOLEN,W,WP,YPOLD,A,QR,ALPHA,TZ,PIVOT, & NFEC,N,IFLAG) CASE (-32:-20) ! TANGENT VECTOR AND NEWTON STEP IF (H > .1_R8) THEN IFLAG = IFLAG - 100 CYCLE TRACK END IF RHOLEN=-1.0 CALL TANGNF(RHOLEN,W,WP,YPOLD,A,QR,ALPHA,TZ,PIVOT, & NFEC,N,IFLAG) CASE (4,6,7) WRITE (6,13) IFLAG 13 FORMAT(/' FATAL ERROR OCCURRED DURING TRACKING WITH', & ' FLAG =',I2,//) STOP END SELECT END DO TRACK C C CLEAN UP WORKING STORAGE. IFLAG=IFLAG - 40 CALL STEPNX (N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR, & ABSERR,S,Y,YP,YOLD,YPOLD,A,TZ,W,WP,RHOLEN,SSPAR) C SAVE YOLD FOR ARC LENGTH CALCULATION LATER. YOLDS = YOLD C C FIND POINT ON HOMOTOPY ZERO CURVE SATISFYING C G(Y(S)) = LAMBDA(S) - 1 = 0 . C ABSERR=ANSAE RELERR=ANSRE END_GAME: DO CALL ROOTNX(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD, & YPOLD,A,GOFW,TZ,W,WP) SELECT CASE (IFLAG) CASE (-42:-10) ! G(W) GOFW = W(1) - 1.0 CASE (-52:-50) ! TANGENT VECTOR AND NEWTON STEP RHOLEN=-1.0 CALL TANGNF(RHOLEN,W,WP,YPOLD,A,QR,ALPHA,TZ,PIVOT, & NFEC,N,IFLAG) CASE (-2:0, 4, 6, 7) EXIT END_GAME END SELECT END DO END_GAME C CALCULATE FINAL ARC LENGTH. W = Y - YOLDS ARCLEN = S - HOLD + SQRT(DOT_PRODUCT(W,W)) C C CALCULATE EXECUTION TIME. C CALL DATE_AND_TIME(VALUES=TIMENEW) IF (TIMENEW(8) .LT. TIMEOLD(8)) THEN TIMENEW(8)=TIMENEW(8)+1000 TIMENEW(7)=TIMENEW(7)-1 ENDIF IF (TIMENEW(7) .LT. TIMEOLD(7)) THEN TIMENEW(7)=TIMENEW(7)+60 TIMENEW(6)=TIMENEW(6)-1 ENDIF IF (TIMENEW(6) .LT. TIMEOLD(6)) THEN TIMENEW(6)=TIMENEW(6)+60 TIMENEW(5)=TIMENEW(5)-1 ENDIF IF (TIMENEW(5) .LT. TIMEOLD(5)) TIMENEW(5)=TIMENEW(5)+24 DTIME=DOT_PRODUCT(TIMENEW(5:8)-TIMEOLD(5:8), & (/3600000,60000,1000,1/) )/1000.0 C WRITE (6,45) 45 FORMAT (//,7X,'TESTING STEPNX AND ROOTNX') WRITE (6,50) Y(1),IFLAG,NFE,DTIME,ARCLEN,(Y(J),J=2,NP1) 50 FORMAT(/' LAMBDA =',F11.8,' FLAG =',I3,I8,' JACOBIAN ', & 'EVALUATIONS',/,1X,'EXECUTION TIME(SECS) =',F10.3,4X, & 'ARCLEN =',F10.3/(1X,4ES16.8)) RETURN END SUBROUTINE MAINX SHAR_EOF fi # end of overwriting check if test -f 'RES3' then echo shar: will not over-write existing file "'RES3'" else cat << \SHAR_EOF > 'RES3' TESTING FIXPQS WITH STORAGE MODE = 1 LAMBDA = 1.00000000 FLAG = 1 33 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.108 ARC LENGTH = 1.274 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 TESTING FIXPNS WITH STORAGE MODE = 1 LAMBDA = 1.00000000 FLAG = 1 20 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.012 ARC LENGTH = 1.275 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 TESTING FIXPDS WITH STORAGE MODE = 1 LAMBDA = 1.00000000 FLAG = 1 70 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.022 ARC LENGTH = 1.281 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 TESTING FIXPQS WITH STORAGE MODE = 2 LAMBDA = 1.00000000 FLAG = 1 33 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.015 ARC LENGTH = 1.274 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 TESTING FIXPNS WITH STORAGE MODE = 2 LAMBDA = 1.00000000 FLAG = 1 20 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.011 ARC LENGTH = 1.275 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 TESTING FIXPDS WITH STORAGE MODE = 2 LAMBDA = 1.00000000 FLAG = 1 70 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.020 ARC LENGTH = 1.281 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 SHAR_EOF fi # end of overwriting check if test -f 'RES2' then echo shar: will not over-write existing file "'RES2'" else cat << \SHAR_EOF > 'RES2' POLSYS1H TEST ROUTINE 7/7/95 TWO QUADRICS, NO SOLUTIONS AT INFINITY, TWO REAL SOLUTIONS. IF IFLGHM=1, HOMOGENEOUS; IF IFLGHM=0, INHOMOGENEOUS; IFLGHM= 1 IF IFLGSC=1, SCLGNP USED; IF IFLGSC=0, NO SCALING; IFLGSC= 1 TOTDG= 4 MAXT= 6 EPSBIG, EPSSML = 1.00000000000000E-04 1.00000000000000E-14 SSPAR(5) = 1.00000000000000E+00 NUMBER OF EQUATIONS = 2 NUMBER OF RECALLS WHEN IFLAG=3: 10 ****** COEFFICIENT TABLEAU ****** NUMT( 1) = 6 KDEG( 1, 1, 1) = 2 KDEG( 1, 2, 1) = 0 COEF( 1, 1) = -9.80000000000000E-04 KDEG( 1, 1, 2) = 0 KDEG( 1, 2, 2) = 2 COEF( 1, 2) = 9.78000000000000E+05 KDEG( 1, 1, 3) = 1 KDEG( 1, 2, 3) = 1 COEF( 1, 3) = -9.80000000000000E+00 KDEG( 1, 1, 4) = 1 KDEG( 1, 2, 4) = 0 COEF( 1, 4) = -2.35000000000000E+02 KDEG( 1, 1, 5) = 0 KDEG( 1, 2, 5) = 1 COEF( 1, 5) = 8.89000000000000E+04 KDEG( 1, 1, 6) = 0 KDEG( 1, 2, 6) = 0 COEF( 1, 6) = -1.00000000000000E+00 NUMT( 2) = 6 KDEG( 2, 1, 1) = 2 KDEG( 2, 2, 1) = 0 COEF( 2, 1) = -1.00000000000000E-02 KDEG( 2, 1, 2) = 0 KDEG( 2, 2, 2) = 2 COEF( 2, 2) = -9.84000000000000E-01 KDEG( 2, 1, 3) = 1 KDEG( 2, 2, 3) = 1 COEF( 2, 3) = -2.97000000000000E+01 KDEG( 2, 1, 4) = 1 KDEG( 2, 2, 4) = 0 COEF( 2, 4) = 9.87000000000000E-03 KDEG( 2, 1, 5) = 0 KDEG( 2, 2, 5) = 1 COEF( 2, 5) = -1.24000000000000E-01 KDEG( 2, 1, 6) = 0 KDEG( 2, 2, 6) = 0 COEF( 2, 6) = -2.50000000000000E-01 IFLG1 = 11 PATH NUMBER = 1 FINAL VALUES FOR PATH ARCLEN = 1.00553319056290E+01 NFE = 53 IFLG2 = 1 REAL, FINITE SOLUTION LAMBDA = 1.00000000000000E+00 X( 1) = 2.34233851959128E+03 1.49779467841657E-11 X( 2) = -7.88344824094141E-01 -5.56115428011477E-15 X( 3) = -9.49359459408655E-03 -1.06447550900257E-03 PATH NUMBER = 2 FINAL VALUES FOR PATH ARCLEN = 1.72112928605711E+00 NFE = 37 IFLG2 = 1 COMPLEX, FINITE SOLUTION LAMBDA = 1.00000000000001E+00 X( 1) = 1.61478579234419E-02 1.68496955498881E+00 X( 2) = 2.67994739614476E-04 4.42802993973660E-03 X( 3) = -3.81948972942403E-01 3.72068943457283E-01 PATH NUMBER = 3 FINAL VALUES FOR PATH ARCLEN = 2.02329527936743E+00 NFE = 35 IFLG2 = 1 COMPLEX, FINITE SOLUTION LAMBDA = 1.00000000000000E+00 X( 1) = 1.61478579234355E-02 -1.68496955498881E+00 X( 2) = 2.67994739614460E-04 -4.42802993973661E-03 X( 3) = -3.29370493847660E-01 5.56619775523013E-01 PATH NUMBER = 4 FINAL VALUES FOR PATH ARCLEN = 4.16326615695899E+00 NFE = 46 IFLG2 = 1 REAL, FINITE SOLUTION LAMBDA = 1.00000000000000E+00 X( 1) = 9.08921229615388E-02 1.90036587651500E-16 X( 2) = -9.11497098197499E-02 4.71849760398007E-18 X( 3) = -5.73673395727962E-02 1.36243663709219E-01 TOTAL NFE OVER ALL PATHS = 171 SHAR_EOF fi # end of overwriting check if test -f 'RES1' then echo shar: will not over-write existing file "'RES1'" else cat << \SHAR_EOF > 'RES1' TESTING FIXPQF LAMBDA = 1.00000000 FLAG = 1 6 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.117 ARCLEN = 2.693 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 TESTING FIXPNF LAMBDA = 1.00000000 FLAG = 1 19 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.004 ARCLEN = 2.676 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 TESTING FIXPDF LAMBDA = 1.00000000 FLAG = 1 71 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.011 ARCLEN = 2.712 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 TESTING STEPNX AND ROOTNX LAMBDA = 1.00000000 FLAG = -1 80 JACOBIAN EVALUATIONS EXECUTION TIME(SECS) = 0.011 ARCLEN = 2.711 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 SHAR_EOF fi # end of overwriting check if test -f 'DATA2' then echo shar: will not over-write existing file "'DATA2'" else cat << \SHAR_EOF > 'DATA2' ' TWO QUADRICS, NO SOLUTIONS AT INFINITY, TWO REAL SOLUTIONS.' &PROBLEM IFLGHM = 1 IFLGSC = 1 TOTDG = 4 MAXT = 6 EPSBIG = 1.D-04 EPSSML = 1.D-14 SSPAR(5) = 1.D+00 NUMRR = 10 N = 2 / 00006 NUMTRM(1) 00002 DEG(1,1,1) 00000 DEG(1,2,1) -.00098D+00 00000 DEG(1,1,2) 00002 DEG(1,2,2) 978000.D+00 00001 DEG(1,1,3) 00001 DEG(1,2,3) -9.8D+00 00001 DEG(1,1,4) 00000 DEG(1,2,4) -235.0D+00 00000 DEG(1,1,5) 00001 DEG(1,2,5) 88900.0D+00 00000 DEG(1,1,6) 00000 DEG(1,2,6) -1.000D+00 00006 NUMTRM(2) 00002 DEG(2,1,1) 00000 DEG(2,2,1) -.0100D+00 00000 DEG(2,1,2) 00002 DEG(2,2,2) -.9840D+00 00001 DEG(2,1,3) 00001 DEG(2,2,3) -29.70D+00 00001 DEG(2,1,4) 00000 DEG(2,2,4) .00987D+00 00000 DEG(2,1,5) 00001 DEG(2,2,5) -.1240D+00 00000 DEG(2,1,6) 00000 DEG(2,2,6) -.2500D+00 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Info' then mkdir 'Info' fi cd 'Info' if test -f 'depend' then echo shar: will not over-write existing file "'depend'" else cat << \SHAR_EOF > 'depend' blas1/dcopy.f blas1/ddot.f blas3/dgemm.f blas2/dgemv.f blas2/dger.f blas1/dnrm2.f blas1/dscal.f blas1/dswap.f blas2/dtpmv.f blas2/dtpsv.f blas3/dtrmm.f blas2/dtrmv.f blas2/dtrsv.f blas1/idamax.f lapack/dgeqpf.f lapack/dgeqr2.f lapack/dgeqrf.f lapack/dlaic1.f lapack/dlamch.f lapack/dlamc1.f lapack/dlamc2.f lapack/dlamc3.f lapack/dlamc4.f lapack/dlamc5.f lapack/dlapy2.f lapack/dlarf.f lapack/dlarfb.f lapack/dlarfg.f lapack/dlarft.f lapack/dorg2r.f lapack/dorgqr.f lapack/dorm2r.f lapack/dormqr.f lapack/ilaenv.f lapack/lsame.f lapack/xerbla.f SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'lapack.f' then echo shar: will not over-write existing file "'lapack.f'" else cat << \SHAR_EOF > 'lapack.f' SUBROUTINE DGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO ) * * -- LAPACK test routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER JPVT( * ) DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * DGEQPF computes a QR factorization with column pivoting of a * real M-by-N matrix A: A*P = Q*R. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0 * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, the upper triangle of the array contains the * min(M,N)-by-N upper triangular matrix R; the elements * below the diagonal, together with the array TAU, * represent the orthogonal matrix Q as a product of * min(m,n) elementary reflectors. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * JPVT (input/output) INTEGER array, dimension (N) * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted * to the front of A*P (a leading column); if JPVT(i) = 0, * the i-th column of A is a free column. * On exit, if JPVT(i) = k, then the i-th column of A*P * was the k-th column of A. * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors. * * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(n) * * Each H(i) has the form * * H = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). * * The matrix P is represented in jpvt as follows: If * jpvt(j) = i * then the jth column of P is the ith canonical unit vector. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, ITEMP, J, MA, MN, PVT DOUBLE PRECISION AII, TEMP, TEMP2 * .. * .. External Subroutines .. EXTERNAL DGEQR2, DLARF, DLARFG, DORM2R, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. External Functions .. INTEGER IDAMAX DOUBLE PRECISION DNRM2 EXTERNAL IDAMAX, DNRM2 * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGEQPF', -INFO ) RETURN END IF * MN = MIN( M, N ) * * Move initial columns up front * ITEMP = 1 DO 10 I = 1, N IF( JPVT( I ).NE.0 ) THEN IF( I.NE.ITEMP ) THEN CALL DSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 ) JPVT( I ) = JPVT( ITEMP ) JPVT( ITEMP ) = I ELSE JPVT( I ) = I END IF ITEMP = ITEMP + 1 ELSE JPVT( I ) = I END IF 10 CONTINUE ITEMP = ITEMP - 1 * * Compute the QR factorization and update remaining columns * IF( ITEMP.GT.0 ) THEN MA = MIN( ITEMP, M ) CALL DGEQR2( M, MA, A, LDA, TAU, WORK, INFO ) IF( MA.LT.N ) THEN CALL DORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU, $ A( 1, MA+1 ), LDA, WORK, INFO ) END IF END IF * IF( ITEMP.LT.MN ) THEN * * Initialize partial column norms. The first n elements of * work store the exact column norms. * DO 20 I = ITEMP + 1, N WORK( I ) = DNRM2( M-ITEMP, A( ITEMP+1, I ), 1 ) WORK( N+I ) = WORK( I ) 20 CONTINUE * * Compute factorization * DO 40 I = ITEMP + 1, MN * * Determine ith pivot column and swap if necessary * PVT = ( I-1 ) + IDAMAX( N-I+1, WORK( I ), 1 ) * IF( PVT.NE.I ) THEN CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) ITEMP = JPVT( PVT ) JPVT( PVT ) = JPVT( I ) JPVT( I ) = ITEMP WORK( PVT ) = WORK( I ) WORK( N+PVT ) = WORK( N+I ) END IF * * Generate elementary reflector H(i) * IF( I.LT.M ) THEN CALL DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) ) ELSE CALL DLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) ) END IF * IF( I.LT.N ) THEN * * Apply H(i) to A(i:m,i+1:n) from the left * AII = A( I, I ) A( I, I ) = ONE CALL DLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ), $ A( I, I+1 ), LDA, WORK( 2*N+1 ) ) A( I, I ) = AII END IF * * Update partial column norms * DO 30 J = I + 1, N IF( WORK( J ).NE.ZERO ) THEN TEMP = ONE - ( ABS( A( I, J ) ) / WORK( J ) )**2 TEMP = MAX( TEMP, ZERO ) TEMP2 = ONE + 0.05D0*TEMP* $ ( WORK( J ) / WORK( N+J ) )**2 IF( TEMP2.EQ.ONE ) THEN IF( M-I.GT.0 ) THEN WORK( J ) = DNRM2( M-I, A( I+1, J ), 1 ) WORK( N+J ) = WORK( J ) ELSE WORK( J ) = ZERO WORK( N+J ) = ZERO END IF ELSE WORK( J ) = WORK( J )*SQRT( TEMP ) END IF END IF 30 CONTINUE * 40 CONTINUE END IF RETURN * * End of DGEQPF * END SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * * -- LAPACK routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * DGEQR2 computes a QR factorization of a real m by n matrix A: * A = Q * R. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the elements on and above the diagonal of the array * contain the min(m,n) by n upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(k), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), * and tau in TAU(i). * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, K DOUBLE PRECISION AII * .. * .. External Subroutines .. EXTERNAL DLARF, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGEQR2', -INFO ) RETURN END IF * K = MIN( M, N ) * DO 10 I = 1, K * * Generate elementary reflector H(i) to annihilate A(i+1:m,i) * CALL DLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, $ TAU( I ) ) IF( I.LT.N ) THEN * * Apply H(i) to A(i:m,i+1:n) from the left * AII = A( I, I ) A( I, I ) = ONE CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), $ A( I, I+1 ), LDA, WORK ) A( I, I ) = AII END IF 10 CONTINUE RETURN * * End of DGEQR2 * END SUBROUTINE DGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * * -- LAPACK routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, M, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( LWORK ) * .. * * Purpose * ======= * * DGEQRF computes a QR factorization of a real M-by-N matrix A: * A = Q * R. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, the elements on and above the diagonal of the array * contain the min(M,N)-by-N upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of min(m,n) elementary reflectors (see Further * Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * For optimum performance LWORK >= N*NB, where NB is * the optimal blocksize. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(k), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), * and tau in TAU(i). * * ===================================================================== * * .. Local Scalars .. INTEGER I, IB, IINFO, IWS, K, LDWORK, NB, NBMIN, NX * .. * .. External Subroutines .. EXTERNAL DGEQR2, DLARFB, DLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGEQRF', -INFO ) RETURN END IF * * Quick return if possible * K = MIN( M, N ) IF( K.EQ.0 ) THEN WORK( 1 ) = 1 RETURN END IF * * Determine the block size. * NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) NBMIN = 2 NX = 0 IWS = N IF( NB.GT.1 .AND. NB.LT.K ) THEN * * Determine when to cross over from blocked to unblocked code. * NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) ) IF( NX.LT.K ) THEN * * Determine if workspace is large enough for blocked code. * LDWORK = N IWS = LDWORK*NB IF( LWORK.LT.IWS ) THEN * * Not enough workspace to use optimal NB: reduce NB and * determine the minimum value of NB. * NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1, $ -1 ) ) END IF END IF END IF * IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN * * Use blocked code initially * DO 10 I = 1, K - NX, NB IB = MIN( K-I+1, NB ) * * Compute the QR factorization of the current block * A(i:m,i:i+ib-1) * CALL DGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, $ IINFO ) IF( I+IB.LE.N ) THEN * * Form the triangular factor of the block reflector * H = H(i) H(i+1) . . . H(i+ib-1) * CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB, $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) * * Apply H' to A(i:m,i+ib:n) from the left * CALL DLARFB( 'Left', 'Transpose', 'Forward', $ 'Columnwise', M-I+1, N-I-IB+1, IB, $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), $ LDA, WORK( IB+1 ), LDWORK ) END IF 10 CONTINUE ELSE I = 1 END IF * * Use unblocked code to factor the last or only block. * IF( I.LE.K ) $ CALL DGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, $ IINFO ) * WORK( 1 ) = IWS RETURN * * End of DGEQRF * END SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER J, JOB DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR * .. * .. Array Arguments .. DOUBLE PRECISION W( J ), X( J ) * .. * * Purpose * ======= * * DLAIC1 applies one step of incremental condition estimation in * its simplest version: * * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j * lower triangular matrix L, such that * twonorm(L*x) = sest * Then DLAIC1 computes sestpr, s, c such that * the vector * [ s*x ] * xhat = [ c ] * is an approximate singular vector of * [ L 0 ] * Lhat = [ w' gamma ] * in the sense that * twonorm(Lhat*xhat) = sestpr. * * Depending on JOB, an estimate for the largest or smallest singular * value is computed. * * Note that [s c]' and sestpr**2 is an eigenpair of the system * * diag(sest*sest, 0) + [alpha gamma] * [ alpha ] * [ gamma ] * * where alpha = x'*w. * * Arguments * ========= * * JOB (input) INTEGER * = 1: an estimate for the largest singular value is computed. * = 2: an estimate for the smallest singular value is computed. * * J (input) INTEGER * Length of X and W * * X (input) DOUBLE PRECISION array, dimension (J) * The j-vector x. * * SEST (input) DOUBLE PRECISION * Estimated singular value of j by j matrix L * * W (input) DOUBLE PRECISION array, dimension (J) * The j-vector w. * * GAMMA (input) DOUBLE PRECISION * The diagonal element gamma. * * SEDTPR (output) DOUBLE PRECISION * Estimated singular value of (j+1) by (j+1) matrix Lhat. * * S (output) DOUBLE PRECISION * Sine needed in forming xhat. * * C (output) DOUBLE PRECISION * Cosine needed in forming xhat. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) DOUBLE PRECISION HALF, FOUR PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS, $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SIGN, SQRT * .. * .. External Functions .. DOUBLE PRECISION DDOT, DLAMCH EXTERNAL DDOT, DLAMCH * .. * .. Executable Statements .. * EPS = DLAMCH( 'Epsilon' ) ALPHA = DDOT( J, X, 1, W, 1 ) * ABSALP = ABS( ALPHA ) ABSGAM = ABS( GAMMA ) ABSEST = ABS( SEST ) * IF( JOB.EQ.1 ) THEN * * Estimating largest singular value * * special cases * IF( SEST.EQ.ZERO ) THEN S1 = MAX( ABSGAM, ABSALP ) IF( S1.EQ.ZERO ) THEN S = ZERO C = ONE SESTPR = ZERO ELSE S = ALPHA / S1 C = GAMMA / S1 TMP = SQRT( S*S+C*C ) S = S / TMP C = C / TMP SESTPR = S1*TMP END IF RETURN ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN S = ONE C = ZERO TMP = MAX( ABSEST, ABSALP ) S1 = ABSEST / TMP S2 = ABSALP / TMP SESTPR = TMP*SQRT( S1*S1+S2*S2 ) RETURN ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN S1 = ABSGAM S2 = ABSEST IF( S1.LE.S2 ) THEN S = ONE C = ZERO SESTPR = S2 ELSE S = ZERO C = ONE SESTPR = S1 END IF RETURN ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN S1 = ABSGAM S2 = ABSALP IF( S1.LE.S2 ) THEN TMP = S1 / S2 S = SQRT( ONE+TMP*TMP ) SESTPR = S2*S C = ( GAMMA / S2 ) / S S = SIGN( ONE, ALPHA ) / S ELSE TMP = S2 / S1 C = SQRT( ONE+TMP*TMP ) SESTPR = S1*C S = ( ALPHA / S1 ) / C C = SIGN( ONE, GAMMA ) / C END IF RETURN ELSE * * normal case * ZETA1 = ALPHA / ABSEST ZETA2 = GAMMA / ABSEST * B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF C = ZETA1*ZETA1 IF( B.GT.ZERO ) THEN T = C / ( B+SQRT( B*B+C ) ) ELSE T = SQRT( B*B+C ) - B END IF * SINE = -ZETA1 / T COSINE = -ZETA2 / ( ONE+T ) TMP = SQRT( SINE*SINE+COSINE*COSINE ) S = SINE / TMP C = COSINE / TMP SESTPR = SQRT( T+ONE )*ABSEST RETURN END IF * ELSE IF( JOB.EQ.2 ) THEN * * Estimating smallest singular value * * special cases * IF( SEST.EQ.ZERO ) THEN SESTPR = ZERO IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN SINE = ONE COSINE = ZERO ELSE SINE = -GAMMA COSINE = ALPHA END IF S1 = MAX( ABS( SINE ), ABS( COSINE ) ) S = SINE / S1 C = COSINE / S1 TMP = SQRT( S*S+C*C ) S = S / TMP C = C / TMP RETURN ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN S = ZERO C = ONE SESTPR = ABSGAM RETURN ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN S1 = ABSGAM S2 = ABSEST IF( S1.LE.S2 ) THEN S = ZERO C = ONE SESTPR = S1 ELSE S = ONE C = ZERO SESTPR = S2 END IF RETURN ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN S1 = ABSGAM S2 = ABSALP IF( S1.LE.S2 ) THEN TMP = S1 / S2 C = SQRT( ONE+TMP*TMP ) SESTPR = ABSEST*( TMP / C ) S = -( GAMMA / S2 ) / C C = SIGN( ONE, ALPHA ) / C ELSE TMP = S2 / S1 S = SQRT( ONE+TMP*TMP ) SESTPR = ABSEST / S C = ( ALPHA / S1 ) / S S = -SIGN( ONE, GAMMA ) / S END IF RETURN ELSE * * normal case * ZETA1 = ALPHA / ABSEST ZETA2 = GAMMA / ABSEST * NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ), $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 ) * * See if root is closer to zero or to ONE * TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 ) IF( TEST.GE.ZERO ) THEN * * root is close to zero, compute directly * B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF C = ZETA2*ZETA2 T = C / ( B+SQRT( ABS( B*B-C ) ) ) SINE = ZETA1 / ( ONE-T ) COSINE = -ZETA2 / T SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST ELSE * * root is closer to ONE, shift by that amount * B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF C = ZETA1*ZETA1 IF( B.GE.ZERO ) THEN T = -C / ( B+SQRT( B*B+C ) ) ELSE T = B - SQRT( B*B+C ) END IF SINE = -ZETA1 / T COSINE = -ZETA2 / ( ONE+T ) SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST END IF TMP = SQRT( SINE*SINE+COSINE*COSINE ) S = SINE / TMP C = COSINE / TMP RETURN * END IF END IF RETURN * * End of DLAIC1 * END DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. CHARACTER CMACH * .. * * Purpose * ======= * * DLAMCH determines double precision machine parameters. * * Arguments * ========= * * CMACH (input) CHARACTER*1 * Specifies the value to be returned by DLAMCH: * = 'E' or 'e', DLAMCH := eps * = 'S' or 's , DLAMCH := sfmin * = 'B' or 'b', DLAMCH := base * = 'P' or 'p', DLAMCH := eps*base * = 'N' or 'n', DLAMCH := t * = 'R' or 'r', DLAMCH := rnd * = 'M' or 'm', DLAMCH := emin * = 'U' or 'u', DLAMCH := rmin * = 'L' or 'l', DLAMCH := emax * = 'O' or 'o', DLAMCH := rmax * * where * * eps = relative machine precision * sfmin = safe minimum, such that 1/sfmin does not overflow * base = base of the machine * prec = eps*base * t = number of (base) digits in the mantissa * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise * emin = minimum exponent before (gradual) underflow * rmin = underflow threshold - base**(emin-1) * emax = largest exponent before overflow * rmax = overflow threshold - (base**emax)*(1-eps) * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, $ RND, SFMIN, SMALL, T * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLAMC2 * .. * .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, $ EMAX, RMAX, PREC * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) BASE = BETA T = IT IF( LRND ) THEN RND = ONE EPS = ( BASE**( 1-IT ) ) / 2 ELSE RND = ZERO EPS = BASE**( 1-IT ) END IF PREC = EPS*BASE EMIN = IMIN EMAX = IMAX SFMIN = RMIN SMALL = ONE / RMAX IF( SMALL.GE.SFMIN ) THEN * * Use SMALL plus a bit, to avoid the possibility of rounding * causing overflow when computing 1/sfmin. * SFMIN = SMALL*( ONE+EPS ) END IF END IF * IF( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF * DLAMCH = RMACH RETURN * * End of DLAMCH * END * ************************************************************************ * SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * DLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * The routine is based on the routine ENVRON by Malcolm and * incorporates suggestions by Gentleman and Marovich. See * * Malcolm M. A. (1972) Algorithms to reveal properties of * floating-point arithmetic. Comms. of the ACM, 15, 949-951. * * Gentleman W. M. and Marovich S. B. (1974) More on algorithms * that reveal properties of floating point arithmetic units. * Comms. of the ACM, 17, 276-277. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * Compute a = 2.0**m with the smallest positive integer m such * that * * fl( a + 1.0 ) = a. * A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 10 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 10 END IF *+ END WHILE * * Now compute b = 2.0**m with the smallest positive integer m * such that * * fl( a + b ) .gt. a. * B = 1 C = DLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = DLAMC3( A, B ) GO TO 20 END IF *+ END WHILE * * Now compute the base. a and c are neighbouring floating point * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so * their difference is beta. Adding 0.25 to c is to ensure that it * is truncated to beta and not ( beta - 1 ). * QTR = ONE / 4 SAVEC = C C = DLAMC3( C, -A ) LBETA = C + QTR * * Now determine whether rounding or chopping occurs, by adding a * bit less than beta/2 and a bit more than beta/2 to a. * B = LBETA F = DLAMC3( B / 2, -B / 100 ) C = DLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = DLAMC3( B / 2, B / 100 ) C = DLAMC3( F, A ) IF( ( LRND ) .AND. ( C.EQ.A ) ) $ LRND = .FALSE. * * Try and decide whether rounding is done in the IEEE 'round to * nearest' style. B/2 is half a unit in the last place of the two * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit * zero, and SAVEC is odd. Thus adding B/2 to A should not change * A, but adding B/2 to SAVEC should change SAVEC. * T1 = DLAMC3( B / 2, A ) T2 = DLAMC3( B / 2, SAVEC ) LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND * * Now find the mantissa, t. It should be the integer part of * log to the base beta of a, however it is safer to determine t * by powering. So we find t as the smallest positive integer for * which * * fl( beta**t + 1.0 ) = 1.0. * LT = 0 A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 30 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 RETURN * * End of DLAMC1 * END * ************************************************************************ * SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T DOUBLE PRECISION EPS, RMAX, RMIN * .. * * Purpose * ======= * * DLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * EPS (output) DOUBLE PRECISION * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) DOUBLE PRECISION * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) DOUBLE PRECISION * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. External Subroutines .. EXTERNAL DLAMC1, DLAMC4, DLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ZERO = 0 ONE = 1 TWO = 2 * * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of * BETA, T, RND, EPS, EMIN and RMIN. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) * * Start to find EPS. * B = LBETA A = B**( -LT ) LEPS = A * * Try some tricks to see whether or not this is the correct EPS. * B = TWO / 3 HALF = ONE / 2 SIXTH = DLAMC3( B, -HALF ) THIRD = DLAMC3( SIXTH, SIXTH ) B = DLAMC3( THIRD, -HALF ) B = DLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = DLAMC3( HALF, -C ) B = DLAMC3( HALF, C ) C = DLAMC3( HALF, -B ) B = DLAMC3( HALF, C ) GO TO 10 END IF *+ END WHILE * IF( A.LT.LEPS ) $ LEPS = A * * Computation of EPS complete. * * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). * Keep dividing A by BETA until (gradual) underflow occurs. This * is detected when we cannot recover the previous A. * RBASE = ONE / LBETA SMALL = ONE DO 20 I = 1, 3 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = DLAMC3( ONE, SMALL ) CALL DLAMC4( NGPMIN, ONE, LBETA ) CALL DLAMC4( NGNMIN, -ONE, LBETA ) CALL DLAMC4( GPMIN, A, LBETA ) CALL DLAMC4( GNMIN, -A, LBETA ) IEEE = .FALSE. * IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN IF( NGPMIN.EQ.GPMIN ) THEN LEMIN = NGPMIN * ( Non twos-complement machines, no gradual underflow; * e.g., VAX ) ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN LEMIN = NGPMIN - 1 + LT IEEE = .TRUE. * ( Non twos-complement machines, with gradual underflow; * e.g., IEEE standard followers ) ELSE LEMIN = MIN( NGPMIN, GPMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) * ( Twos-complement machines, no gradual underflow; * e.g., CYBER 205 ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. $ ( GPMIN.EQ.GNMIN ) ) THEN IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT * ( Twos-complement machines with gradual underflow; * no known machine ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )LEMIN END IF *** * * Assume IEEE arithmetic if we found denormalised numbers above, * or if arithmetic seems to round in the IEEE style, determined * in routine DLAMC1. A true IEEE machine should have both things * true; however, faulty machines may have one or the other. * IEEE = IEEE .OR. LIEEE1 * * Compute RMIN by successive division by BETA. We could compute * RMIN as BASE**( EMIN - 1 ), but some machines underflow during * this computation. * LRMIN = 1 DO 30 I = 1, 1 - LEMIN LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 30 CONTINUE * * Finally, call DLAMC5 to compute EMAX and RMAX. * CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) END IF * BETA = LBETA T = LT RND = LRND EPS = LEPS EMIN = LEMIN RMIN = LRMIN EMAX = LEMAX RMAX = LRMAX * RETURN * 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', $ ' EMIN = ', I8, / $ ' If, after inspection, the value EMIN looks', $ ' acceptable please comment out ', $ / ' the IF block as marked within the code of routine', $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) * * End of DLAMC2 * END * ************************************************************************ * DOUBLE PRECISION FUNCTION DLAMC3( A, B ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. DOUBLE PRECISION A, B * .. * * Purpose * ======= * * DLAMC3 is intended to force A and B to be stored prior to doing * the addition of A and B , for use in situations where optimizers * might hold one of these in a register. * * Arguments * ========= * * A, B (input) DOUBLE PRECISION * The values A and B. * * ===================================================================== * * .. Executable Statements .. * DLAMC3 = A + B * RETURN * * End of DLAMC3 * END * ************************************************************************ * SUBROUTINE DLAMC4( EMIN, START, BASE ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER BASE, EMIN DOUBLE PRECISION START * .. * * Purpose * ======= * * DLAMC4 is a service routine for DLAMC2. * * Arguments * ========= * * EMIN (output) EMIN * The minimum exponent before (gradual) underflow, computed by * setting A = START and dividing by BASE until the previous A * can not be recovered. * * START (input) DOUBLE PRECISION * The starting point for determining EMIN. * * BASE (input) INTEGER * The base of the machine. * * ===================================================================== * * .. Local Scalars .. INTEGER I DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Executable Statements .. * A = START ONE = 1 RBASE = ONE / BASE ZERO = 0 EMIN = 1 B1 = DLAMC3( A*RBASE, ZERO ) C1 = A C2 = A D1 = A D2 = A *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 10 CONTINUE IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. $ ( D2.EQ.A ) ) THEN EMIN = EMIN - 1 A = B1 B1 = DLAMC3( A / BASE, ZERO ) C1 = DLAMC3( B1*BASE, ZERO ) D1 = ZERO DO 20 I = 1, BASE D1 = D1 + B1 20 CONTINUE B2 = DLAMC3( A*RBASE, ZERO ) C2 = DLAMC3( B2 / RBASE, ZERO ) D2 = ZERO DO 30 I = 1, BASE D2 = D2 + B2 30 CONTINUE GO TO 10 END IF *+ END WHILE * RETURN * * End of DLAMC4 * END * ************************************************************************ * SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL IEEE INTEGER BETA, EMAX, EMIN, P DOUBLE PRECISION RMAX * .. * * Purpose * ======= * * DLAMC5 attempts to compute RMAX, the largest machine floating-point * number, without overflow. It assumes that EMAX + abs(EMIN) sum * approximately to a power of 2. It will fail on machines where this * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, * EMAX = 28718). It will also fail if the value supplied for EMIN is * too large (i.e. too close to zero), probably with overflow. * * Arguments * ========= * * BETA (input) INTEGER * The base of floating-point arithmetic. * * P (input) INTEGER * The number of base BETA digits in the mantissa of a * floating-point value. * * EMIN (input) INTEGER * The minimum exponent before (gradual) underflow. * * IEEE (input) LOGICAL * A logical flag specifying whether or not the arithmetic * system is thought to comply with the IEEE standard. * * EMAX (output) INTEGER * The largest exponent before overflow * * RMAX (output) DOUBLE PRECISION * The largest machine floating-point number. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP DOUBLE PRECISION OLDY, RECBAS, Y, Z * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * * First compute LEXP and UEXP, two powers of 2 that bound * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum * approximately to the bound that is closest to abs(EMIN). * (EMAX is the exponent of the required number RMAX). * LEXP = 1 EXBITS = 1 10 CONTINUE TRY = LEXP*2 IF( TRY.LE.( -EMIN ) ) THEN LEXP = TRY EXBITS = EXBITS + 1 GO TO 10 END IF IF( LEXP.EQ.-EMIN ) THEN UEXP = LEXP ELSE UEXP = TRY EXBITS = EXBITS + 1 END IF * * Now -LEXP is less than or equal to EMIN, and -UEXP is greater * than or equal to EMIN. EXBITS is the number of bits needed to * store the exponent. * IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN EXPSUM = 2*LEXP ELSE EXPSUM = 2*UEXP END IF * * EXPSUM is the exponent range, approximately equal to * EMAX - EMIN + 1 . * EMAX = EXPSUM + EMIN - 1 NBITS = 1 + EXBITS + P * * NBITS is the total number of bits needed to store a * floating-point number. * IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN * * Either there are an odd number of bits used to store a * floating-point number, which is unlikely, or some bits are * not used in the representation of numbers, which is possible, * (e.g. Cray machines) or the mantissa has an implicit bit, * (e.g. IEEE machines, Dec Vax machines), which is perhaps the * most likely. We have to assume the last alternative. * If this is true, then we need to reduce EMAX by one because * there must be some way of representing zero in an implicit-bit * system. On machines like Cray, we are reducing EMAX by one * unnecessarily. * EMAX = EMAX - 1 END IF * IF( IEEE ) THEN * * Assume we are on an IEEE machine which reserves one exponent * for infinity and NaN. * EMAX = EMAX - 1 END IF * * Now create RMAX, the largest machine number, which should * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . * * First compute 1.0 - BETA**(-P), being careful that the * result is less than 1.0 . * RECBAS = ONE / BETA Z = BETA - ONE Y = ZERO DO 20 I = 1, P Z = Z*RECBAS IF( Y.LT.ONE ) $ OLDY = Y Y = DLAMC3( Y, Z ) 20 CONTINUE IF( Y.GE.ONE ) $ Y = OLDY * * Now multiply by BETA**EMAX to get RMAX. * DO 30 I = 1, EMAX Y = DLAMC3( Y*BETA, ZERO ) 30 CONTINUE * RMAX = Y RETURN * * End of DLAMC5 * END DOUBLE PRECISION FUNCTION DLAPY2( X, Y ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. DOUBLE PRECISION X, Y * .. * * Purpose * ======= * * DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary * overflow. * * Arguments * ========= * * X (input) DOUBLE PRECISION * Y (input) DOUBLE PRECISION * X and Y specify the values x and y. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION W, XABS, YABS, Z * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * XABS = ABS( X ) YABS = ABS( Y ) W = MAX( XABS, YABS ) Z = MIN( XABS, YABS ) IF( Z.EQ.ZERO ) THEN DLAPY2 = W ELSE DLAPY2 = W*SQRT( ONE+( Z / W )**2 ) END IF RETURN * * End of DLAPY2 * END SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER SIDE INTEGER INCV, LDC, M, N DOUBLE PRECISION TAU * .. * .. Array Arguments .. DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) * .. * * Purpose * ======= * * DLARF applies a real elementary reflector H to a real m by n matrix * C, from either the left or the right. H is represented in the form * * H = I - tau * v * v' * * where tau is a real scalar and v is a real vector. * * If tau = 0, then H is taken to be the unit matrix. * * Arguments * ========= * * SIDE (input) CHARACTER*1 * = 'L': form H * C * = 'R': form C * H * * M (input) INTEGER * The number of rows of the matrix C. * * N (input) INTEGER * The number of columns of the matrix C. * * V (input) DOUBLE PRECISION array, dimension * (1 + (M-1)*abs(INCV)) if SIDE = 'L' * or (1 + (N-1)*abs(INCV)) if SIDE = 'R' * The vector v in the representation of H. V is not used if * TAU = 0. * * INCV (input) INTEGER * The increment between elements of v. INCV <> 0. * * TAU (input) DOUBLE PRECISION * The value tau in the representation of H. * * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) * On entry, the m by n matrix C. * On exit, C is overwritten by the matrix H * C if SIDE = 'L', * or C * H if SIDE = 'R'. * * LDC (input) INTEGER * The leading dimension of the array C. LDC >= max(1,M). * * WORK (workspace) DOUBLE PRECISION array, dimension * (N) if SIDE = 'L' * or (M) if SIDE = 'R' * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. External Subroutines .. EXTERNAL DGEMV, DGER * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Executable Statements .. * IF( LSAME( SIDE, 'L' ) ) THEN * * Form H * C * IF( TAU.NE.ZERO ) THEN * * w := C' * v * CALL DGEMV( 'Transpose', M, N, ONE, C, LDC, V, INCV, ZERO, $ WORK, 1 ) * * C := C - v * w' * CALL DGER( M, N, -TAU, V, INCV, WORK, 1, C, LDC ) END IF ELSE * * Form C * H * IF( TAU.NE.ZERO ) THEN * * w := C * v * CALL DGEMV( 'No transpose', M, N, ONE, C, LDC, V, INCV, $ ZERO, WORK, 1 ) * * C := C - w * v' * CALL DGER( M, N, -TAU, WORK, 1, V, INCV, C, LDC ) END IF END IF RETURN * * End of DLARF * END SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, $ T, LDT, C, LDC, WORK, LDWORK ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER DIRECT, SIDE, STOREV, TRANS INTEGER K, LDC, LDT, LDV, LDWORK, M, N * .. * .. Array Arguments .. DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), $ WORK( LDWORK, * ) * .. * * Purpose * ======= * * DLARFB applies a real block reflector H or its transpose H' to a * real m by n matrix C, from either the left or the right. * * Arguments * ========= * * SIDE (input) CHARACTER*1 * = 'L': apply H or H' from the Left * = 'R': apply H or H' from the Right * * TRANS (input) CHARACTER*1 * = 'N': apply H (No transpose) * = 'T': apply H' (Transpose) * * DIRECT (input) CHARACTER*1 * Indicates how H is formed from a product of elementary * reflectors * = 'F': H = H(1) H(2) . . . H(k) (Forward) * = 'B': H = H(k) . . . H(2) H(1) (Backward) * * STOREV (input) CHARACTER*1 * Indicates how the vectors which define the elementary * reflectors are stored: * = 'C': Columnwise * = 'R': Rowwise * * M (input) INTEGER * The number of rows of the matrix C. * * N (input) INTEGER * The number of columns of the matrix C. * * K (input) INTEGER * The order of the matrix T (= the number of elementary * reflectors whose product defines the block reflector). * * V (input) DOUBLE PRECISION array, dimension * (LDV,K) if STOREV = 'C' * (LDV,M) if STOREV = 'R' and SIDE = 'L' * (LDV,N) if STOREV = 'R' and SIDE = 'R' * The matrix V. See further details. * * LDV (input) INTEGER * The leading dimension of the array V. * If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); * if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); * if STOREV = 'R', LDV >= K. * * T (input) DOUBLE PRECISION array, dimension (LDT,K) * The triangular k by k matrix T in the representation of the * block reflector. * * LDT (input) INTEGER * The leading dimension of the array T. LDT >= K. * * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) * On entry, the m by n matrix C. * On exit, C is overwritten by H*C or H'*C or C*H or C*H'. * * LDC (input) INTEGER * The leading dimension of the array C. LDA >= max(1,M). * * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K) * * LDWORK (input) INTEGER * The leading dimension of the array WORK. * If SIDE = 'L', LDWORK >= max(1,N); * if SIDE = 'R', LDWORK >= max(1,M). * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. CHARACTER TRANST INTEGER I, J * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DCOPY, DGEMM, DTRMM * .. * .. Executable Statements .. * * Quick return if possible * IF( M.LE.0 .OR. N.LE.0 ) $ RETURN * IF( LSAME( TRANS, 'N' ) ) THEN TRANST = 'T' ELSE TRANST = 'N' END IF * IF( LSAME( STOREV, 'C' ) ) THEN * IF( LSAME( DIRECT, 'F' ) ) THEN * * Let V = ( V1 ) (first K rows) * ( V2 ) * where V1 is unit lower triangular. * IF( LSAME( SIDE, 'L' ) ) THEN * * Form H * C or H' * C where C = ( C1 ) * ( C2 ) * * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) * * W := C1' * DO 10 J = 1, K CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 10 CONTINUE * * W := W * V1 * CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N, $ K, ONE, V, LDV, WORK, LDWORK ) IF( M.GT.K ) THEN * * W := W + C2'*V2 * CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K, $ ONE, C( K+1, 1 ), LDC, V( K+1