This is the mail archive of the
gcc-help@gcc.gnu.org
mailing list for the GCC project.
I need some help with an stability analysis subroutine package (second part)
- To: help-gcc at gnu dot org
- Subject: I need some help with an stability analysis subroutine package (second part)
- From: Sonia Cortassa <intech at cpsarg dot com>
- Date: Thu, 30 Nov 2000 09:59:23 -0300
Hello!
I am writing because I am currently using a software for continuation and
bifurcation analysis in ordinary differential equations. Its name is AUTO
and was created in 1986 by E. Doedle at Caltech. I must admit I am not a
program developer, neither am I familiar with fortran 77 synthaxis. I am
just able to run this software as a tool for my studies of models of
biological systems. Some people have helped me to install this software
in a GNU/Linux operating system that I am now using with g77.
My specific question is the following: I need to use it to evaluate
stability in a large model (11 variables) and the problem is that the
output files only display the first six variables. I know it is
computing all the eleven variables since in the unit where the eigen
values are stored there are 11 of them but the units 6, 7, 8 and 9 where
the actual value of the variables is printed does not display them all.
My question is how can I modify the software, since I know it is able to
deal with up to 25 variables, for it to display all the variables I need?
Attached you'll find the source code of the preprocessor of the AUTO
package called "AUTPRPD.F" and the subroutines library
"AUTLIBD.F" (due to size reasons this file was divided in two
files AUTLIBDa .F and AUTLIBDb.F this second file is attached to a second
part of the message)
Thank you very much in advance for your help
Sonia Cortassa
############################################
Sonia Cortassa, PhD
Instituto Tecnologico de Chascomus,
Casilla Correo 164
7130 - Chascomus, Pvcia Buenos Aires
Tel: 54 2241 424049
Fax: 54 2241 424048
E-mail: maaon@criba.edu.ar
############################################
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C this document contains the second part of the library file AUTLIBD.F
C
C A U T O 8 6
C
C
C A Subroutine Package for the Bifurcation Analysis of
C Autonomous Systems of Ordinary Differential Equations.
C
C
C Author : Eusebius Doedel
C Applied Mathematics 217-50
C California Institute of Technology
C Pasadena, California 91125
C
C (Further distribution requires notification of the author.)
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Version January 1986 (Partially Vectorized).
C Updated : August 1987
C
C For Documentation see the AUTO 86 User Manual.
C
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C Detection and Location of Bifurcations in Boundary Value Problems
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C
C ---------- ------
SUBROUTINE LCSPBV(FNCS,FUNI,BCNI,ICNI,ISTOP,ITP,SP1,NITPS,IBR,
* NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,M1U,UPS,
* UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,DFDU,DFDP,
* UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD,P0,P1,POIN,
* EV,WKEV)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C This subroutine uses the Secant method to accurately locate limit
C points, bifurcation points, and zero(es) of function(s) from USZR.
C Such points are located as points on a solution branch where the
C EXTERNAL function FNCS changes sign.
C It involves calling the basic solution subroutines CONTBV and SOLVBV
C with decreasing values of RDS (stepsize along branch).
C The point is assumed to have been found with sufficient accuracy if
C the ratio between RDS and the user supplied value of DS is less than
C the user-supplied tolerance EPSS.
C This subroutine is called from CNRLB, which controls the computation
C of branches of solutions to general boundary value problems.
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLTHT/ THETAL(20),THETAU
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLEPS/ EPSL(20),EPSU,EPSS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDET/ DETGE,NINS
C
EXTERNAL FUNI,BCNI,ICNI
C
COMPLEX*16 EV(NDIM)
CSGLE COMPLEX*8 EV(NDIM)
C
LOGICAL CHNG
C
DIMENSION AA(M1AA,M2AA,1),BB(M1BB,M2BB,1),CC(M1CC,1),DD(M1DD,1)
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),UOLDPS(M1U,1)
DIMENSION RHSA(M1U,1),RHSD(1),DUPS(M1U,1),TM(1),DTM(1),U(1),F(1)
DIMENSION DFDU(M1DF,1),DFDP(M1DF,1),WBRBD(1),IWBRBD(1),IR(1),IC(1)
DIMENSION UBC0(1),UBC1(1),DBC(M1BC,1),UICD(1),FICD(1),DICD(M1IC,1)
DIMENSION P0(NDIM,1),P1(NDIM,1),POIN(NDIM,1),WKEV(NDIM)
C
SP10=SP1
C
C Check for zero.
C
SP1=FNCS(CHNG,FUNI,BCNI,ICNI,ISTOP,ITP,NITPS,P0,P1,POIN,EV,WKEV,
* IBR,NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,
* M1U,UPS,UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,
* DFDU,DFDP,UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
C
PSP1=SP10*SP1
NTOTP1=NTOT+1
IF(PSP1.GE.ZERO .OR. (.NOT. CHNG))RETURN
C
C Compute next RDS by a perturbed Secant method :
C
RDS=RDSOLD
NITSP1=0
1 DSP1=SP10-SP1
IF(DSP1.EQ.ZERO)RDS=ZERO
IF(DSP1.NE.ZERO)RDS=SP1/DSP1*RDS
RDS=(ONE+HMACH)*RDS
C
C If requested write additional output on unit 9 :
C
IF(IID.GE.2)WRITE(9,102)NITSP1,RDS
C
C Return if tolerance has been met :
C
RRDS=DABS(RDS)/(ONE+DABS(DS))
CSGLE RRDS= ABS(RDS)/(ONE+ ABS(DS))
IF(RRDS.LT.EPSS) THEN
ITP=-1
RETURN
ENDIF
C
CALL CONTBV(FUNI,RDS,M1U,UPS,UOLDPS,UDOTPS,UPOLDP,U,F,M1DF,DFDU,
* DFDP,DTM)
CALL SOLVBV(FUNI,BCNI,ICNI,ISTOP,RDS,NITPS,IBR,NTOT,M1AA,M2AA,AA,
* M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,M1U,UPS,UOLDPS,UDOTPS,UPOLDP,
* RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,DFDU,DFDP,UBC0,UBC1,M1BC,DBC,
* UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
IF(ISTOP.NE.0)THEN
SP1=ZERO
RETURN
ENDIF
C
C Check for zero.
C
SP10=SP1
SP1=FNCS(CHNG,FUNI,BCNI,ICNI,ISTOP,ITP,NITPS,P0,P1,POIN,EV,WKEV,
* IBR,NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,
* M1U,UPS,UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,
* DFDU,DFDP,UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
C
NITSP1=NITSP1+1
IF(NITSP1.LE.ITMX)GOTO 1
C
WRITE(9,101)IBR,NTOTP1
SP1=ZERO
101 FORMAT(' *** POSSIBLE SINGULAR POINT (BRANCH ',I3,' POINT ',
* I4,')')
102 FORMAT(' * DETECTION OF SINGULAR POINT : ITERATION ',I3,
* ' STEPSIZE =',E11.3)
C
RETURN
END
C
C ------ ---------
DOUBLE PRECISION
CSGLE REAL
*FUNCTION FNLPBV(CHNG,FUNI,BCNI,ICNI,ISTOP,ITP,NITPS,P0,P1,POIN,EV,
* WKEV,IBR,NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,
* M1U,UPS,UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,
* DFDU,DFDP,UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C RETURNS A QUANTITY THAT CHANGES SIGN AT A LIMIT POINT (BVP)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLEPS/ EPSL(20),EPSU,EPSS
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
COMPLEX*16 EV(NDIM)
CSGLE COMPLEX*8 EV(NDIM)
C
LOGICAL CHNG
C
EXTERNAL FUNI,BCNI,ICNI
C
DIMENSION AA(M1AA,M2AA,1),BB(M1BB,M2BB,1),CC(M1CC,1),DD(M1DD,1)
DIMENSION UPS(M1U,1),UDOTPS(M1U,NROW),UPOLDP(M1U,1),UOLDPS(M1U,1)
DIMENSION RHSA(M1U,NROW),RHSD(NDRHS),DUPS(M1U,1),TM(1),DTM(1),U(1)
DIMENSION UBC0(1),UBC1(1),DBC(M1BC,1),UICD(1),FICD(1),DICD(M1IC,1)
DIMENSION IR(1),IC(1),IWBRBD(1),F(1),DFDU(M1DF,1),DFDP(M1DF,1)
DIMENSION P0(NDIM,1),P1(NDIM,1),POIN(NDIM,1),WKEV(NDIM),WBRBD(1)
C
C Find the direction vector.
C
IFST=0
IF(IID.LT.4)THEN
IIBR=0
ELSE IF(IID.EQ.4)THEN
IIBR=1
ELSE
IIBR=2
ENDIF
CALL BRBD(NTST,NROW,NCLM,M1AA,M2AA,AA,NFPAR,M1BB,M2BB,BB,NRC,
* M1CC,CC,M1DD,DD,M1U,RHSA,RHSD,WBRBD,IIBR,IFST,IR,IC,IWBRBD,-1)
C
DO 1 I=1,NDIM
UDOTPS(NTST+1,I)=RHSD(I)
1 CONTINUE
DO 2 I=1,NFPAR
RLDOT(I)=RHSD(NDIM+I)
2 CONTINUE
C
DO 4 J=1,NTST
DO 3 I=1,NROW
UDOTPS(J,I)=RHSA(J,I)
3 CONTINUE
4 CONTINUE
C
C Scale the direction vector.
C
CALL SCALEB(M1U,UDOTPS,RLDOT,DTM)
IF(IID.GE.2)WRITE(9,101)RLDOT(1)
C
C Set the quantity to be returned.
C
FNLPBV=RLDOT(1)
IF(IABS(ISP).EQ.3)FNLPBV=RLDOT(2)
CHNG=.TRUE.
C
101 FORMAT(' LIMIT POINT FUNCTION = ',E11.3)
C
RETURN
END
C
C ------ ---------
DOUBLE PRECISION
CSGLE REAL
*FUNCTION FNBPBV(CHNG,FUNI,BCNI,ICNI,ISTOP,ITP,NITPS,P0,P1,POIN,EV,
* WKEV,IBR,NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,
* M1U,UPS,UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,
* DFDU,DFDP,UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLEPS/ EPSL(20),EPSU,EPSS
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDET/ DETGE,NINS
C
COMPLEX*16 EV(NDIM)
CSGLE COMPLEX*8 EV(NDIM)
C
LOGICAL CHNG
C
EXTERNAL FUNI,BCNI,ICNI
C
DIMENSION AA(M1AA,M2AA,1),BB(M1BB,M2BB,1),CC(M1CC,1),DD(M1DD,1)
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),UOLDPS(M1U,1)
DIMENSION RHSA(M1U,1),RHSD(1),DUPS(M1U,1),TM(1),DTM(1),U(1)
DIMENSION UBC0(1),UBC1(1),DBC(M1BC,1),UICD(1),FICD(1),DICD(M1IC,1)
DIMENSION IR(1),IC(1),IWBRBD(1),F(1),DFDU(M1DF,1),DFDP(M1DF,1)
DIMENSION P0(NDIM,1),P1(NDIM,1),POIN(NDIM,1),WKEV(NDIM),WBRBD(1)
C
C Save the determinant of the reduced system.
C
DET=DETGE
C
C Compute the determinant of P1.
CALL GE(NDIM,NDIM,P1,0,1,U,1,F,IR,IC)
C
C Set the determinant of the normalized reduced system.
C
IF(DETGE.NE.ZERO)THEN
FNBPBV=DET/DETGE
CHNG=.TRUE.
ELSE
FNBPBV=ZERO
CHNG=.FALSE.
ENDIF
C
IF(IID.GE.2)WRITE(9,101)FNBPBV
101 FORMAT(' BIFURCATION FUNCTION = ',E11.3)
C
RETURN
END
C
C ------ ---------
DOUBLE PRECISION
CSGLE REAL
*FUNCTION FNSPBV(CHNG,FUNI,BCNI,ICNI,ISTOP,ITP,NITPS,P0,P1,POIN,EV,
* WKEV,IBR,NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,
* M1U,UPS,UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,
* DFDU,DFDP,UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C This function returns a quantity that changes sign when a complex
C pair of eigenvalues of the linearized Poincare map moves in or out
C of the unit circle.
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLEPS/ EPSL(20),EPSU,EPSS
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDET/ DETGE,NINS
C
COMPLEX*16 EV(NDIM),ZTMP
CSGLE COMPLEX*8 EV(NDIM),ZTMP
C
LOGICAL CHNG
C
EXTERNAL FUNI,BCNI,ICNI
C
DIMENSION AA(M1AA,M2AA,1),BB(M1BB,M2BB,1),CC(M1CC,1),DD(M1DD,1)
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),UOLDPS(M1U,1)
DIMENSION RHSA(M1U,1),RHSD(1),DUPS(M1U,1),TM(1),DTM(1),U(1)
DIMENSION UBC0(1),UBC1(1),DBC(M1BC,1),UICD(1),FICD(1),DICD(M1IC,1)
DIMENSION IR(1),IC(1),IWBRBD(1),F(1),DFDU(M1DF,1),DFDP(M1DF,1)
DIMENSION P0(NDIM,1),P1(NDIM,1),POIN(NDIM,1),WKEV(NDIM),WBRBD(1)
C
C Initialize.
C
FNSPBV=ZERO
D=ZERO
CHNG=.FALSE.
C
C Compute the linearization of the Poincare map.
C
CALL POINC(NDIM,P0,P1,POIN,IID,IR,IC)
C
C Compute the Floquet multipliers.
C
CALL EIG(NDIM,NDIM,POIN,EV,WKEV,IER)
C
C Order the Floquet multipliers by distance from z=1.
C
DO 2 I=1,NDIM-1
AMIN=RLARGE
DO 1 J=I,NDIM
AZM1=CDABS( EV(J) - ONE )
CSGLE AZM1= CABS( EV(J) - ONE )
IF(AZM1.GT.AMIN)GOTO 1
AMIN=AZM1
LOC=J
1 CONTINUE
IF(LOC.NE.I) THEN
ZTMP=EV(LOC)
EV(LOC)=EV(I)
EV(I)=ZTMP
ENDIF
2 CONTINUE
C
C Find eigenvalue closest to the unit circle
C (excluding the eigenvalue at z=1).
C
NINS1=IABS(ISP)-1
IF(NINS1.LT.1)NINS1=1
IF(NINS1.GT.NDIM)NINS1=NDIM
C
IF(NINS1.LT.NDIM)THEN
AMIN=RLARGE
DO 3 I=NINS1+1,NDIM
D=CDABS(EV(I)) - ONE
CSGLE D= CABS(EV(I)) - ONE
AD=DABS(D)
CSGLE AD= ABS(D)
IF(AD.GT.AMIN)GOTO 3
AMIN=AD
LOC=I
3 CONTINUE
C Interchange, to put eigenvalue in ISP'th position.
ISP1=ISP
IF(ISP.EQ.1)ISP1=2
IF(LOC.NE.ISP1) THEN
ZTMP=EV(LOC)
EV(LOC)=EV(ISP1)
EV(ISP1)=ZTMP
ENDIF
ENDIF
C
C Print error message if the Floquet multiplier at z=1 is inaccurate.
C (ISP is set to negative and detection of bifurations is discontinued)
C
AMIN=CDABS( EV(1) - ONE )
CSGLE AMIN= CABS( EV(1) - ONE )
C
IF(AMIN.GT.5.0E-2 .AND. ISP.GT.1) THEN
NTOT1=NTOT+1
IF(IID.GE.2)WRITE(9,104)
WRITE(9,105)IBR,NTOT1,(EV(I),I=1,NDIM)
NINS=0
WRITE(9,101)
ISP=-ISP
RETURN
ENDIF
C
C Restart automatic detection if the Floquet multiplier at z=1 is
C sufficiently accurate again.
C
IF(ISP.LT.0)THEN
IF(AMIN.LT.2.0E-2)THEN
WRITE(9,102)
ISP=-ISP
ELSE
NTOT1=NTOT+1
IF(IID.GE.2)WRITE(9,104)
WRITE(9,105)IBR,NTOT1,(EV(I),I=1,NDIM)
RETURN
ENDIF
ENDIF
C
C Count the number of Floquet multipliers inside the unit circle.
C
IF(NINS1.EQ.NDIM) THEN
D=ZERO
FNSPBV=D
GOTO 5
ENDIF
C
NINS2=NINS1
DO 4 I=NINS2+1,NDIM
IF(CDABS(EV(I)).LE.ONE)NINS1=NINS1+1
CSGLE IF( CABS(EV(I)).LE.ONE)NINS1=NINS1+1
4 CONTINUE
C
IF(ISP.GE.2) THEN
D=CDABS(EV(ISP)) - ONE
CSGLE D= CABS(EV(ISP)) - ONE
FNSPBV=D
IF(NINS1.NE.NINS)CHNG=.TRUE.
ENDIF
5 NINS=NINS1
IF(IID.GE.2.AND.ISP.GE.1)WRITE(9,103)D,NINS
C
C Print the Floquet multipliers.
C
NTOT1=NTOT+1
IF(NINS.EQ.NDIM)NTOT1=-NTOT1
IF(IID.GE.2)WRITE(9,104)
WRITE(9,105)IBR,NTOT1,(EV(I),I=1,NDIM)
C
101 FORMAT(' *** FLOQUET MULTIPLIER AT 1 INACCURATE')
102 FORMAT(' *** FLOQUET MULTIPLIER AT 1 ACCURATE AGAIN')
103 FORMAT(' BIFURCATION FUNCTION = ',E11.3,
* ' # OF MULTIPLIERS IN UNIT CIRCLE =',I3)
104 FORMAT(' FLOQUET MULTIPLIERS :')
105 FORMAT(' BRANCH ',I3,' POINT ',I4,2(2X,2E12.5),
* 50(/,23X,2(2X,2E12.5)))
C
RETURN
END
C
C ------ ---------
DOUBLE PRECISION
CSGLE REAL
*FUNCTION FNUZBV(CHNG,FUNI,BCNI,ICNI,ISTOP,ITP,NITPS,P0,P1,POIN,EV,
* WKEV,IBR,NTOT,M1AA,M2AA,AA,M1BB,M2BB,BB,M1CC,CC,M1DD,DD,WBRBD,
* M1U,UPS,UOLDPS,UDOTPS,UPOLDP,RHSA,RHSD,DUPS,TM,DTM,U,F,M1DF,
* DFDU,DFDP,UBC0,UBC1,M1BC,DBC,UICD,FICD,M1IC,DICD,IR,IC,IWBRBD)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLEPS/ EPSL(20),EPSU,EPSS
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDET/ DETGE,NINS
COMMON /BLUSZ/ IUZR
C
COMPLEX*16 EV(NDIM)
CSGLE COMPLEX*8 EV(NDIM)
C
LOGICAL CHNG
C
EXTERNAL FUNI,BCNI,ICNI
C
DIMENSION AA(M1AA,M2AA,1),BB(M1BB,M2BB,1),CC(M1CC,1),DD(M1DD,1)
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),UOLDPS(M1U,1)
DIMENSION RHSA(M1U,1),RHSD(1),DUPS(M1U,1),TM(1),DTM(1),U(1)
DIMENSION UBC0(1),UBC1(1),DBC(M1BC,1),UICD(1),FICD(1),DICD(M1IC,1)
DIMENSION IR(1),IC(1),IWBRBD(1),F(1),DFDU(M1DF,1),DFDP(M1DF,1)
DIMENSION P0(NDIM,1),P1(NDIM,1),POIN(NDIM,1),WKEV(NDIM),WBRBD(1)
C
FNUZBV=USZR(IUZR,NUZR,PAR)
CHNG=.TRUE.
C
IF(IID.GE.2)WRITE(9,101)FNUZBV
101 FORMAT(' USZR FUNCTION = ',E11.3)
C
RETURN
END
C
C ---------- -----
SUBROUTINE POINC(NDIM,P0,P1,POIN,IID,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Computes the linearized Poincare map. This map is extracted from the
C decomposition of the Jacobian matrix as generated by BRBD.
C
DIMENSION P0(NDIM,1),P1(NDIM,1),POIN(NDIM,1)
DIMENSION IR(1),IC(1)
C
DO 2 I=1,NDIM
DO 1 J=1,NDIM
P0(I,J)=-P0(I,J)
1 CONTINUE
2 CONTINUE
C
CALL GE(NDIM,NDIM,P1,NDIM,NDIM,POIN,NDIM,P0,IR,IC)
C
IF(IID.GT.2) THEN
WRITE(9,101)
DO 3 I=1,NDIM
WRITE(9,102)(POIN(I,J),J=1,NDIM)
3 CONTINUE
ENDIF
C
101 FORMAT(' LINEARIZED POINCARE MAP')
102 FORMAT(1X,6E21.14)
C
RETURN
END
C
C ---------- ------
SUBROUTINE TPSPBV(EV,ITP)
C
C Determines type of secondary periodic bifurcation.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLEPS/ EPSL(20),EPSU,EPSS
COMMON /BLITP/ ITPST,ITPSP,IBRSP
C
COMPLEX*16 EV(NDIM)
CSGLE COMPLEX*8 EV(NDIM)
C
C Find the eigenvalue closest to z=1.
C
LOC=1
AMIN=RLARGE
DO 1 I=1,NDIM
AZM1=CDABS( EV(I) - ONE )
CSGLE AZM1= CABS( EV(I) - ONE )
IF(AZM1.GT.AMIN)GOTO 1
AMIN=AZM1
LOC=I
1 CONTINUE
C
C Find the eigenvalue closest to the unit circle
C (excluding the eigenvalue at z=1).
C
LOC1=1
AMIN=RLARGE
DO 2 I=1,NDIM
IF(I.EQ.LOC)GOTO 2
D=CDABS(EV(I)) - ONE
CSGLE D= CABS(EV(I)) - ONE
AD=DABS(D)
CSGLE AD= ABS(D)
IF(AD.GT.AMIN)GOTO 2
AMIN=AD
LOC1=I
2 CONTINUE
C
IF(DABS(DIMAG(EV(LOC1))).GT.DSQRT(EPSS))THEN
CSGLE IF( ABS(AIMAG(EV(LOC1))).GT. SQRT(EPSS))THEN
C ** torus bifurcation
ITP=8+10*ITPST
PAR(12)=DASIN(DIMAG(EV(LOC1)))
CSGLE PAR(12)= ASIN(AIMAG(EV(LOC1)))
ELSE IF(DREAL(EV(LOC1)).LT.-HALF)THEN
CSGLE ELSE IF( REAL(EV(LOC1)).LT.-HALF)THEN
C ** period doubling
ITP=7+10*ITPST
ELSE
C ** ordinary bifurcation or something else...
ITP=6+10*ITPST
ENDIF
C
RETURN
END
C
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C Output (Boundary Value Problems)
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C
C ---------- ------
SUBROUTINE STPLBV(ISTOP,ITP,NITPS,NTOT,LAB,IBR,M1U,UPS,UDOTPS,TM,
* DTM,M1DF)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Writes the bifurcation diagram on unit 7 (Differential Equations)
C (Also controls the writing of complete solutions on unit 8).
C Every line written contains, in order, the following:
C
C IBR : The label of the branch.
C NTOT : The index of the point on the branch.
C (Points are numbered consecutively along a branch).
C If IPS=2 or 3, then the sign of NTOT indicates stability :
C - = stable , + = unstable, or unknown.
C ITP : An integer indicating the type of point :
C
C 4 ( ) : Output point (Every NPR steps along branch).
C -4 (UZ) : Output point (Zero of user function USZR).
C 5 (LP) : Limit point (fold).
C 6 (BP) : Bifurcation point.
C 7 (PD) : Period doubling bifurcation.
C 8 (TR) : Bifurcation to an invariant torus.
C 9 (EP) : End point of branch, normal termination.
C -9 (MX) : End point of branch, abnormal termination.
C
C LAB : The label of a special point.
C PAR(ICP(1)): The principal parameter.
C A : The L2-norm of the solution vector, or other measure of
C the solution (see the user-supplied parameter IPLT).
C MAX U(*) : The maxima of the first few solution components.
C PAR(ICP(*)): Further free parameters (if any).
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLITP/ ITPST,ITPSP,IBRSP
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLDET/ DETGE,NINS
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
C
DIMENSION UPS(M1U,NROW),UDOTPS(M1U,NROW)
DIMENSION TM(NTSTP1),DTM(NTSTP1),UMX(7)
C
NTOT=NTOT+1
C
C ITP is set to 4 every NPR steps along a branch of solns and the entire
C solution is written on unit 8.
C
IF(MOD(NTOT,NPR).EQ.0 .AND. MOD(ITP,10).EQ.0)ITP=4+10*ITPST
C
C Check whether limits of the bifurcation diagram have been reached :
C
IAB=IABS(IPLT)
IF(IAB.EQ.0.OR.IAB.GT.3*NDM)A=DSQRT(RNRMSQ(NDM,M1U,UPS,DTM))
CSGLE IF(IAB.EQ.0.OR.IAB.GT.3*NDM)A= SQRT(RNRMSQ(NDM,M1U,UPS,DTM))
IF(IPLT.GT.0.AND.IAB.LE.NDM)A=RMXUPS(M1U,IAB,UPS)
IF(IPLT.GT.NDM.AND.IAB.LE.2*NDM)A=RINTG(M1U,IAB-NDM,UPS,DTM)
IF(IPLT.GT.2*NDM.AND.IAB.LE.3*NDM)A=RNRM2(M1U,IAB-2*NDM,UPS,DTM)
IF(IPLT.LT.0.AND.IAB.LE.NDM)A=RMNUPS(M1U,IAB,UPS)
C
IF(ISTOP.EQ.1)THEN
C ** Maximum number of iterations reached somewhere.
ITP=-9-10*ITPST
ELSE
IF(PAR(ICP(1)).LT.RL0.OR.PAR(ICP(1)).GT.RL1
* .OR. A.LT.A0.OR.A.GT.A1 .OR. NTOT.GE.NMX)THEN
ISTOP=1
ITP=9+10*ITPST
ENDIF
ENDIF
C
C All special points receive label:
C
LAB1=0
IF(MOD(ITP,10).NE.0) THEN
LAB=LAB+1
LAB1=LAB
ENDIF
C
C Compute maxima of solution components.
C
N2=NDM
IF(N2.GT.7)N2=7
DO 1 I=1,N2
ITMP=I
UMX(I)=RMXUPS(M1U,ITMP,UPS)
1 CONTINUE
C
C Determine stability, and write output on units 7 and 8.
C
IBR1=IBR
NTOT1=NTOT
IF( (IPS.EQ.2.OR.IPS.EQ.3.OR.IPS.EQ.6.OR.IPS.EQ.12.OR.IPS.EQ.13)
* .AND.IABS(ISW).NE.2)THEN
IBR1=-IBR
IF(NINS.EQ.NDIM)NTOT1=-NTOT
ENDIF
CALL WRLINE(IBR1,NTOT1,ITP,LAB1,A,UMX)
C
C Write plotting and restart data on unit 8.
C
IF(IRS.NE.0.AND.NTOT.EQ.1)RETURN
IF(MOD(ITP,10).NE.0)
* CALL WRTBV8(ITP,NTOT,LAB,IBR,M1U,UPS,UDOTPS,TM,DTM)
C
RETURN
END
C
C ---------- ------
SUBROUTINE WRTBV8(ITP,NTOT,LAB,IBR,M1U,UPS,UDOTPS,TM,DTM)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Writes plotting and restart data on unit 8, viz.:
C (1) data identifying the corresponding point on unit 7,
C (2) the complete solution,
C (3) the direction of the branch.
C
C Specifically the following is written:
C
C IBR : The index of the branch.
C NTOT : The index of the point.
C ITP : The type of point (see STPLBV above).
C LAB : The label of the point.
C NFPAR : The number of free parameters used in the computation.
C ISW : The value of ISW used in the computation.
C NTPL : The number of points in the time interval [0,1] for which
C solution values are wriiten.
C NAR : The number of values written per point.
C (NAR=NDIM+1, since T and U(i), i=1,..,NDIM are written).
C NROWPR: The number of lines printed following the identifying line
C and before the next data set or the end of the file.
C (Used for quickly skipping a data set when searching).
C NTST : The number of time intervals used in the discretization.
C NCOL : The number of collocation points used.
C ICP : The indices of the free parameters in PAR(.).
C
C Following the above described identifying line there are NTPL lines
C containing :
C T , U-1(T) , U-2(T) , ... , U-NDIM(T),
C where NDIM is the dimension of the system of differential equations.
C
C Following this is a line containing
C RL-dot(i) , i=1,NFPAR,
C
C and following this are NTPL lines each containing
C U-dot-1(T), U-dot-2(T), ... , U-dot-NDIM(T).
C
C Finally the parameter values PAR(i) , i=1,NPAR, are written.
C
C Above, RL-dot(.) and U-dot(.) specify the direction of the branch.
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLDET/ DETGE,NINS
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLTIM/ TSETUB,TCONPA,TCONRH,TINFPA,TREDUC,TWR8
C
DIMENSION UPS(M1U,NROW),UDOTPS(M1U,NROW),TM(NTSTP1),DTM(NTSTP1)
C
C Set initial time (for timing of this subroutine).
C
CALL AUTIM0(TIME0)
C
C Write information identifying the solution :
C
NTSTP1=NTST+1
NTPL=NCOL*NTST+1
NAR=NDIM+1
NRD=2+NDIM/7+NDIM/8
NROWPR=NRD*(NCOL*NTST+1) + NFPAR/8+1 + NPAR/8+1
WRITE(8,101)IBR,NTOT,ITP,LAB,NFPAR,ISW,NTPL,NAR,NROWPR,NTST,NCOL
* ,(ICP(I),I=1,NFPAR)
C
C Write the entire solution on unit 8 :
C
DO 2 J=1,NTST
RN=ONE/NCOL
DO 1 I=1,NCOL
K1=(I-1)*NDIM+1
K2=I*NDIM
T=TM(J)+(I-1)*RN*DTM(J)
WRITE(8,102)T,(UPS(J,K),K=K1,K2)
1 CONTINUE
2 CONTINUE
WRITE(8,102)TM(NTSTP1),(UPS(NTSTP1,I),I=1,NDIM)
C
C Store the direction of the branch:
C
WRITE(8,102)(RLDOT(I),I=1,NFPAR)
DO 4 J=1,NTST
DO 3 I=1,NCOL
K1=(I-1)*NDIM+1
K2=I*NDIM
WRITE(8,102)(UDOTPS(J,K),K=K1,K2)
3 CONTINUE
4 CONTINUE
WRITE(8,102)(UDOTPS(NTSTP1,K),K=1,NDIM)
C
C Write the parameter values.
C
WRITE(8,102)(PAR(I),I=1,NPAR)
101 FORMAT(11I5,20I3)
102 FORMAT(4X,1P7E18.10)
C
C Determine the time spent in this subroutine.
C
CALL AUTIM1(TIME1)
TWR8=TWR8 + TIME1-TIME0
C
RETURN
END
C
C ---------- ------
SUBROUTINE WRTBV9(NITPS,IBR,NTOT,M1U,UPS,TM,DTM)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Writes additional output on unit 9.
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION UPS(M1U,NROW),TM(NTSTP1),DTM(NTSTP1)
C
IAB=IABS(IPLT)
IF(IAB.EQ.0.OR.IAB.GT.NDIM)A=DSQRT(RNRMSQ(NDM,M1U,UPS,DTM))
CSGLE IF(IAB.EQ.0.OR.IAB.GT.NDIM)A= SQRT(RNRMSQ(NDM,M1U,UPS,DTM))
IF(IPLT.GT.0.AND.IAB.LE.NDIM)A=RMXUPS(M1U,IAB,UPS)
IF(IPLT.LT.0.AND.IAB.LE.NDIM)A=RMNUPS(M1U,IAB,UPS)
IF(IID.GE.2)WRITE(9,101)IBR,NTOT+1,NITPS,A,(RL(I),I=1,NFPAR)
C
IF(.NOT.(IID.GE.5))RETURN
C
WRITE(9,103)
DO 2 J=1,NTST
RN=ONE/NCOL
DO 1 I=1,NCOL
T=TM(J)+(I-1)*RN*DTM(J)
K1=(I-1)*NDIM+1
K2=I*NDIM
WRITE(9,102)T,(UPS(J,K),K=K1,K2)
1 CONTINUE
2 CONTINUE
C
NTSTP1=NTST+1
WRITE(9,102)TM(NTSTP1),(UPS(NTSTP1,I),I=1,NDIM)
C
101 FORMAT(' BRANCH ',I2,' N=',I4,1X,'IT=',I2,1X,'A=',1PE15.8,
* 1X,' PAR= ',1P5E15.8,3(/,47X,1P5E15.8))
102 FORMAT(1X,1P7E18.10)
103 FORMAT(' UPS :')
C
RETURN
END
C
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C Linear Equation Solver
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE BRBD(NA,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* NRC,MC1,C,MD1,D,MFA1,FA,FC,WKDR,IDB,IFST,IR,IC,IWKDR,NLLV)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Solves linear systems with matrix profile:
C
C -----------------------------------------------
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C !XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX !XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C ! XXXXXXXXXX!XX!
C -----------------------------------------------
C !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX!
C !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX!
C !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX!
C !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX!
C -----------------------------------------------
C
C
C partioned as
C
C
C ---------
C ! ! ! ! ! ! !
C ! A !B! ! XA ! ! FA !
C ! ! ! . ! ! = ! ! .
C !-----!-! !----! !----!
C ! C !D! ! XC ! ! FC !
C ! ! ! ! XC ! ! FC !
C ---------
C
C
C Input parameters :
C
C NA number of blocks in A,
C NRA number of rows in each block of A,
C NCA number of columns in each block of A,
C MA1 first dimension of A from DIMENSION statement,
C MA2 second dimension of A from DIMENSION statement,
C A the matrix in the schematic representation above,
C
C NCB number of columns in each block of B,
C (note that B is also three dimensional),
C MB1 first dimension of B from DIMENSION statement,
C MB2 second dimension of B from DIMENSION statement,
C B the matrix in the schema above,
C
C NRC the number of rows of the two dimensional matrix C,
C MC1 the first dimension of C from DIMENSION statement,
C C the matrix C in the schema above,
C
C MD1 the first dimension of D from DIMENSION statement,
C D the matrix D above,
C
C MFA1 the first dimension of FA from DIMENSION statement,
C FA part of the right hand side vector,
C (note that FA is also two dimensional),
C FC part of the right hand side vector.
C
C WKDR: A one dimensional array used as workspace.
C This array should be dimensioned at least
C
C (4*NOV+NCB+NRC+1)*NOV*NA+(NRC+NOV)**2+2*(NRC+2*NOV)+NRC*NOV
C
C with NOV defined by
C
C NA*(NCA-NRA)+NCB-NRC
C NOV = -------------------- .
C NA-1
C
C IWKDR: Integer workspace array of dimension at least 3*NOV*(NA-1)+NA.
C
C IFST = 1 on first call,
C = 0 on subsequent calls with the same right hand side.
C (WKDR,IWKDR should not be modified between such calls).
C
C IDB = 0 no debug output,
C = 1,2 debug output on unit 9,
C = 3 print residuals for test problem (see PRINT2),
C
C IR, IC: Two integer arrays of dimension at least NRC+NOV.
C
C NLLV : If NLLV>0 then the system is assumed to have a NLLV-
C dimensional nullspace.
C In this case a null vector will be returned.
C If NLLVC = -1 then the system will be solved with zero right
C hand side, except for the last equation, for which the right
C hand side entry will be set to 1 (i.e., the last entry of FC
C will be set to 1, otherwise FA and FC are zero).
C If the linear system is the same as in the preceding call
C then IFST=0 may be used even if NLLV is nonzero.
C
C Returned values :
C
C FA Part of solution vector corresponding to XA in the diagram.
C FC Part of solution vector corresponding to XC in the diagram.
C
C Note : The number of columns of overlap for every two consecutive
C blocks should be equal to the number NOV defined above.
C
LOGICAL ERBRBD
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLBND/ NAM1,NAP1,NXE
COMMON /BLCNT/ NDECOM,NBCKSB
C
DIMENSION A(MA1,MA2,NCA),B(MB1,MB2,NCB),C(MC1,NRC)
DIMENSION D(MD1,NCB),WKDR(NWBR)
DIMENSION FA(MFA1,NRA),FC(NDRHS)
DIMENSION IR(NDIRC),IC(NDIRC),IWKDR(NIWBR)
C
C Check for consistency of data.
C
IF(ERBRBD(NA,NRA,NCA,NCB,NRC))STOP
NOV=(NA*(NCA-NRA)+NCB-NRC)/(NA-1)
NOV2=2*NOV
C
C Print debug output.
C
IF(IDB.GE.2)THEN
WRITE(9,101)IFST,NLLV
CALL PRINT1(NOV,NA,NRA,NCA,NCB,NRC,
* MA1,MA2,A,MB1,MB2,B,MC1,C,MD1,D,MFA1,FA,FC)
ENDIF
C
C Eliminate local variables by "Condensation of Parameters".
C
IF(IFST.NE.0)THEN
NDECOM=NDECOM+1
CALL CONPAR(NOV,NA,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* NRC,MC1,C,MD1,D,WKDR,IWKDR)
C PRINT DEBUG OUTPUT
IF(IDB.GE.2)CALL PRINT1(NOV,NA,NRA,NCA,NCB,NRC,
* MA1,MA2,A,MB1,MB2,B,MC1,C,MD1,D,MFA1,FA,FC)
ENDIF
C
C Allocate workspace.
C
LA1=1
LA2=LA1+NOV**2*NA
LB=LA2+NOV**2*NA
LC=LB+NOV*NCB*NA
LFA=LC+NRC*NOV*(NA+1)
LS=LFA+NOV*(NA+2)
LE=LS+NOV**2*NA
LRHSE=LE+(NRC+NOV)**2
LXE=LRHSE+NRC+NOV
LT=LXE+NOV+NRC
LNEXT=LT+NOV**2*NA
C
LIR=1
LIC=LIR+2*NOV*(NA-1)
LLC=LIC+NOV*(NA-1)
LNEXT=LLC+NA
C
NAM1=NA-1
NAP1=NA+1
NAP2=NA+2
NXE=NOV+NRC
C
C Copy the reduced system generated by CONPAR into the workspace area
C for further processing by REDUCE.
C
IF(IFST.NE.0)THEN
CALL COPYCP(NA,NOV,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* NRC,MC1,C,WKDR(LA1),WKDR(LA2),WKDR(LB),WKDR(LC))
C Reduction of the system.
CALL REDUCE(NA,NOV,NOV2,NCB,NRC,MD1,D,WKDR(LA1),WKDR(LA2),
* WKDR(LB),WKDR(LC),WKDR(LS),WKDR(LT),IWKDR(LIR),IWKDR(LIC))
ENDIF
C
C Condensation of the right hand side following CONPAR.
C
IF(NLLV.EQ.0)THEN
CALL CONRHS(NOV,NA,NRA,NCA,MA1,MA2,A,NRC,MC1,C,MFA1,
* FA,FC,IWKDR(LLC))
C Copy the reduced right hand side into workspace.
CALL CPYRHS(NA,NRA,NOV,MFA1,FA,WKDR(LFA))
C Print debug output.
IF(IDB.GE.2)CALL PRINT3(NA,NOV,WKDR(LFA))
C Reduction of the right hand side following REDUCE.
CALL REDRHS(NA,NOV,NOV2,NRC,WKDR(LA1),WKDR(LA2),WKDR(LC),
* WKDR(LFA),FC,IWKDR(LIR))
ELSE
C Set right hand sides to zero.
DO 2 I=1,NA
DO 1 J=1,NRA
FA(I,J)=ZERO
1 CONTINUE
2 CONTINUE
DO 3 I=1,NRC
FC(I)=ZERO
3 CONTINUE
DO 4 I=LFA,LS-1
WKDR(I)=ZERO
4 CONTINUE
ENDIF
C
C Print debug output
C
IF(IDB.GE.2)THEN
CALL PRINT2(IDB,NA,NOV,NCB,NRC,WKDR(LS),WKDR(LA1),
* WKDR(LA2),WKDR(LT),WKDR(LB),WKDR(LC),WKDR(LFA))
CALL PRINT3(NA,NOV,WKDR(LFA))
ENDIF
C
NE=NRC+NOV
C
C Solve the system generated by REDUCE
C by Gauss elimination with complete pivoting.
C
CALL DIMRGE(NA,NOV,NCA,WKDR(LS),WKDR(LA2),NCB,WKDR(LB),
* WKDR(LC),NRC,MC1,C,MD1,D,WKDR(LFA),FC,
* NE,WKDR(LE),WKDR(LRHSE),WKDR(LXE),IDB,IR,IC,NLLV)
C
C Backsubstitution in the reduction process.
C
CALL BCKSUB(NOV,NOV2,NA,NCB,WKDR(LS),WKDR(LA1),WKDR(LA2),
* WKDR(LT),WKDR(LB),WKDR(LFA),WKDR(LXE),FC,IWKDR(LIR),
* IWKDR(LIC))
C
C Backsubstitution in the condensation of parameters process.
C
NBCKSB=NBCKSB+1
CALL INFPAR(NA,NOV,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* MFA1,FA,NDRHS,FC,NAP2,WKDR(LFA))
C
101 FORMAT(' SUBROUTINE BRBD : IFST=',I1,1X,'NLLV=',I2)
C
RETURN
END
C
C ------- -------- ------
LOGICAL FUNCTION ERBRBD(NA,NRA,NCA,NCB,NRC)
C
C Checks correctness of dimensions.
C
ERBRBD=.FALSE.
C
IN=NA*(NCA-NRA)+NCB-NRC
ID=NA-1
IF(MOD(IN,ID).NE.0)THEN
WRITE(9,101)
ERBRBD=.TRUE.
RETURN
ENDIF
C
NOV=IN/ID
NEX=NCA-2*NOV
IF(NEX.GE.0)GOTO 1
C
WRITE(9,101)
ERBRBD=.TRUE.
RETURN
101 FORMAT(' ERROR IN DATA IN SUBROUTINE -BRBD-',/,
* ' (LINEAR EQUATION SOLVER)')
C
1 RETURN
END
C
C ---------- ------
SUBROUTINE COPYCP(NA,NOV,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* NRC,MC1,C,A1,A2,BC,CC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBND/ NAM1,NAP1,NXE
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
C
DIMENSION A(MA1,MA2,NCA),B(MB1,MB2,NCB),C(MC1,NDCC)
DIMENSION A1(NOV,NOV,NA),A2(NOV,NOV,NA)
DIMENSION BC(NOV,NCB,NA),CC(NRC,NOV,NAP1)
C
C Copies the condensed system generated by CONPAR into the workspace.
C
DO 4 I=1,NA
C
DO 3 IR=1,NOV
DO 1 IC=1,NOV
A1(IR,IC,I)=A(I,NRA-NOV+IR,IC)
A2(IR,IC,I)=A(I,NRA-NOV+IR,NCA-NOV+IC)
1 CONTINUE
C
DO 2 IC=1,NCB
BC(IR,IC,I)=B(I,NRA-NOV+IR,IC)
2 CONTINUE
C
3 CONTINUE
C
4 CONTINUE
C
DO 7 I=1,NAP1
DO 6 IR=1,NRC
DO 5 IC=1,NOV
CC(IR,IC,I)=C((I-1)*(NCA-NOV)+IC,IR)
5 CONTINUE
6 CONTINUE
7 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE REDUCE(NA,NOV,NOV2,NCB,NRC,MD1,D,A1,A2,B,C,S,T,
* IPR,IPC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBND/ NAM1,NAP1,NXE
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLTIM/ TSETUB,TCONPA,TCONRH,TINFPA,TREDUC,TWR8
C
DIMENSION A1(NOV,NOV,NA),A2(NOV,NOV,NA),S(NOV,NOV,NA)
DIMENSION T(NOV,NOV,NA),IPR(NOV2,NAM1),IPC(NOV,NAM1)
DIMENSION B(NOV,NCB,NA),C(NRC,NOV,NAP1),D(MD1,NCB)
C
C Set initial time (For timing of this subroutine).
C
CALL AUTIM0(TIME0)
C
NOVM1=NOV-1
NAM1=NA-1
C
C CLEAR S AND T
C
DO 3 I1=1,NA
DO 2 IR=1,NOV
DO 1 IC=1,NOV
T(IR,IC,I1)=ZERO
S(IR,IC,I1)=ZERO
1 CONTINUE
2 CONTINUE
3 CONTINUE
C
C Equate first block of S with first block of A1.
C
DO 5 IR=1,NOV
DO 4 IC=1,NOV
S(IR,IC,1)=A1(IR,IC,1)
4 CONTINUE
5 CONTINUE
C
C Initialize pivot arrays.
C
DO 7 I1=1,NAM1
DO 6 I=1,NOV
IPC(I,I1)=I
6 CONTINUE
7 CONTINUE
C
DO 32 I1=1,NAM1
C
I2=I1+1
C
DO 25 IC=1,NOV
ICP1=IC+1
C
C Pivoting :
C
RMXA2=ZERO
DO 9 IR=IC,NOV
DO 8 K=IC,IC
DABA2=DABS(A2(IR,K,I1))
CSGLE DABA2= ABS(A2(IR,K,I1))
IF(DABA2.GT.RMXA2)THEN
IRA2=IR
ICA2=K
RMXA2=DABA2
ENDIF
8 CONTINUE
9 CONTINUE
C
RMXA1=ZERO
DO 11 IR=1,NOV
DO 10 K=IC,IC
DABA1=DABS(A1(IR,K,I2))
CSGLE DABA1= ABS(A1(IR,K,I2))
IF(DABA1.GT.RMXA1)THEN
IRA1=IR
ICA1=K
RMXA1=DABA1
ENDIF
10 CONTINUE
11 CONTINUE
C
IF(RMXA2.GT.RMXA1)THEN
IPR(IC,I1)=IRA2
ITMP=IPC(IC,I1)
IPC(IC,I1)=IPC(ICA2,I1)
IPC(ICA2,I1)=ITMP
DO 12 K=1,NOV
IF(K.GE.IC)THEN
TMP=A2(IC,K,I1)
A2(IC,K,I1)=A2(IRA2,K,I1)
A2(IRA2,K,I1)=TMP
ENDIF
TMP=S(IC,K,I1)
S(IC,K,I1)=S(IRA2,K,I1)
S(IRA2,K,I1)=TMP
TMP=T(IC,K,I1)
T(IC,K,I1)=T(IRA2,K,I1)
T(IRA2,K,I1)=TMP
12 CONTINUE
DO 13 K=1,NCB
TMP=B(IC,K,I1)
B(IC,K,I1)=B(IRA2,K,I1)
B(IRA2,K,I1)=TMP
13 CONTINUE
C
ELSE
C
IPR(IC,I1)=IRA1+NOV
ITMP=IPC(IC,I1)
IPC(IC,I1)=IPC(ICA1,I1)
IPC(ICA1,I1)=ITMP
DO 14 K=1,NOV
IF(K.GE.IC)THEN
TMP=A2(IC,K,I1)
A2(IC,K,I1)=A1(IRA1,K,I2)
A1(IRA1,K,I2)=TMP
ENDIF
TMP=S(IC,K,I1)
S(IC,K,I1)=S(IRA1,K,I2)
S(IRA1,K,I2)=TMP
TMP=T(IC,K,I1)
T(IC,K,I1)=A2(IRA1,K,I2)
A2(IRA1,K,I2)=TMP
14 CONTINUE
DO 15 K=1,NCB
TMP=B(IC,K,I1)
B(IC,K,I1)=B(IRA1,K,I2)
B(IRA1,K,I2)=TMP
15 CONTINUE
C
C
ENDIF
C
C End of pivoting.
C
IF(IC.NE.NOV)THEN
DO 19 IR=ICP1,NOV
RM=A2(IR,IC,I1)/A2(IC,IC,I1)
A2(IR,IC,I1)=RM
DO 16 K=ICP1,NOV
A2(IR,K,I1)=A2(IR,K,I1)-RM*A2(IC,K,I1)
16 CONTINUE
DO 17 K=1,NOV
S(IR,K,I1)=S(IR,K,I1)-RM*S(IC,K,I1)
T(IR,K,I1)=T(IR,K,I1)-RM*T(IC,K,I1)
17 CONTINUE
DO 18 K=1,NCB
B(IR,K,I1)=B(IR,K,I1)-RM*B(IC,K,I1)
18 CONTINUE
19 CONTINUE
ENDIF
C
DO 24 IR=1,NOV
RM=A1(IR,IC,I2)/A2(IC,IC,I1)
A1(IR,IC,I2)=RM
IF(ICP1.GT.NOV)GOTO 21
DO 20 K=ICP1,NOV
A1(IR,K,I2)=A1(IR,K,I2)-RM*A2(IC,K,I1)
20 CONTINUE
21 DO 22 K=1,NOV
S(IR,K,I2)=S(IR,K,I2)-RM*S(IC,K,I1)
A2(IR,K,I2)=A2(IR,K,I2)-RM*T(IC,K,I1)
22 CONTINUE
DO 23 K=1,NCB
B(IR,K,I2)=B(IR,K,I2)-RM*B(IC,K,I1)
23 CONTINUE
24 CONTINUE
C
25 CONTINUE
C
DO 31 IC=1,NOV
ICP1=IC+1
DO 30 IR=1,NRC
RM=C(IR,IC,I2)/A2(IC,IC,I1)
C(IR,IC,I2)=RM
IF(ICP1.GT.NOV)GOTO 27
DO 26 K=ICP1,NOV
C(IR,K,I2)=C(IR,K,I2)-RM*A2(IC,K,I1)
26 CONTINUE
27 DO 28 K=1,NOV
C(IR,K,1)=C(IR,K,1)-RM*S(IC,K,I1)
C(IR,K,I2+1)=C(IR,K,I2+1)-RM*T(IC,K,I1)
28 CONTINUE
DO 29 K=1,NCB
D(IR,K)=D(IR,K)-RM*B(IC,K,I1)
29 CONTINUE
30 CONTINUE
31 CONTINUE
C
32 CONTINUE
C
C Determine the time spent in this subroutine.
C
CALL AUTIM1(TIME1)
TREDUC=TREDUC + TIME1-TIME0
C
RETURN
END
C
C ---------- ------
SUBROUTINE CPYRHS(NA,NRA,NOV,MFA1,FA,FAC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBND/ NAM1,NAP1,NXE
C
DIMENSION FA(MFA1,NRA),FAC(NOV,NAP1)
C
DO 2 I=1,NA
DO 1 IR=1,NOV
FAC(IR,I)=FA(I,NRA-NOV+IR)
1 CONTINUE
2 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE REDRHS(NA,NOV,NOV2,NRC,A1,A2,C,FA,FC,IPR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBND/ NAM1,NAP1,NXE
C
DIMENSION A1(NOV,NOV,NA),A2(NOV,NOV,NA),C(NRC,NOV,NAP1)
DIMENSION FA(NOV,NAP1),FC(NDRHS)
DIMENSION IPR(NOV2,NAM1)
C
NOVM1=NOV-1
C
DO 6 I1=1,NAM1
I2=I1+1
C
DO 3 IC=1,NOV
ICP1=IC+1
C
C Interchanges.
C
IF(IPR(IC,I1).LE.NOV)THEN
TMP=FA(IC,I1)
FA(IC,I1)=FA(IPR(IC,I1),I1)
FA(IPR(IC,I1),I1)=TMP
ELSE
TMP=FA(IC,I1)
FA(IC,I1)=FA(IPR(IC,I1)-NOV,I2)
FA(IPR(IC,I1)-NOV,I2)=TMP
ENDIF
C
C End interchanges.
C
IF(IC.NE.NOV)THEN
DO 1 IR=ICP1,NOV
RM=A2(IR,IC,I1)
FA(IR,I1)=FA(IR,I1)-RM*FA(IC,I1)
1 CONTINUE
ENDIF
C
DO 2 IR=1,NOV
RM=A1(IR,IC,I2)
FA(IR,I2)=FA(IR,I2)-RM*FA(IC,I1)
2 CONTINUE
3 CONTINUE
C
DO 5 IC=1,NOV
DO 4 IR=1,NRC
RM=C(IR,IC,I2)
FC(IR)=FC(IR)-RM*FA(IC,I1)
4 CONTINUE
5 CONTINUE
C
6 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE DIMRGE(NA,NOV,NCA,S,A2,NCB,B,CC,NRC,MC1,C,
* MD1,D,FA,FC,NE,E,RHSE,XE,IDB,IR,IC,NLLV)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBND/ NAM1,NAP1,NXE
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION S(NOV,NOV,NA),A2(NOV,NOV,NA)
DIMENSION B(NOV,NCB,NA),CC(NRC,NOV,NAP1),FA(NOV,NAP1)
DIMENSION IR(NDIRC),IC(NDIRC)
DIMENSION C(MC1,NDCC),D(MD1,NCB),FC(NDRHS)
DIMENSION E(NE,NE),RHSE(NE),XE(NE)
C
DO 3 I=1,NOV
C
DO 1 J=1,NOV
E(I,J)=S(I,J,NA)
E(I,NOV+J)=A2(I,J,NA)
1 CONTINUE
C
DO 2 J=1,NCB
E(I,2*NOV+J)=B(I,J,NA)
2 CONTINUE
C
RHSE(I)=FA(I,NA)
C
3 CONTINUE
C
DO 6 I=1,NRC
C
DO 4 J=1,NOV
E(NOV+I,J)=CC(I,J,1)
E(NOV+I,NOV+J)=CC(I,J,NA+1)
4 CONTINUE
C
DO 5 J=1,NCB
E(NOV+I,2*NOV+J)=D(I,J)
5 CONTINUE
C
RHSE(NOV+I)=FC(I)
C
6 CONTINUE
C
IF(IDB.NE.0)THEN
WRITE(9,101)
DO 7 I=1,NE
WRITE(9,102)(E(I,J),J=1,NE),RHSE(I)
7 CONTINUE
ENDIF
C
IF(NLLV.EQ.0)THEN
CALL GE(NE,NE,E,1,NE,XE,NE,RHSE,IR,IC)
ELSEIF(NLLV.GT.0)THEN
CALL NLVC(NE,NE,IABS(NLLV),E,XE,IR,IC)
ELSE
DO 8 I=1,NE-1
RHSE(I)=ZERO
8 CONTINUE
RHSE(NE)=ONE
CALL GE(NE,NE,E,1,NE,XE,NE,RHSE,IR,IC)
ENDIF
C
IF(IDB.NE.0)WRITE(9,103)(XE(I),I=1,NE)
C
DO 9 I=1,NRC
FC(I)=XE(NOV+I)
9 CONTINUE
C
101 FORMAT(' REDUCED SYSTEM:')
102 FORMAT(1X,11E10.3)
103 FORMAT(' SOLUTION VECTOR :',/,10E10.3)
C
RETURN
END
C
C ---------- ----
SUBROUTINE BCKSUB(NOV,NOV2,NA,NCB,S,A1,A2,T,B,FA,XE,FC,IPR,IPC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBND/ NAM1,NAP1,NXE
C
C Backsubstitution in the reduction process.
C
DIMENSION S(NOV,NOV,NA),A1(NOV,NOV,NA),A2(NOV,NOV,NA)
DIMENSION T(NOV,NOV,NA),B(NOV,NCB,NA)
DIMENSION FA(NOV,NAP1),XE(NXE),FC(NDRHS)
DIMENSION IPR(NOV2,NAM1),IPC(NOV,NAM1)
C
DO 1 I=1,NOV
FA(I,NA+1)=XE(NOV+I)
1 CONTINUE
C
DO 8 I=1,NAM1
I1=NA-I
I2=I1+1
DO 7 J1=1,NOV
J=NOV+1-J1
FA(J,I2)=FA(J,I1)
DO 3 K=1,NOV
FA(J,I2)=FA(J,I2)-S(J,K,I1)*XE(K)
FA(J,I2)=FA(J,I2)-T(J,K,I1)*FA(K,I2+1)
3 CONTINUE
DO 4 K=1,NCB
FA(J,I2)=FA(J,I2)-B(J,K,I1)*FC(NOV+K)
4 CONTINUE
IF(J.EQ.NOV)GOTO 6
K1=J+1
DO 5 K=K1,NOV
FA(J,I2)=FA(J,I2)-A2(J,K,I1)*FA(K,I2)
5 CONTINUE
6 FA(J,I2)=FA(J,I2)/A2(J,J,I1)
7 CONTINUE
8 CONTINUE
C
DO 9 K=1,NOV
FA(K,1)=XE(K)
9 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE PRINT1(NOV,NA,NRA,NCA,NCB,NRC,
* MA1,MA2,A,MB1,MB2,B,MC1,C,MD1,D,MFA1,FA,FC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
C
DIMENSION A(MA1,MA2,NCA),B(MB1,MB2,NCB),C(MC1,NDCC),D(MD1,NCB)
DIMENSION FA(MFA1,NRA),FC(NDRHS)
C
WRITE(9,101)
DO 2 I=1,NA
WRITE(9,102)I
DO 1 IR=1,NRA
WRITE(9,103)(A(I,IR,IC),IC=1,NCA),(B(I,IR,IC),IC=1,NCB)
* ,FA(I,IR)
1 CONTINUE
2 CONTINUE
C
WRITE(9,104)
C
DO 4 I=1,NA
WRITE(9,102)I
DO 3 IR=1,NRC
IC1=(I-1)*(NCA-NOV)+1
IC2=IC1+NCA-NOV-1
WRITE(9,103)(C(IC,IR),IC=IC1,IC2)
3 CONTINUE
4 CONTINUE
C
WRITE(9,106)
IC1=NA*(NCA-NOV)+1
IC2=IC1+NOV-1
DO 5 IR=1,NRC
WRITE(9,103)(C(IC,IR),IC=IC1,IC2)
5 CONTINUE
C
WRITE(9,105)
C
DO 6 IR=1,NRC
WRITE(9,103)(D(IR,IC),IC=1,NCB),FC(IR)
6 CONTINUE
C
101 FORMAT(' A , B , FA (FULL DIMENSION) :')
102 FORMAT(' I=',I2)
103 FORMAT(1X,12E10.3)
104 FORMAT(' C (FULL DIMENSION) :')
105 FORMAT(' D , FC')
106 FORMAT(' LAST NOV COLUMNS OF C :')
C
RETURN
END
C
C ---------- ------
SUBROUTINE PRINT2(IDB,NA,NOV,NCB,NRC,S,A1,A2,T,B,C,F)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLBND/ NAM1,NAP1,NXE
C
DIMENSION A1(NOV,NOV,NA),A2(NOV,NOV,NA),B(NOV,NCB,NA)
DIMENSION S(NOV,NOV,NA),T(NOV,NOV,NA),F(NOV,NAP1)
DIMENSION C(NRC,NOV,NAP1)
C
WRITE(9,101)
C
DO 2 I=1,NA
WRITE(9,104)I
DO 1 IR=1,NOV
WRITE(9,102)(A1(IR,IC,I),IC=1,NOV),(A2(IR,IC,I),IC=1,NOV),
* ( B(IR,IC,I),IC=1,NCB)
1 CONTINUE
2 CONTINUE
C
WRITE(9,105)
DO 4 I=1,NA
WRITE(9,104)I
DO 3 IR=1,NOV
WRITE(9,102)(S(IR,IC,I),IC=1,NOV),(T(IR,IC,I),IC=1,NOV)
3 CONTINUE
4 CONTINUE
C
WRITE(9,103)
C
DO 6 I=1,NA
WRITE(9,104)I
DO 5 IR=1,NRC
WRITE(9,102)(C(IR,IC,I),IC=1,NOV)
5 CONTINUE
6 CONTINUE
C
C The following can be used if a linear system is given
C for which all all solution components are equal to 1.
C Compute residuals for test problem :
C
IF(IDB.EQ.3)THEN
WRITE(9,106)
DO 10 I1=1,NA
DO 9 IR=1,NOV
R=F(IR,I1)
DO 7 IC=1,NOV
R=R-S(IR,IC,I1)
R=R-T(IR,IC,I1)
IF(IC.GE.IR)R=R-A2(IR,IC,I1)
7 CONTINUE
DO 8 IC=1,NCB
R=R-B(IR,IC,I1)
8 CONTINUE
WRITE(9,107)R
9 CONTINUE
10 CONTINUE
ENDIF
C
101 FORMAT(' A1 , A2 , B :')
102 FORMAT(1X,9E10.3)
103 FORMAT(' C :')
104 FORMAT(' I=',I2)
105 FORMAT(' S AND T : ')
106 FORMAT(' RESIDUALS IN PRINT2')
107 FORMAT(1X,E10.3)
C
RETURN
END
C
C --------- ------
SUBROUTINE PRINT3(NA,NOV,FA)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
DIMENSION FA(NOV,NA)
C
WRITE(9,101)
DO 1 I=1,NA
WRITE(9,102)I
WRITE(9,103)(FA(IR,I),IR=1,NOV)
1 CONTINUE
C
101 FORMAT(' FA : ')
102 FORMAT(' I=',I3)
103 FORMAT(1X,9E10.3)
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Vectorized Subroutines for the Linear Equation Solver BRBD
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ------
SUBROUTINE CONPAR(NOV,NA,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* NRC,MC1,C,MD1,D,RM,LC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLTIM/ TSETUB,TCONPA,TCONRH,TINFPA,TREDUC,TWR8
C
DIMENSION A(MA1,MA2,NCA),B(MB1,MB2,NCB),C(MC1,NRC)
* ,D(MD1,NCB),RM(NA),LC(NA)
C
C Set initial time (for timing of this subroutine).
C
CALL AUTIM0(TIME0)
C
NEX=NCA-2*NOV
IF(NEX.EQ.0)RETURN
C
C Condensation of parameters ( Elimination of "local" variables ).
C
M1=NOV+1
M2=NOV+NEX
C
DO 9 IC=M1,M2
IR1=IC-NOV+1
IRP=IR1-1
ICP1=IC+1
C
DO 4 IR=IR1,NRA
DO 10 I=1,NA
RM(I)=A(I,IR,IC)/A(I,IRP,IC)
A(I,IR,IC)=RM(I)
10 CONTINUE
DO 1 L=1,NOV
DO 11 I=1,NA
A(I,IR,L)=A(I,IR,L)-RM(I)*A(I,IRP,L)
11 CONTINUE
1 CONTINUE
DO 2 L=ICP1,NCA
DO 12 I=1,NA
A(I,IR,L)=A(I,IR,L)-RM(I)*A(I,IRP,L)
12 CONTINUE
2 CONTINUE
DO 3 L=1,NCB
DO 13 I=1,NA
B(I,IR,L)=B(I,IR,L)-RM(I)*B(I,IRP,L)
13 CONTINUE
3 CONTINUE
4 CONTINUE
C
DO 8 IR=1,NRC
DO 40 I=1,NA
LC(I)=(I-1)*(NCA-NOV)+IC
40 CONTINUE
DO 50 I=1,NA
RM(I)=C(LC(I),IR)/A(I,IRP,IC)
C(LC(I),IR)=RM(I)
50 CONTINUE
DO 5 L=1,NOV
DO 51 I=1,NA
LC(I)=(I-1)*(NCA-NOV)+L
51 CONTINUE
DO 52 I=1,NA
C(LC(I),IR)=C(LC(I),IR)-RM(I)*A(I,IRP,L)
52 CONTINUE
5 CONTINUE
DO 6 L=ICP1,NCA
DO 60 I=1,NA
LC(I)=(I-1)*(NCA-NOV)+L
60 CONTINUE
DO 61 I=1,NA
C(LC(I),IR)=C(LC(I),IR)-RM(I)*A(I,IRP,L)
61 CONTINUE
6 CONTINUE
DO 7 L=1,NCB
DO 71 I=1,NA
D(IR,L)=D(IR,L)-RM(I)*B(I,IRP,L)
71 CONTINUE
7 CONTINUE
8 CONTINUE
C
9 CONTINUE
C
C Determine the time spent in this subroutine.
C
CALL AUTIM1(TIME1)
TCONPA=TCONPA + TIME1-TIME0
C
RETURN
END
C
C ---------- ------
SUBROUTINE CONRHS(NOV,NA,NRA,NCA,MA1,MA2,A,NRC,MC1,C,
* MFA1,FA,FC,LC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLTIM/ TSETUB,TCONPA,TCONRH,TINFPA,TREDUC,TWR8
C
DIMENSION A(MA1,MA2,NCA),C(MC1,NDCC),FA(MFA1,NRA),FC(NDRHS),LC(NA)
C
C Set initial time (for timing of this subroutine).
C
CALL AUTIM0(TIME0)
C
C Condensation of the right hand side.
C
NEX=NCA-2*NOV
IF(NEX.EQ.0)RETURN
C
M1=NOV+1
M2=NOV+NEX
C
DO 3 IC=M1,M2
IR1=IC-NOV+1
IRP=IR1-1
C
DO 1 IR=IR1,NRA
C Note that RM=A(I,IR,IC) is the multiplier.
DO 11 I=1,NA
FA(I,IR)=FA(I,IR)-A(I,IR,IC)*FA(I,IRP)
11 CONTINUE
1 CONTINUE
C
DO 2 IR=1,NRC
DO 20 I=1,NA
LC(I)=(I-1)*(NCA-NOV)+IC
20 CONTINUE
C Note that RM=C(LC(I),IR) is the multiplier.
DO 21 I=1,NA
FC(IR)=FC(IR)-C(LC(I),IR)*FA(I,IRP)
21 CONTINUE
2 CONTINUE
C
3 CONTINUE
C
C Determine the time spent in this subroutine,
C
CALL AUTIM1(TIME1)
TCONRH=TCONRH + TIME1-TIME0
C
RETURN
END
C
C ---------- ------
SUBROUTINE INFPAR(NA,NOV,NRA,NCA,MA1,MA2,A,NCB,MB1,MB2,B,
* MFA1,FA,NC,FC,NFADR,FADR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLDIM/ NDIMP1,NDIRC,NTSTP1,NDCC,NDRHS,NDBC,NUICD,NDICD
* ,NWBR,NIWBR
COMMON /BLBND/ NAM1,NAP1,NXE
COMMON /BLTIM/ TSETUB,TCONPA,TCONRH,TINFPA,TREDUC,TWR8
C
DIMENSION A(MA1,MA2,NCA),B(MB1,MB2,NCB),FA(MFA1,NRA),FC(NC)
DIMENSION FADR(NOV,NFADR)
C
C Set initial time (for timing of this subroutine).
C
CALL AUTIM0(TIME0)
C
C Determine the "local" variables by backsubstitution.
C
NRAM=NRA-NOV
DO 9 IR1=1,NRAM
IR=NRAM+1-IR1
DO 1 K=1,NOV
DO 11 I=1,NA
FA(I,IR)=FA(I,IR)-A(I,IR,K)*FADR(K,I)
11 CONTINUE
1 CONTINUE
DO 2 K=1,NOV
FADR(K,NA+2)=FC(K)
2 CONTINUE
DO 3 K=1,NOV
DO 31 I=1,NA
FA(I,IR)=FA(I,IR)-A(I,IR,NCA-NOV+K)*FADR(K,I+1)
31 CONTINUE
3 CONTINUE
DO 6 K=1,NCB
DO 61 I=1,NA
FA(I,IR)=FA(I,IR)-B(I,IR,K)*FC(NOV+K)
61 CONTINUE
6 CONTINUE
IF(IR1.EQ.1)GOTO 8
K1=NCA-NOV-IR1+2
K2=NCA-NOV
DO 7 K=K1,K2
DO 71 I=1,NA
FA(I,IR)=FA(I,IR)-A(I,IR,K)*FA(I,K)
71 CONTINUE
7 CONTINUE
8 CONTINUE
DO 91 I=1,NA
FA(I,IR+NOV)=FA(I,IR)/A(I,IR,NCA-NOV-IR1+1)
91 CONTINUE
9 CONTINUE
C
C Copy the solution generated by DRBSUB (stored in FADR)
C into the final solution vector (stored in FA).
C
DO 22 J=1,NOV
DO 21 I=1,NA
FA(I,J)=FADR(J,I)
21 CONTINUE
22 CONTINUE
C
C Determine the time spent in this subroutine.
C
CALL AUTIM1(TIME1)
TINFPA=TINFPA + TIME1-TIME0
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Limit Points (Algebraic Problems)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNLP(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-par continuation of limit points.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1),ICP(1)
C
C Generate the function.
C
CALL FFLP(NDIM,U,UOLD,ICP,PAR,F,NDM,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFLP(NDIM,U1XX,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
CALL FFLP(NDIM,U2XX,UOLD,ICP,PAR,F2XX,NDM,DFUXX,DFPXX)
DO 3 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))+EP
C
CALL FFLP(NDIM,U,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
C
DO 5 J=1,NDIM
DFDP(J,ICP(1))=(F1XX(J)-F(J))/EP
5 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))-EP
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFLP(NDIM1,U,UOLD,ICP1,PAR1,F,NDM,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
C
DIMENSION U(NDIM)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
C
IJAC=1
PAR(ICP(2))=U(NDIM)
IF(IPS.EQ.-1) THEN
CALL FNDS(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
ELSE
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
ENDIF
C
DO 2 I=1,NDM
F(NDM+I)=ZERO
DO 1 J=1,NDM
F(NDM+I)=F(NDM+I)+DFDU(I,J)*U(NDM+J)
1 CONTINUE
2 CONTINUE
C
F(NDIM)=-ONE
CSGLE F(NDIM)=-ONE
C
DO 3 I=1,NDM
F(NDIM)=F(NDIM)+U(NDM+I)*U(NDM+I)
3 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNLP(IBR,U,NDM2,SMAT,DFDU,DFUXX,DFDP,V,F,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
LOGICAL FOUND
C
C Generates starting data for the continuation of limit points.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
C
DIMENSION U(NDIM),V(NDMP1),F(NDM),SMAT(NDM2,NDM2)
DIMENSION DFDU(NDM,NDM),DFUXX(NDMP1,NDMP1),DFDP(NDM,1)
DIMENSION IR(1),IC(1)
C
CALL FINDL3(IRS,ITP,NFPAR1,FOUND)
CALL READL3(IPS,IBR,U,PAR)
C
IJAC=1
IF(IPS.EQ.-1)THEN
CALL FNDS(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
ELSE
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
ENDIF
CALL NLVC(NDM,NDM,1,DFDU,V,IR,IC)
CALL NRMLZ(NDM,V)
DO 1 I=1,NDM
U(NDM+I)=V(I)
1 CONTINUE
U(NDIM)=PAR(ICP(2))
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Optimization of Algebraic Systems
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNC1(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generate the equations for the continuation scheme used for
C the optimization of algebraic systems (one parameter).
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,NPAR)
C
PAR(ICP(2))=U(NDIM)
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
C Rearrange (Since dimensions in FNC1 and FUNI differ).
C
IF(IJAC.NE.0)THEN
DO 2 J=NDM,1,-1
DO 1 I=NDM,1,-1
DFDU(I,J)=DFDU( (J-1)*NDM+I ,1 )
1 CONTINUE
2 CONTINUE
C
DO 4 J=NPAR,1,-1
DO 3 I=NDM,1,-1
DFDP(I,J)=DFDP( (J-1)*NDM+I , 1 )
3 CONTINUE
4 CONTINUE
ENDIF
C
CALL FOPI(NDM,U,ICP,PAR,IJAC,F(NDIM),DDUXX,DDPXX)
F(NDIM)=PAR(ICP(1))-F(NDIM)
C
IF(IJAC.NE.0)THEN
DO 5 I=1,NDM
DFDU(NDIM,I)=-DDUXX(I)
DFDU(I,NDIM)=DFDP(I,ICP(2))
DFDP(I,ICP(1))=0
5 CONTINUE
DFDU(NDIM,NDIM)=-DDPXX(ICP(2))
DFDP(NDIM,ICP(1))=1
ENDIF
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNC1(IBR,U,NDM2,SMAT,DFDU,DFUXX,DFDP,V,F,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
LOGICAL FOUND
C
C Generate starting data for optimization problems (one parameter).
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
C
DIMENSION U(NDIM),V(NDMP1),F(NDM),SMAT(NDM2,NDM2)
DIMENSION DFDU(NDM,NDM),DFUXX(NDMP1,NDMP1),DFDP(NDM,1)
DIMENSION IR(1),IC(1)
C
CALL STPNT(NDIM,U,PAR)
CALL FOPI(NDM,U,ICP,PAR,0,FOP,DFDU,DFDP)
PAR(ICP(1))=FOP
U(NDIM)=PAR(ICP(2))
C
RETURN
END
C
C ---------- ----
SUBROUTINE FNC2(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generate the equations for the continuation scheme used for the
C optimization of algebraic systems (more than one parameter).
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFC2(NDIM,U,UOLD,ICP,PAR,F,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFC2(NDIM,U1XX,UOLD,ICP,PAR,F1XX,DFUXX,DFPXX)
CALL FFC2(NDIM,U2XX,UOLD,ICP,PAR,F2XX,DFUXX,DFPXX)
DO 3 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 5 I=1,NDIM
DFDP(I,ICP(1))=0
5 CONTINUE
DFDP(NDIM,ICP(1))=1
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFC2(NDIM,U,UOLD,ICP,PAR,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
C
IJAC=1
DO 1 I=2,NFPAR
PAR(ICP(I))=U(2*NDM+I)
1 CONTINUE
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
CALL FOPI(NDM,U,ICP,PAR,IJAC,FOP,DDUXX,DDPXX)
C
DO 3 I=1,NDM
F(NDM+I)=DDUXX(I)*U(2*NDM+1)
DO 2 J=1,NDM
F(NDM+I)=F(NDM+I)+DFDU(J,I)*U(NDM+J)
2 CONTINUE
3 CONTINUE
C
NDM2=2*NDM
ICPM=NFPAR-2
DO 4 I=1,ICPM
F(NDM2+I)=DDPXX(ICP(I+1))*U(NDM2+1)
4 CONTINUE
C
DO 6 I=1,ICPM
DO 5 J=1,NDM
F(NDM2+I)=F(NDM2+I)+U(NDM+J)*DFDP(J,ICP(I+1))
5 CONTINUE
6 CONTINUE
C
F(NDIM-1)=U(NDM2+1)*U(NDM2+1)-1
DO 7 J=1,NDM
F(NDIM-1)=F(NDIM-1)+U(NDM+J)*U(NDM+J)
7 CONTINUE
F(NDIM)=PAR(ICP(1))-FOP
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNC2(IBR,U,NDM2,SMAT,DFDU,DFU,DFDP,V,F,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
LOGICAL FOUND
C
C Generates starting data for the continuation equations for
C optimization of algebraic systems (More than one parameter).
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),V(NDMP1),F(NDM),SMAT(NDM2,NDM2)
DIMENSION DFDU(NDM,NDM),DFU(NDMP1,NDMP1),DFDP(NDM,1)
DIMENSION IR(1),IC(1)
C
CALL FINDL3(IRS,ITP,NFPAR,FOUND)
NFPAR=NFPAR+1
CALL READL3(IPS,IBR,U,PAR)
C
IF(NFPAR.EQ.3)THEN
IJAC=1
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
CALL FOPI(NDM,U,ICP,PAR,IJAC,FOP,DDUXX,DDPXX)
C TRANSPOSE
DO 2 I=1,NDM
DO 1 J=1,NDM
DFU(I,J)=DFDU(J,I)
1 CONTINUE
2 CONTINUE
DO 3 I=1,NDM
DFU(I,NDMP1)=DDUXX(I)
DFU(NDMP1,I)=DFDP(I,ICP(2))
3 CONTINUE
DFU(NDMP1,NDMP1)=DDPXX(ICP(2))
CALL NLVC(NDMP1,NDMP1,1,DFU,V,IR,IC)
CALL NRMLZ(NDMP1,V)
DO 4 I=1,NDMP1
U(NDM+I)=V(I)
4 CONTINUE
PAR(ICP(1))=FOP
ENDIF
C
DO 5 I=1,NFPAR-1
U(NDIM-NFPAR+1+I)=PAR(ICP(I+1))
5 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for Discrete Dynamical Systems
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNDS(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generate the equations for continuing fixed points.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
CALL FUNI(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 1 I=1,NDIM
F(I)=F(I)-U(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
DO 2 I=1,NDIM
DFDU(I,I)=DFDU(I,I)-ONE
2 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Hopf Bifurcation Points (Maps)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNHD(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-parameter continuation of Hopf
C bifurcation points for maps.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),UOLD(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFHD(NDIM,U,UOLD,ICP,PAR,F,NDM,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 2 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
2 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 5 I=1,NDIM
DO 3 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
3 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFHD(NDIM,U1XX,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
CALL FFHD(NDIM,U2XX,UOLD,ICP,PAR,F2XX,NDM,DFUXX,DFPXX)
DO 4 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
4 CONTINUE
5 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))+EP
C
CALL FFHD(NDIM,U,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
C
DO 6 J=1,NDIM
DFDP(J,ICP(1))=(F1XX(J)-F(J))/EP
6 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))-EP
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFHD(NDIM,U,UOLD,ICP,PAR,F,NDM,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
DIMENSION UOLD(NDIM)
C
NDM2=2*NDM
C
THETA=U(NDIM-1)
S1=DSIN(THETA)
CSGLE S1= SIN(THETA)
C1=DCOS(THETA)
CSGLE C1= COS(THETA)
PAR(ICP(2))=U(NDIM)
IJAC=1
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
DO 1 I=1,NDM
F(I)=F(I)-U(I)
DFDU(I,I)=DFDU(I,I)-C1
1 CONTINUE
C
DO 3 I=1,NDM
F(NDM+I)=S1*U(NDM2+I)
F(NDM2+I)=-S1*U(NDM+I)
DO 2 J=1,NDM
F(NDM+I)=F(NDM+I)+DFDU(I,J)*U(NDM+J)
F(NDM2+I)=F(NDM2+I)+DFDU(I,J)*U(NDM2+J)
2 CONTINUE
3 CONTINUE
C
F(NDIM-1)=-ONE
C
DO 4 I=1,NDM
F(NDIM-1)=F(NDIM-1)+U(NDM+I)*U(NDM+I)+U(NDM2+I)*U(NDM2+I)
4 CONTINUE
C
F(NDIM)=ZERO
C
DO 5 I=1,NDM
F(NDIM)=F(NDIM)+UOLD(NDM2+I)*U(NDM+I)-UOLD(NDM+I)*U(NDM2+I)
5 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNHD(IBR,U,NDM2,SMAT,DFDU,DFUXX,DFDP,V,F,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
LOGICAL FOUND
C
C Generates starting data for the continuation of Hopf bifurcation
C points for maps.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),V(NDMP1),F(NDM),SMAT(NDM2,NDM2)
DIMENSION DFDU(NDM,NDM),DFUXX(NDMP1,NDMP1),DFDP(NDM,1)
DIMENSION IR(1),IC(1)
C
CALL FINDL3(IRS,ITP,NFPAR1,FOUND)
CALL READL3(IPS,IBR,U,PAR)
C
THETA=PI(TWO)/PAR(11)
S1=DSIN(THETA)
CSGLE S1= SIN(THETA)
C1=DCOS(THETA)
CSGLE C1= COS(THETA)
IJAC=1
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
NDM2=2*NDM
DO 2 I=1,NDM2
DO 1 J=1,NDM2
SMAT(I,J)=ZERO
1 CONTINUE
2 CONTINUE
C
DO 3 I=1,NDM
SMAT(I,NDM+I)=S1
3 CONTINUE
C
DO 4 I=1,NDM
SMAT(NDM+I,I)=-S1
4 CONTINUE
C
DO 6 I=1,NDM
DO 5 J=1,NDM
SMAT(I,J)=DFDU(I,J)
SMAT(NDM+I,NDM+J)=DFDU(I,J)
5 CONTINUE
SMAT(I,I)=SMAT(I,I)-C1
SMAT(NDM+I,NDM+I)=SMAT(NDM+I,NDM+I)-C1
6 CONTINUE
CALL NLVC(NDM2,NDM2,2,SMAT,V,IR,IC)
CALL NRMLZ(NDM2,V)
C
DO 7 I=1,NDM2
U(NDM+I)=V(I)
7 CONTINUE
C
U(NDIM-1)=THETA
U(NDIM)=PAR(ICP(2))
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Hopf Bifurcation Points (ODE)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNHB(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-parameter continuation of Hopf
C bifurcation points in ODE.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),UOLD(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFHB(NDIM,U,UOLD,ICP,PAR,F,NDM,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 2 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
2 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 5 I=1,NDIM
DO 3 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
3 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFHB(NDIM,U1XX,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
CALL FFHB(NDIM,U2XX,UOLD,ICP,PAR,F2XX,NDM,DFUXX,DFPXX)
DO 4 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
4 CONTINUE
5 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))+EP
C
CALL FFHB(NDIM,U,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
C
DO 6 J=1,NDIM
DFDP(J,ICP(1))=(F1XX(J)-F(J))/EP
6 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))-EP
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFHB(NDIM,U,UOLD,ICP,PAR,F,NDM,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
DIMENSION UOLD(NDIM)
C
NDM2=2*NDM
C
ROM=U(NDIM-1)
PAR(11)=ROM*PI(TWO)
PAR(ICP(2))=U(NDIM)
IJAC=1
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 2 I=1,NDM
F(NDM+I)=U(NDM2+I)
F(NDM2+I)=-U(NDM+I)
DO 1 J=1,NDM
F(NDM+I)=F(NDM+I)+ROM*DFDU(I,J)*U(NDM+J)
F(NDM2+I)=F(NDM2+I)+ROM*DFDU(I,J)*U(NDM2+J)
1 CONTINUE
2 CONTINUE
C
F(NDIM-1)=-ONE
C
DO 3 I=1,NDM
F(NDIM-1)=F(NDIM-1)+U(NDM+I)*U(NDM+I)+U(NDM2+I)*U(NDM2+I)
3 CONTINUE
C
F(NDIM)=ZERO
C
DO 4 I=1,NDM
F(NDIM)=F(NDIM)+UOLD(NDM2+I)*(U(NDM+I)-UOLD(NDM+I)) -
* UOLD(NDM+I)*(U(NDM2+I)-UOLD(NDM2+I))
4 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNHB(IBR,U,NDM2,SMAT,DFDU,DFUXX,DFDP,V,F,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
LOGICAL FOUND
C
C Generates starting data for the 2-parameter continuation of
C Hopf bifurcation point (ODE).
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),V(NDMP1),F(NDM),SMAT(NDM2,NDM2)
DIMENSION DFDU(NDM,NDM),DFUXX(NDMP1,NDMP1),DFDP(NDM,1)
DIMENSION IR(1),IC(1)
C
CALL FINDL3(IRS,ITP,NFPAR1,FOUND)
CALL READL3(IPS,IBR,U,PAR)
C
IJAC=1
RHO=PAR(11)
ROM=RHO/PI(TWO)
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
NDM2=2*NDM
DO 2 I=1,NDM2
DO 1 J=1,NDM2
SMAT(I,J)=ZERO
1 CONTINUE
2 CONTINUE
C
DO 3 I=1,NDM
SMAT(I,NDM+I)=ONE
3 CONTINUE
C
DO 4 I=1,NDM
SMAT(NDM+I,I)=-ONE
4 CONTINUE
C
DO 6 I=1,NDM
DO 5 J=1,NDM
SMAT(I,J)=ROM*DFDU(I,J)
SMAT(NDM+I,NDM+J)=ROM*DFDU(I,J)
5 CONTINUE
6 CONTINUE
CALL NLVC(NDM2,NDM2,2,SMAT,V,IR,IC)
CALL NRMLZ(NDM2,V)
C
DO 7 I=1,NDM2
U(NDM+I)=V(I)
7 CONTINUE
C
U(NDIM-1)=ROM
U(NDIM)=PAR(ICP(2))
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Hopf Bifurcation Points (Waves)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNHW(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-parameter continuation of a
C bifurcation to a traveling wave.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),UOLD(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFHW(NDIM,U,UOLD,ICP,PAR,F,NDM,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 2 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
2 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 5 I=1,NDIM
DO 3 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
3 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFHW(NDIM,U1XX,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
CALL FFHW(NDIM,U2XX,UOLD,ICP,PAR,F2XX,NDM,DFUXX,DFPXX)
DO 4 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
4 CONTINUE
5 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))+EP
C
CALL FFHW(NDIM,U,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
C
DO 6 J=1,NDIM
DFDP(J,ICP(1))=(F1XX(J)-F(J))/EP
6 CONTINUE
C
PAR(ICP(1))=PAR(ICP(1))-EP
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFHW(NDIM,U,UOLD,ICP,PAR,F,NDM,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
DIMENSION UOLD(NDIM)
C
NDM2=2*NDM
C
ROM=U(NDIM-1)
PAR(ICP(2))=U(NDIM)
IJAC=1
CALL FNWS(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 2 I=1,NDM
F(NDM+I)=U(NDM2+I)
F(NDM2+I)=-U(NDM+I)
DO 1 J=1,NDM
F(NDM+I)=F(NDM+I)+ROM*DFDU(I,J)*U(NDM+J)
F(NDM2+I)=F(NDM2+I)+ROM*DFDU(I,J)*U(NDM2+J)
1 CONTINUE
2 CONTINUE
C
F(NDIM-1)=-ONE
C
DO 3 I=1,NDM
F(NDIM-1)=F(NDIM-1)+U(NDM+I)*U(NDM+I)+U(NDM2+I)*U(NDM2+I)
3 CONTINUE
C
F(NDIM)=ZERO
C
DO 4 I=1,NDM
F(NDIM)=F(NDIM)+UOLD(NDM2+I)*(U(NDM+I)-UOLD(NDM+I)) -
* UOLD(NDM+I)*(U(NDM2+I)-UOLD(NDM2+I))
4 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNHW(IBR,U,NDM2,SMAT,DFDU,DFUXX,DFDP,V,F,IR,IC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
LOGICAL FOUND
C
C Generates starting data for the continuation of a bifurcation to a
C traveling wave.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),V(NDMP1),F(NDM),SMAT(NDM2,NDM2)
DIMENSION DFDU(NDM,NDM),DFUXX(NDMP1,NDMP1),DFDP(NDM,1)
DIMENSION IR(1),IC(1)
C
CALL FINDL3(IRS,ITP,NFPAR1,FOUND)
CALL READL3(IPS,IBR,U,PAR)
C
IJAC=1
RHO=PAR(11)
ROM=RHO/PI(TWO)
CALL FNWS(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
NDM2=2*NDM
DO 2 I=1,NDM2
DO 1 J=1,NDM2
SMAT(I,J)=ZERO
1 CONTINUE
2 CONTINUE
C
DO 3 I=1,NDM
SMAT(I,NDM+I)=ONE
3 CONTINUE
C
DO 4 I=1,NDM
SMAT(NDM+I,I)=-ONE
4 CONTINUE
C
DO 6 I=1,NDM
DO 5 J=1,NDM
SMAT(I,J)=ROM*DFDU(I,J)
SMAT(NDM+I,NDM+J)=ROM*DFDU(I,J)
5 CONTINUE
6 CONTINUE
CALL NLVC(NDM2,NDM2,2,SMAT,V,IR,IC)
CALL NRMLZ(NDM2,V)
C
DO 7 I=1,NDM2
U(NDM+I)=V(I)
7 CONTINUE
C
U(NDIM-1)=ROM
U(NDIM)=PAR(ICP(2))
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Periodic Limit Points
C and Period Doubling Bifurcations
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNPL(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Subroutines for the 2-parameter continuation of of limit points
C on branches of periodic solutions.
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFPL(NDIM,U,UOLD,ICP,PAR,F,NDM,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFPL(NDIM,U1XX,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
CALL FFPL(NDIM,U2XX,UOLD,ICP,PAR,F2XX,NDM,DFUXX,DFPXX)
DO 3 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 6 I=1,NFPAR
PAR(ICP(I))=PAR(ICP(I))+EP
C
CALL FFPL(NDIM,U,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
C
DO 5 J=1,NDIM
DFDP(J,ICP(I))=(F1XX(J)-F(J))/EP
5 CONTINUE
C
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFPL(NDIM,U,UOLD,ICP,PAR,F,NDM,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
C
IJAC=1
C1=PAR(ICP(3))
CSGLE C1=PAR(ICP(3))
C2=PAR(ICP(4))
CSGLE C2=PAR(ICP(4))
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 2 I=1,NDM
F(NDM+I)=ZERO
DO 1 J=1,NDM
F(NDM+I)=F(NDM+I)+DFDU(I,J)*U(NDM+J)
1 CONTINUE
F(NDM+I)=C1*F(NDM+I)+C2*F(I)
F(I)=C1*F(I)
2 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE BCPL(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION PAR(1),ICP(1),U0(NBC),U1(NBC),F(NBC)
DIMENSION DBC(NBC,1)
C
DO 1 I=1,NDIM
F(I)=U0(I)-U1(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=2*NDIM+NPAR
DO 3 I=1,NBC
DO 2 J=1,NN
DBC(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDIM
DBC(I,I)=ONE
DBC(I,NDIM+I)=-ONE
4 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE BCPD(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C
C Generate boundary conditions for the 2-parameter continuation
C of period doubling bifurcations.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION PAR(1),ICP(1),U0(NBC),U1(NBC),F(NBC)
DIMENSION DBC(NBC,1)
C
DO 1 I=1,NDM
F(I)=U0(I)-U1(I)
F(NDM+I)=U0(NDM+I)+U1(NDM+I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=2*NDIM+NPAR
DO 3 I=1,NBC
DO 2 J=1,NN
DBC(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDIM
DBC(I,I)=ONE
IF(I.LE.NDM) THEN
DBC(I,NDIM+I)=-ONE
ELSE
DBC(I,NDIM+I)=ONE
ENDIF
4 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE ICPL(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,
* DINT)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(1),UOLD(1),UDOT(1),UPOLD(1),F(3),DINT(NINT,1)
DIMENSION ICP(20),PAR(20)
C
F(1)=ZERO
F(2)=ZERO
F(3)=PAR(ICP(4))**2-ONE
C
DO 1 I=1,NDM
F(1)=F(1)+U(I)*UPOLD(I)
F(2)=F(2)+U(NDM+I)*UPOLD(I)
F(3)=F(3)+U(NDM+I)*U(NDM+I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=NDIM+NPAR
DO 3 I=1,NINT
DO 2 J=1,NN
DINT(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDM
DINT(1,I)=UPOLD(I)
DINT(2,NDM+I)=UPOLD(I)
DINT(3,NDM+I)=TWO*U(NDM+I)
4 CONTINUE
C
DINT(3,NDIM+ICP(4))=TWO*PAR(ICP(4))
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNPL(NTSTRS,NCOLRS,LAB,IBR,M1U,U,UPS,UDOTPS,UPOLDP,
* TM,DTM,NDIM2,SMAT,RNLLV,IR,IC,F,DFDU,DFDP,NODIR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates starting data for the 2-parameter continuation of limit
C points on a branch of periodic solutions.
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),TM(1),DTM(1)
DIMENSION IR(NDIM2),IC(NDIM2),SMAT(NDIM2,NDIM2),RNLLV(NDIM2)
DIMENSION U(NDIM),F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,NPAR)
DIMENSION TEMP(7),ICPRS(20)
C
LOGICAL FOUND
C
CALL FINDL3(IRS,ITP1,NFPAR1,FOUND)
READ(3,*)IBR,NTOT1,ITP1,LAB1,NFPAR1,ISW1,NTPL1,NAR1,NSKIP1,
* NTSTRS,NCOLRS,ICPRS(1),ICPRS(2)
NRSP1=NTSTRS+1
C
DO 2 J=1,NTSTRS
DO 1 I=1,NCOLRS
K1=(I-1)*NDIM+1
K2=K1+NDM-1
READ(3,101)TEMP(I),(UPS(J,K),K=K1,K2)
1 CONTINUE
TM(J)=TEMP(1)
2 CONTINUE
READ(3,101)TM(NRSP1),(UPS(NRSP1,K),K=1,NDM)
C
READ(3,101)RLDOT(1),RLDOT(2)
C
C Read U-dot (derivative with respect to arclength).
C
DO 4 J=1,NTSTRS
DO 3 I=1,NCOLRS
K1=(I-1)*NDIM+NDM+1
K2=I*NDIM
READ(3,101)(UPS(J,K),K=K1,K2)
3 CONTINUE
4 CONTINUE
K1=NDM+1
READ(3,101)(UPS(NRSP1,K),K=K1,NDIM)
C
C Read the parameter values.
C
READ(3,101)(PAR(I),I=1,NPAR)
PAR(ICP(4))=RLDOT(2)
DO 5 I=1,NFPAR
RL(I)=PAR(ICP(I))
5 CONTINUE
C
NODIR=1
C
101 FORMAT(4X,1P7E18.10)
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Limit Points for BVP.
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNBL(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-parameter continuation
C of limit points (BVP).
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),ICP(1),PAR(1),F(NDIM)
DIMENSION DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFBL(NDIM,U,UOLD,ICP,PAR,F,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFBL(NDIM,U1XX,UOLD,ICP,PAR,F1XX,DFUXX,DFPXX)
CALL FFBL(NDIM,U2XX,UOLD,ICP,PAR,F2XX,DFUXX,DFPXX)
DO 3 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 6 I=1,NFPAR
PAR(ICP(I))=PAR(ICP(I))+EP
CALL FFBL(NDIM,U,UOLD,ICP,PAR,F1XX,DFUXX,DFPXX)
DO 5 J=1,NDIM
DFDP(J,ICP(I))=(F1XX(J)-F(J))/EP
5 CONTINUE
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFBL(NDIM,U,UOLD,ICP,PAR,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(1),PAR(1),F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
C
IJAC=1
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
NFPX=NFPAR/2-1
DO 3 I=1,NDM
F(NDM+I)=ZERO
DO 1 J=1,NDM
F(NDM+I)=F(NDM+I)+DFDU(I,J)*U(NDM+J)
1 CONTINUE
IF(NFPX.GT.0)THEN
DO 2 J=1,NFPX
F(NDM+I)=F(NDM+I)
* + DFDP(I,ICP(1+J))*PAR(ICP(NFPAR-NFPX+J))
2 CONTINUE
ENDIF
3 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE BCBL(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the boundary conditions for the 2-parameter continuation
C of limit points (BVP).
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U0(NDIM),U1(NDIM),F(NBC),ICP(1),PAR(1),DBC(NBC,1)
C
C Generate the function.
C
CALL FBBL(NDIM,PAR,ICP,NBC,U0,U1,F,DFUXX)
C
IF(IJAC.EQ.0)RETURN
C
C Derivatives with respect to U0.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U0(I)).GT.UMX)UMX=DABS(U0(I))
CSGLE IF( ABS(U0(I)).GT.UMX)UMX= ABS(U0(I))
1 CONTINUE
EP=HMACH*(ONE+UMX)
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U0(J)
U2XX(J)=U0(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FBBL(NDIM,PAR,ICP,NBC,U1XX,U1,F1XX,DFUXX)
CALL FBBL(NDIM,PAR,ICP,NBC,U2XX,U1,F2XX,DFUXX)
DO 3 J=1,NBC
DBC(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
C Derivatives with respect to U1.
C
UMX=ZERO
DO 5 I=1,NDIM
IF(DABS(U1(I)).GT.UMX)UMX=DABS(U1(I))
CSGLE IF( ABS(U1(I)).GT.UMX)UMX= ABS(U1(I))
5 CONTINUE
EP=HMACH*(ONE+UMX)
DO 8 I=1,NDIM
DO 6 J=1,NDIM
U1XX(J)=U1(J)
U2XX(J)=U1(J)
6 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FBBL(NDIM,PAR,ICP,NBC,U0,U1XX,F1XX,DFUXX)
CALL FBBL(NDIM,PAR,ICP,NBC,U0,U2XX,F2XX,DFUXX)
DO 7 J=1,NBC
DBC(J,NDIM+I)=(F2XX(J)-F1XX(J))/(2*EP)
7 CONTINUE
8 CONTINUE
C
DO 10 I=1,NFPAR
PAR(ICP(I))=PAR(ICP(I))+EP
CALL FBBL(NDIM,PAR,ICP,NBC,U0,U1,F2XX,DFUXX)
DO 9 J=1,NBC
DBC(J,2*NDIM+ICP(I))=(F2XX(J)-F(J))/EP
9 CONTINUE
PAR(ICP(I))=PAR(ICP(I))-EP
10 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FBBL(NDIM,PAR,ICP,NBC,U0,U1,F,DBC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION PAR(1),ICP(1),U0(NDIM),U1(NDIM),F(NBC)
DIMENSION DBC(NBC0,1)
C
IJAC=1
NFPX=NFPAR/2-1
CALL BCNI(NDM,PAR,ICP,NBC0,U0,U1,F,IJAC,DBC)
DO 3 I=1,NBC0
F(NBC0+I)=ZERO
DO 1 J=1,NDM
F(NBC0+I)=F(NBC0+I)+DBC(I,J)*U0(NDM+J)
F(NBC0+I)=F(NBC0+I)+DBC(I,NDM+J)*U1(NDM+J)
1 CONTINUE
IF(NFPX.NE.0) THEN
DO 2 J=1,NFPX
F(NBC0+I)=F(NBC0+I)
* + DBC(I,NDIM+ICP(1+J))*PAR(ICP(NFPAR-NFPX+J))
2 CONTINUE
ENDIF
3 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE ICBL(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,DINT)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates integral conditions for the 2-parameter continuation of
C limit points (BVP).
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(1),UOLD(1),UDOT(1),UPOLD(1),F(NINT),DINT(NINT,1)
DIMENSION ICP(1),PAR(1)
C
C Generate the function.
C
CALL FIBL(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,DFUXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FIBL(NDIM,PAR,ICP,NINT,U1XX,UOLD,UDOT,UPOLD,F1XX,DFUXX)
CALL FIBL(NDIM,PAR,ICP,NINT,U2XX,UOLD,UDOT,UPOLD,F2XX,DFUXX)
DO 3 J=1,NINT
DINT(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 6 I=1,NFPAR
PAR(ICP(I))=PAR(ICP(I))+EP
CALL FIBL(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F1XX,DFUXX)
DO 5 J=1,NINT
DINT(J,NDIM+ICP(I))=(F1XX(J)-F(J))/EP
5 CONTINUE
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FIBL(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,DINT)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(1),UOLD(1),UDOT(1),UPOLD(1),F(NINT),DINT(NINT0,1)
DIMENSION ICP(1),PAR(1)
C
IF(NINT0.GT.0) THEN
NFPX=NFPAR/2-1
IJAC=1
CALL ICNI(NDM,PAR,ICP,NINT0,U,UOLD,UDOT,UPOLD,F,IJAC,DINT)
DO 3 I=1,NINT0
F(NINT0+I)=ZERO
DO 1 J=1,NDM
F(NINT0+I)=F(NINT0+I)+DINT(I,J)*U(NDM+J)
1 CONTINUE
IF(NFPX.NE.0) THEN
DO 2 J=1,NFPX
F(NINT0+I)=F(NINT0+I)
* + DINT(I,NDM+ICP(1+J))*PAR(ICP(NFPAR-NFPX+J))
2 CONTINUE
ENDIF
3 CONTINUE
ENDIF
C
F(NINT)=-ONE
DO 4 I=1,NDM
F(NINT)=F(NINT)+U(NDM+I)*U(NDM+I)
4 CONTINUE
IF(NFPX.NE.0) THEN
DO 5 I=1,NFPX
F(NINT)=F(NINT)+PAR(ICP(NFPAR-NFPX+I))**2
5 CONTINUE
ENDIF
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNBL(NTSTRS,NCOLRS,LAB,IBR,M1U,U,UPS,UDOTPS,UPOLDP,
* TM,DTM,NDIM2,SMAT,RNLLV,IR,IC,F,DFDU,DFDP,NODIR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates starting data for the 2-parameter continuation of limit
C points (BVP).
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),TM(1),DTM(1)
DIMENSION IR(NDIM2),IC(NDIM2),SMAT(NDIM2,NDIM2),RNLLV(NDIM2)
DIMENSION U(NDIM),F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,NPAR)
DIMENSION TEMP(7),ICPRS(20)
C
LOGICAL FOUND
C
CALL FINDL3(IRS,ITP1,NFPAR1,FOUND)
READ(3,*)IBR,NTOT1,ITP1,LAB1,NFPAR1,ISW1,NTPL1,NAR1,NSKIP1,
* NTSTRS,NCOLRS,ICPRS(1)
NRSP1=NTSTRS+1
C
DO 2 J=1,NTSTRS
DO 1 I=1,NCOLRS
K1=(I-1)*NDIM+1
K2=K1+NDM-1
READ(3,101)TEMP(I),(UPS(J,K),K=K1,K2)
1 CONTINUE
TM(J)=TEMP(1)
2 CONTINUE
READ(3,101)TM(NRSP1),(UPS(NRSP1,K),K=1,NDM)
C
NFPAR0=NFPAR/2
READ(3,101)(RLDOT(I),I=1,NFPAR0)
C
C Read U-dot (Derivative with respect to arclength).
C
DO 4 J=1,NTSTRS
DO 3 I=1,NCOLRS
K1=(I-1)*NDIM+NDM+1
K2=I*NDIM
READ(3,101)(UPS(J,K),K=K1,K2)
3 CONTINUE
4 CONTINUE
K1=NDM+1
READ(3,101)(UPS(NRSP1,K),K=K1,NDIM)
C
C Read the parameter values.
C
READ(3,101)(PAR(I),I=1,NPAR)
C
NFPX=NFPAR/2-1
IF(NFPX.GT.0) THEN
DO 5 I=1,NFPX
PAR(ICP(NFPAR0+1+I))=RLDOT(I+1)
5 CONTINUE
ENDIF
C
DO 6 I=1,NFPAR
RL(I)=PAR(ICP(I))
6 CONTINUE
C
NODIR=1
C
101 FORMAT(4X,1P7E18.10)
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Subroutines for the Continuation of Bifurcations to Tori
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNTR(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-parameter continuation of
C bifurcations to tori.
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),ICP(1),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FFTR(NDIM,U,UOLD,ICP,PAR,F,NDM,DFUXX,DFPXX)
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1XX(J)=U(J)
U2XX(J)=U(J)
2 CONTINUE
U1XX(I)=U1XX(I)-EP
U2XX(I)=U2XX(I)+EP
CALL FFTR(NDIM,U1XX,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
CALL FFTR(NDIM,U2XX,UOLD,ICP,PAR,F2XX,NDM,DFUXX,DFPXX)
DO 3 J=1,NDIM
DFDU(J,I)=(F2XX(J)-F1XX(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 6 I=1,NFPAR
PAR(ICP(I))=PAR(ICP(I))+EP
C
CALL FFTR(NDIM,U,UOLD,ICP,PAR,F1XX,NDM,DFUXX,DFPXX)
C
DO 5 J=1,NDIM
DFDP(J,ICP(I))=(F1XX(J)-F(J))/EP
5 CONTINUE
C
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFTR(NDIM,U,UOLD,ICP,PAR,F,NDM,DFDU,DFDP)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDM,NDM),DFDP(NDM,1)
C
IJAC=1
C1=PAR(11)
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
NDM2=2*NDM
DO 2 I=1,NDM
F(NDM+I)=ZERO
F(NDM2+I)=ZERO
DO 1 J=1,NDM
F(NDM+I)=F(NDM+I)+DFDU(I,J)*U(NDM+J)
F(NDM2+I)=F(NDM2+I)+DFDU(I,J)*U(NDM2+J)
1 CONTINUE
F(NDM+I)=C1*F(NDM+I)
F(NDM2+I)=C1*F(NDM2+I)
F(I)=C1*F(I)
2 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE BCTR(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION PAR(20),ICP(20),U0(NBC),U1(NBC),F(NBC)
DIMENSION DBC(NBC,1)
C
NDM2=2*NDM
SIG=PAR(12)
C
SS=DSIN(SIG)
CSGLE SS= SIN(SIG)
CS=DCOS(SIG)
CSGLE CS= COS(SIG)
C
DO 1 I=1,NDM
F(I)=U0(I)-U1(I)
F(NDM+I)=U1(NDM+I)-CS*U0(NDM+I)+SS*U0(NDM2+I)
F(NDM2+I)=U1(NDM2+I)-CS*U0(NDM2+I)-SS*U0(NDM+I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=2*NDIM+NPAR
DO 3 I=1,NBC
DO 2 J=1,NN
DBC(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDM
DBC(I,I)=1
DBC(I,NDIM+I)=-1
DBC(NDM+I,NDM+I)=-CS
DBC(NDM+I,NDM2+I)=SS
DBC(NDM+I,NDIM+NDM+I)=1
DBC(NDM+I,2*NDIM+ICP(4))=CS*U0(NDM2+I)+SS*U0(NDM+I)
DBC(NDM2+I,NDM+I)=-SS
DBC(NDM2+I,NDM2+I)=-CS
DBC(NDM2+I,NDIM+NDM2+I)=1
DBC(NDM2+I,2*NDIM+ICP(4))=SS*U0(NDM2+I)-CS*U0(NDM+I)
4 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE ICTR(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,
* DINT)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(1),UOLD(1),UDOT(1),UPOLD(1),F(3),DINT(NINT,1)
DIMENSION ICP(1),PAR(20)
C
F(1)=ZERO
F(2)=ZERO
F(3)=-PAR(13)
C
NDM2=2*NDM
DO 1 I=1,NDM
F(1)=F(1)+U(I)*UPOLD(I)
F(2)=F(2)+U(NDM+I)*UOLD(NDM2+I)-U(NDM2+I)*UOLD(NDM+I)
F(3)=F(3)+U(NDM+I)*U(NDM+I)+U(NDM2+I)*U(NDM2+I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=NDIM+NPAR
DO 3 I=1,NINT
DO 2 J=1,NN
DINT(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDM
DINT(1,I)=UPOLD(I)
DINT(2,NDM+I)=UOLD(NDM2+I)
DINT(2,NDM2+I)=-UOLD(NDM+I)
DINT(3,NDM+I)=TWO*U(NDM+I)
DINT(3,NDM2+I)=TWO*U(NDM2+I)
4 CONTINUE
C
DINT(3,NDIM+13)=-ONE
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNTR(NTSTRS,NCOLRS,LAB,IBR,M1U,U,UPS,UDOTPS,UPOLDP,
* TM,DTM,NDIM2,SMAT,RNLLV,IR,IC,F,DFDU,DFDP,NODIR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates starting data for the 2-parameter continuation of torus
C bifurcations.
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),TM(1),DTM(1)
DIMENSION IR(NDIM2),IC(NDIM2),SMAT(NDIM2,NDIM2),RNLLV(NDIM2)
DIMENSION U(NDIM),F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,NPAR)
DIMENSION TEMP(7),ICPRS(20)
C
LOGICAL FOUND,EOF3
C
CALL FINDL3(IRS,ITP1,NFPAR1,FOUND)
READ(3,*)IBR,NTOT1,ITP1,LAB1,NFPAR1,ISW1,NTPL1,NAR1,NSKIP1,
* NTSTRS,NCOLRS,ICPRS(1),ICPRS(2)
NRSP1=NTSTRS+1
C
DO 3 J=1,NTSTRS
DO 2 I=1,NCOLRS
K1=(I-1)*NDIM+1
K2=K1+NDM-1
READ(3,101)TEMP(I),(UPS(J,K),K=K1,K2)
C
K2P1=K2+1
K3=K2+NDM
DO 1 K=K2P1,K3
UPS(J,K)=RSMALL*DSIN(TEMP(I))
CSGLE UPS(J,K)=RSMALL* SIN(TEMP(I))
UPS(J,K+NDM)=ZERO
1 CONTINUE
C
2 CONTINUE
TM(J)=TEMP(1)
3 CONTINUE
READ(3,101)TM(NRSP1),(UPS(NRSP1,K),K=1,NDM)
C
DO 4 I=1,NDM
UPS(NRSP1,NDM+I)=ZERO
UPS(NRSP1,2*NDM+I)=ZERO
4 CONTINUE
C
C
READ(3,101)RLDOT(1),RLDOT(3)
RLDOT(2)=ZERO
RLDOT(4)=ZERO
C
C Read the direction vector and initialize the starting direction.
C
DO 7 J=1,NTSTRS
DO 6 I=1,NCOLRS
K1=(I-1)*NDIM+1
K2=K1+NDM-1
READ(3,101)(UDOTPS(J,K),K=K1,K2)
C
K2P1=K2+1
K3=K2+NDM
DO 5 K=K2P1,K3
UDOTPS(J,K)=ZERO
UDOTPS(J,K+NDM)=ZERO
5 CONTINUE
C
6 CONTINUE
7 CONTINUE
READ(3,101)(UDOTPS(NRSP1,K),K=1,NDM)
C
DO 8 I=1,NDM
UDOTPS(NRSP1,NDM+I)=ZERO
UDOTPS(NRSP1,2*NDM+I)=ZERO
8 CONTINUE
C
C Read the parameter values.
C
READ(3,101)(PAR(I),I=1,NPAR)
PAR(13)=ZERO
DO 9 I=1,NFPAR
RL(I)=PAR(ICP(I))
9 CONTINUE
C
NODIR=0
C
101 FORMAT(4X,1P7E18.10)
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Periodic Solutions and Fixed Period Orbits
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNPS(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the continuation of periodic orbits.
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
CALL FUNI(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
RHO=PAR(ICP(2))
DO 1 I=1,NDIM
DFDP(I,ICP(2))=F(I)
F(I)=RHO*DFDP(I,ICP(2))
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
DO 3 I=1,NDIM
DO 2 J=1,NDIM
DFDU(I,J)=RHO*DFDU(I,J)
2 CONTINUE
DFDP(I,ICP(1))=RHO*DFDP(I,ICP(1))
3 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FNFP(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the 2-parameter continuation of
C periodic orbits of fixed period.
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
C
C Generate the function.
C
RHO=PAR(ICP(3))
CALL FUNI(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 1 I=1,NDIM
F(I)=RHO*F(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
C Generate the Jacobian.
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
DFDU(I,J)=RHO*DFDU(I,J)
2 CONTINUE
DO 3 J=1,2
DFDP(I,ICP(J))=RHO*DFDP(I,ICP(J))
3 CONTINUE
4 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE BCPS(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
C
DIMENSION PAR(1),ICP(1),U0(NBC),U1(NBC),F(NBC)
DIMENSION DBC(NBC,1)
C
DO 1 I=1,NDIM
F(I)=U0(I)-U1(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=2*NDIM+NPAR
DO 3 I=1,NBC
DO 2 J=1,NN
DBC(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDIM
DBC(I,I)=ONE
DBC(I,NDIM+I)=-ONE
4 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE ICPS(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,
* DINT)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(1),UOLD(1),UDOT(1),UPOLD(1),F(1),DINT(NINT,1)
DIMENSION ICP(1),PAR(1)
C
F(1)=ZERO
DO 1 I=1,NDIM
F(1)=F(1)+U(I)*UPOLD(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
NN=NDIM+NPAR
DO 3 I=1,NINT
DO 2 J=1,NN
DINT(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDIM
DINT(1,I)=UPOLD(I)
4 CONTINUE
C
RETURN
END
C
C ---------- -----
SUBROUTINE PDBLE(NDIM,NTST,NCOL,M1U,UPS,UDOTPS,TM,RHO)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Preprocesses restart data for switching branches at a period doubling
C bifurcation.
C
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
DIMENSION TM(1),UPS(M1U,1),UDOTPS(M1U,1)
C
RHO=TWO*RHO
C
DO 1 I=1,NTST
TM(I)=HALF*TM(I)
TM(NTST+I)=HALF+TM(I)
1 CONTINUE
C
TM(2*NTST+1)=ONE
C
DO 3 J=1,NTST+1
DO 2 I1=1,NDIM
DO 2 I2=1,NCOL
I=(I2-1)*NDIM+I1
UPS(NTST+J,I)= UPS(NTST+1,I1)+ UPS(J,I)- UPS(1,I1)
UDOTPS(NTST+J,I)=UDOTPS(NTST+1,I1)+UDOTPS(J,I)-UDOTPS(1,I1)
2 CONTINUE
3 CONTINUE
C
NTST=2*NTST
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNPS(NTSTRS,NCOLRS,LAB,IBR,M1U,U,UPS,UDOTPS,UPOLDP,
* TM,DTM,NDIM2,SMAT,RNLLV,IR,IC,F,DFDU,DFDP,NODIR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates starting data for the continuation of a branch of periodic
C solutions from a Hopf bifurcation point.
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),TM(1),DTM(1)
DIMENSION IR(NDIM2),IC(NDIM2),SMAT(NDIM2,NDIM2),RNLLV(NDIM2)
DIMENSION U(NDIM),F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,NPAR)
DIMENSION ICPRS(20)
C
LOGICAL FOUND
C
CALL FINDL3(IRS,ITP,NFPAR1,FOUND)
CALL READL3(IPS,IBR,U,PAR)
C
DO 1 I=1,NFPAR
RL(I)=PAR(ICP(I))
1 CONTINUE
C
RHO=PAR(11)
TPI=PI(TWO)
RIMHB=TPI/RHO
NTSTRS=NTST
NCOLRS=NCOL
C
DO 3 I=1,NDIM2
DO 2 J=1,NDIM2
SMAT(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDIM
SMAT(I,I)=-RIMHB
SMAT(NDIM+I,NDIM+I)=RIMHB
4 CONTINUE
C
IJAC=1
CALL FUNI(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 6 I=1,NDIM
DO 5 J=1,NDIM
SMAT(I,NDIM+J)=DFDU(I,J)
SMAT(NDIM+I,J)=DFDU(I,J)
5 CONTINUE
6 CONTINUE
C
CALL NLVC(NDIM2,NDIM2,2,SMAT,RNLLV,IR,IC)
CALL NRMLZ(NDIM2,RNLLV)
C
C Generate the (initially uniform) mesh.
C
CALL MSH(TM)
DT=ONE/NTST
C
DO 8 J=1,NTST+1
T=TM(J)
S=DSIN(TPI*T)
CSGLE S=SIN(TPI*T)
C=DCOS(TPI*T)
CSGLE C=COS(TPI*T)
DO 7 K=1,NDIM
UDOTPS(J,K)=S*RNLLV(K)+C*RNLLV(NDIM+K)
UPOLDP(J,K)=C*RNLLV(K)-S*RNLLV(NDIM+K)
UPS(J,K)=U(K)
7 CONTINUE
8 CONTINUE
C
DO 11 I=1,NCOL-1
DO 10 J=1,NTST
T=TM(J)+I*( TM(J+1)-TM(J) )/NCOL
S=DSIN(TPI*T)
CSGLE S=SIN(TPI*T)
C=DCOS(TPI*T)
CSGLE C=COS(TPI*T)
DO 9 K=1,NDIM
K1=I*NDIM+K
UDOTPS(J,K1)=S*RNLLV(K)+C*RNLLV(NDIM+K)
UPOLDP(J,K1)=C*RNLLV(K)-S*RNLLV(NDIM+K)
UPS(J,K1)=U(K)
9 CONTINUE
10 CONTINUE
11 CONTINUE
C
RLDOT(1)=ZERO
RLDOT(2)=ZERO
C
DO 12 I=1,NTST
DTM(I)=DT
12 CONTINUE
C
CALL SCALEB(M1U,UDOTPS,RLDOT,DTM)
C
NODIR=-1
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Travelling Wave Solutions to Diffusive Systems
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNWS(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Sets up equations for the continuation of spatially homogeneous
C solutions to parabolic systems, for the purpose of finding
C bifurcations to travelling wave solutions.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1),ICP(1)
C
C Generate the function.
C
NDM2=NDM/2
CALL FFWS(NDIM,U,UOLD,ICP,NFPAR,PAR,IJAC,F,DFDU,DFDP,NDM2,
* DFUXX,DFPXX)
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFWS(NDIM,U,UOLD,ICP,NFPAR,PAR,IJAC,F,DFDU,DFDP,
* NDM,DFUXX,DFPXX)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
DIMENSION DFUXX(NDM,NDM),DFPXX(NDM,20)
C
C=PAR(10)
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F,DFUXX,DFPXX)
C
DO 1 I=1,NDM
F(NDM+I)=-( C*U(NDM+I) + F(I) )/PAR(14+I)
F(I)=U(NDM+I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
DO 3 I=1,NDM
DO 2 J=1,NDM
DFDU(I,J) =ZERO
DFDU(I,J+NDM) =ZERO
DFDU(I+NDM,J) =-DFUXX(I,J)/PAR(14+I)
DFDU(I+NDM,J+NDM)=ZERO
2 CONTINUE
DFDU(I,I+NDM) =ONE
DFDU(I+NDM,I+NDM)=-C/PAR(14+I)
IF(ICP(1).LT.10)THEN
DFDP(I,ICP(1)) =ZERO
DFDP(I+NDM,ICP(1))=-DFPXX(I,ICP(1))/PAR(14+I)
ENDIF
IF(NFPAR.GT.1.AND.ICP(2).LT.10)THEN
DFDP(I,ICP(2)) =ZERO
DFDP(I+NDM,ICP(2))=-DFPXX(I,ICP(2))/PAR(14+I)
ENDIF
3 CONTINUE
C
C Derivative with respect to the wave speed.
C
DO 4 I=1,NDM
DFDP(I,10) =ZERO
DFDP(I+NDM,10)=-U(NDM+I)/PAR(14+I)
4 CONTINUE
C
C Derivatives with respect to the diffusion coefficients.
C
DO 6 J=1,NDM
DO 5 I=1,NDM
DFDP(I,14+J) =ZERO
DFDP(I+NDM,14+J)=ZERO
5 CONTINUE
DFDP(J+NDM,14+J)=-F(J+NDM)/PAR(14+J)
6 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FNWP(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Equations for the continuation of traveling waves.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1),ICP(20)
C
C Generate the function and Jacobian.
C
CALL FNWS(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
RHO=PAR(ICP(2))
DO 1 I=1,NDIM
DFDP(I,ICP(2))=F(I)
F(I)=RHO*F(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
DO 2 I=1,NDIM
DO 2 J=1,NDIM
DFDU(I,J)=RHO*DFDU(I,J)
2 CONTINUE
C
DO 3 I=1,NDIM
DFDP(I,ICP(1))=RHO*DFDP(I,ICP(1))
3 CONTINUE
C
RETURN
END
C
C ---------- ----
SUBROUTINE FNWF(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for the two parameter continuation
C of wavetrains of fixed wave length.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION U(NDIM),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1),ICP(20)
C
C Generate the function and Jacobian.
C
CALL FNWS(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
RHO=PAR(ICP(3))
DO 1 I=1,NDIM
F(I)=RHO*F(I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
DO 2 I=1,NDIM
DO 2 J=1,NDIM
DFDU(I,J)=RHO*DFDU(I,J)
2 CONTINUE
C
DO 3 I=1,NDIM
DO 3 J=1,2
DFDP(I,ICP(J))=RHO*DFDP(I,ICP(J))
3 CONTINUE
C
RETURN
END
C
C ---------- ------
SUBROUTINE STPNWP(NTSTRS,NCOLRS,LAB,IBR,M1U,U,UPS,UDOTPS,UPOLDP,
* TM,DTM,NDIM2,SMAT,RNLLV,IR,IC,F,DFDU,DFDP,NODIR)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates starting data for the continuation of a branch of periodic
C solutions starting from a Hopf bifurcation point (Waves).
C
COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20)
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLLIM/ NMX,NUZR,RL0,RL1,A0,A1
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
C
DIMENSION UPS(M1U,1),UDOTPS(M1U,1),UPOLDP(M1U,1),TM(1),DTM(1)
DIMENSION IR(NDIM2),IC(NDIM2),SMAT(NDIM2,NDIM2),RNLLV(NDIM2)
DIMENSION U(NDIM),F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,NPAR)
DIMENSION ICPRS(20)
C
LOGICAL FOUND
C
CALL FINDL3(IRS,ITP,NFPAR1,FOUND)
CALL READL3(IPS,IBR,U,PAR)
C
DO 1 I=1,NFPAR
RL(I)=PAR(ICP(I))
1 CONTINUE
C
RHO=PAR(11)
TPI=PI(TWO)
RIMHB=TPI/RHO
NTSTRS=NTST
NCOLRS=NCOL
C
DO 3 I=1,NDIM2
DO 2 J=1,NDIM2
SMAT(I,J)=ZERO
2 CONTINUE
3 CONTINUE
C
DO 4 I=1,NDIM
SMAT(I,I)=-RIMHB
SMAT(NDIM+I,NDIM+I)=RIMHB
4 CONTINUE
C
IJAC=1
CALL FNWS(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
DO 6 I=1,NDIM
DO 5 J=1,NDIM
SMAT(I,NDIM+J)=DFDU(I,J)
SMAT(NDIM+I,J)=DFDU(I,J)
5 CONTINUE
6 CONTINUE
C
CALL NLVC(NDIM2,NDIM2,2,SMAT,RNLLV,IR,IC)
CALL NRMLZ(NDIM2,RNLLV)
C
C Generate the (initially uniform) mesh.
C
CALL MSH(TM)
DT=ONE/NTST
C
DO 8 J=1,NTST+1
T=TM(J)
S=DSIN(TPI*T)
CSGLE S=SIN(TPI*T)
C=DCOS(TPI*T)
CSGLE C=COS(TPI*T)
DO 7 K=1,NDIM
UDOTPS(J,K)=S*RNLLV(K)+C*RNLLV(NDIM+K)
UPOLDP(J,K)=C*RNLLV(K)-S*RNLLV(NDIM+K)
UPS(J,K)=U(K)
7 CONTINUE
8 CONTINUE
C
DO 11 I=1,NCOL-1
DO 10 J=1,NTST
T=TM(J)+I*( TM(J+1)-TM(J) )/NCOL
S=DSIN(TPI*T)
CSGLE S=SIN(TPI*T)
C=DCOS(TPI*T)
CSGLE C=COS(TPI*T)
DO 9 K=1,NDIM
K1=I*NDIM+K
UDOTPS(J,K1)=S*RNLLV(K)+C*RNLLV(NDIM+K)
UPOLDP(J,K1)=C*RNLLV(K)-S*RNLLV(NDIM+K)
UPS(J,K1)=U(K)
9 CONTINUE
10 CONTINUE
11 CONTINUE
C
RLDOT(1)=ZERO
RLDOT(2)=ZERO
C
DO 12 I=1,NTST
DTM(I)=DT
12 CONTINUE
C
CALL SCALEB(M1U,UDOTPS,RLDOT,DTM)
C
NODIR=-1
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Diffusive Systems : Time Evolution
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FNPE(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Generates the equations for taking one time step (Implicit Euler).
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLWIF/ U0XX(152),U1XX(152),U2XX(152),F1XX(152),F2XX(152),
* DFUXX(50,120),DFPXX(50,20),DDUXX(50),DDPXX(20)
C
DIMENSION U(NDIM),UOLD(NDIM),PAR(1)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1),ICP(1)
C
C Generate the function and Jacobian.
C
CALL FFPE(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP,
* NDM,DFUXX,DFPXX)
C
RETURN
END
C
C ---------- ----
SUBROUTINE FFPE(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP,
* NDM,DFUXX,DFPXX)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLCRL/ RDSOLD,A,RL(20),RLOLD(20),RLDOT(20)
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS
C
DIMENSION U(NDIM),UOLD(NDIM),ICP(20),PAR(20)
DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,1)
DIMENSION DFUXX(NDM,NDM),DFPXX(NDM,20)
C
RHO=PAR(11)
T=PAR(ICP(1))
DT=T-RLOLD(1)
IF(DABS(DT).LT.DSMIN)DT=DS
CSGLE IF( ABS(DT).LT.DSMIN)DT=DS
C
CALL FUNI(NDM,U,UOLD,ICP,PAR,IJAC,F(NDM+1),DFUXX,DFPXX)
C
DO 1 I=1,NDM
F(I)=RHO*U(NDM+I)
F(NDM+I)=RHO*( (U(I)-UOLD(I))/DT - F(NDM+I) )/PAR(14+I)
1 CONTINUE
C
IF(IJAC.EQ.0)RETURN
C
DO 3 I=1,NDM
DO 2 J=1,NDM
DFDU(I,J) =ZERO
DFDU(I,J+NDM) =ZERO
DFDU(I+NDM,J) =-RHO*DFUXX(I,J)/PAR(14+I)
DFDU(I+NDM,J+NDM)=ZERO
2 CONTINUE
DFDU(I,I+NDM) =RHO
DFDU(I+NDM,I) =DFDU(I+NDM,I) + RHO/(DT*PAR(14+I))
DFDP(I,ICP(1)) =ZERO
DFDP(I+NDM,ICP(1))=-RHO*(U(I)-UOLD(I))/(DT**2*PAR(14+I))
3 CONTINUE
C
RETURN
END
C
SUBROUTINE ICPE(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,DINT)
C ---------- ----
C
C Dummy integral condition subroutine for parabolic systems.
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Routines for Interface with User Supplied Routines
C (To generate Jacobian by differencing, if not supplied analytically)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C ---------- ----
SUBROUTINE FUNI(NDIM,U,UOLD,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Interface subroutine to user supplied FUNC.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDIF/ U1ZZ(50),U2ZZ(50),F1ZZ(50),F2ZZ(50)
C
DIMENSION U(NDIM),UOLD(NDIM),ICP(20),PAR(20),F(NDIM)
DIMENSION DFDU(NDIM,NDIM),DFDP(NDIM,20)
C
C Generate the function.
C
CALL FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IF(JAC.EQ.1 .OR. IJAC.EQ.0)RETURN
C
C Generate the Jacobian by differencing.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1ZZ(J)=U(J)
U2ZZ(J)=U(J)
2 CONTINUE
U1ZZ(I)=U1ZZ(I)-EP
U2ZZ(I)=U2ZZ(I)+EP
CALL FUNC(NDIM,U1ZZ,ICP,PAR,0,F1ZZ,DFDU,DFDP)
CALL FUNC(NDIM,U2ZZ,ICP,PAR,0,F2ZZ,DFDU,DFDP)
DO 3 J=1,NDIM
DFDU(J,I)=(F2ZZ(J)-F1ZZ(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 6 I=1,NFPAR
EP=HMACH*( ONE +DABS(PAR(ICP(I))) )
CSGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) )
PAR(ICP(I))=PAR(ICP(I))+EP
CALL FUNC(NDIM,U,ICP,PAR,0,F1ZZ,DFDU,DFDP)
DO 5 J=1,NDIM
DFDP(J,ICP(I))=(F1ZZ(J)-F(J))/EP
5 CONTINUE
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
SUBROUTINE BCNI(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C ---------- ----
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Interface subroutine to the user supplied BCND.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDIF/ U1ZZ(50),U2ZZ(50),F1ZZ(50),F2ZZ(50)
C
DIMENSION PAR(20),ICP(20),U0(NDIM),U1(NDIM),F(NBC),DBC(NBC,1)
C
C Generate the function.
C
CALL BCND(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
C
IF(JAC.EQ.1 .OR. IJAC.EQ.0)RETURN
C
C Generate the Jacobian by differencing.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U0(I)).GT.UMX)UMX=DABS(U0(I))
CSGLE IF( ABS(U0(I)).GT.UMX)UMX= ABS(U0(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1ZZ(J)=U0(J)
U2ZZ(J)=U0(J)
2 CONTINUE
U1ZZ(I)=U1ZZ(I)-EP
U2ZZ(I)=U2ZZ(I)+EP
CALL BCND(NDIM,PAR,ICP,NBC,U1ZZ,U1,F1ZZ,0,DBC)
CALL BCND(NDIM,PAR,ICP,NBC,U2ZZ,U1,F2ZZ,0,DBC)
DO 3 J=1,NBC
DBC(J,I)=(F2ZZ(J)-F1ZZ(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
UMX=ZERO
DO 5 I=1,NDIM
IF(DABS(U1(I)).GT.UMX)UMX=DABS(U1(I))
CSGLE IF( ABS(U1(I)).GT.UMX)UMX= ABS(U1(I))
5 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 8 I=1,NDIM
DO 6 J=1,NDIM
U1ZZ(J)=U1(J)
U2ZZ(J)=U1(J)
6 CONTINUE
U1ZZ(I)=U1ZZ(I)-EP
U2ZZ(I)=U2ZZ(I)+EP
CALL BCND(NDIM,PAR,ICP,NBC,U0,U1ZZ,F1ZZ,0,DBC)
CALL BCND(NDIM,PAR,ICP,NBC,U0,U2ZZ,F2ZZ,0,DBC)
DO 7 J=1,NBC
DBC(J,NDIM+I)=(F2ZZ(J)-F1ZZ(J))/(2*EP)
7 CONTINUE
8 CONTINUE
C
DO 10 I=1,NFPAR
EP=HMACH*( ONE +DABS(PAR(ICP(I))) )
CSGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) )
PAR(ICP(I))=PAR(ICP(I))+EP
CALL BCND(NDIM,PAR,ICP,NBC,U0,U1,F1ZZ,0,DBC)
DO 9 J=1,NBC
DBC(J,2*NDIM+ICP(I))=(F1ZZ(J)-F(J))/EP
9 CONTINUE
PAR(ICP(I))=PAR(ICP(I))-EP
10 CONTINUE
C
RETURN
END
C
SUBROUTINE ICNI(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,DINT)
C ---------- ----
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDIF/ U1ZZ(50),U2ZZ(50),F1ZZ(50),F2ZZ(50)
C
C Interface subroutine to user supplied ICND.
C
DIMENSION U(NDIM),UOLD(NDIM),UDOT(NDIM),UPOLD(NDIM)
DIMENSION F(NINT),DINT(NINT,1),ICP(20),PAR(20)
C
C Generate the integrand.
C
CALL ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F,IJAC,DINT)
C
IF(JAC.EQ.1 .OR. IJAC.EQ.0)RETURN
C
C Generate the Jacobian by differencing.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1ZZ(J)=U(J)
U2ZZ(J)=U(J)
2 CONTINUE
U1ZZ(I)=U1ZZ(I)-EP
U2ZZ(I)=U2ZZ(I)+EP
CALL ICND(NDIM,PAR,ICP,NINT,U1ZZ,UOLD,UDOT,UPOLD,F1ZZ,0,DINT)
CALL ICND(NDIM,PAR,ICP,NINT,U2ZZ,UOLD,UDOT,UPOLD,F2ZZ,0,DINT)
DO 3 J=1,NINT
DINT(J,I)=(F2ZZ(J)-F1ZZ(J))/(2*EP)
3 CONTINUE
4 CONTINUE
C
DO 6 I=1,NFPAR
EP=HMACH*( ONE +DABS(PAR(ICP(I))) )
CSGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) )
PAR(ICP(I))=PAR(ICP(I))+EP
CALL ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,F1ZZ,0,DINT)
DO 5 J=1,NINT
DINT(J,NDIM+ICP(I))=(F1ZZ(J)-F(J))/EP
5 CONTINUE
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
SUBROUTINE FOPI(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
C ---------- ----
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Interface subroutine to user supplied FOPT.
C
COMMON /BLICN/ NDM,NDMP1,NROW,NCLM,NRC,NCC,NPAR,NFPAR,NBC0,NINT0
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC
COMMON /BLRCN/ HALF,ZERO,ONE,TWO,HMACH,RSMALL,RLARGE
COMMON /BLDIF/ U1ZZ(50),U2ZZ(50),F1ZZ(50),F2ZZ(50)
C
DIMENSION U(NDIM),ICP(20),PAR(20),DFDU(NDIM),DFDP(20)
C
C Generate the objective function.
C
CALL FOPT(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
C
IF(JAC.EQ.1 .OR. IJAC.EQ.0)RETURN
C
C Generate the Jacobian by differencing.
C
UMX=ZERO
DO 1 I=1,NDIM
IF(DABS(U(I)).GT.UMX)UMX=DABS(U(I))
CSGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I))
1 CONTINUE
C
EP=HMACH*(ONE+UMX)
C
DO 4 I=1,NDIM
DO 2 J=1,NDIM
U1ZZ(J)=U(J)
U2ZZ(J)=U(J)
2 CONTINUE
U1ZZ(I)=U1ZZ(I)-EP
U2ZZ(I)=U2ZZ(I)+EP
CALL FOPT(NDIM,U1ZZ,ICP,PAR,0,F1,DFDU,DFDP)
CALL FOPT(NDIM,U2ZZ,ICP,PAR,0,F2,DFDU,DFDP)
DFDU(I)=(F2-F1)/(2*EP)
4 CONTINUE
C
DO 6 I=1,NFPAR
EP=HMACH*( ONE +DABS(PAR(ICP(I))) )
CSGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) )
PAR(ICP(I))=PAR(ICP(I))+EP
CALL FOPT(NDIM,U,ICP,PAR,0,F1,DFDU,DFDP)
DFDP(ICP(I))=(F1-F)/EP
PAR(ICP(I))=PAR(ICP(I))-EP
6 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Installation Dependent Subroutines for Timing AUTO
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
SUBROUTINE AUTIM0(T)
C ---------- ------
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Set initial time for measuring CPU time used.
C
T=0
C
RETURN
END
C
SUBROUTINE AUTIM1(T)
C ---------- ------
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CSGLE IMPLICIT REAL (A-H,O-Z)
C
C Set final time for measuring CPU time used.
C
T=0
C
RETURN
END
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------