SUBROUTINE DBCMAV(N,MA,MB,BCV,IHOM,NRTI,Q,NTI,PHIREC,NU,Z,KKK, 1 R,BB,WS) C -------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION MA(N,N),MB(N,N) DIMENSION BCV(N),Q(N,N,NTI),PHIREC(NU,NTI),Z(N,NTI), 1 R(N,N),BB(N),WS(N,N) INTEGER KKK(N) LOGICAL HOM0,HOM1 C HOM0 = IHOM.EQ.0 HOM1 = IHOM.NE.0 DO 1000 I = 1 , N KKK(I) = I*(I+1)/2 1000 CONTINUE DO 1010 I = 1 , N BB(I) = BCV(I) 1010 CONTINUE C COMPUTATION OF MA.Q(1) STORED IN WS DO 1120 K = 1 , N CALL DMATVC(MA,N,N,N,N,Q(1,K,1),WS(1,K)) 1120 CONTINUE IF (HOM1) THEN C COMPUTATION OF MA.Q(1).Z(1) AND BCV-MA.Q(1).Z(1) STORED IN BB. DO 1210 I = 1 , N SOM = 0.D0 DO 1220 J = 1 , N SOM = SOM + WS(I,J) * Z(J,1) 1220 CONTINUE BB(I) = BB(I) - SOM 1210 CONTINUE ENDIF C COMPUTATION OF MA.Q(1).PHIREC(1) STORED IN R DO 1420 K = 1 , N NK = KKK(K) - K DO 1410 I = 1 , N SOM = 0.D+0 DO 1400 J = 1 , K NJ = NK + J SOM = SOM + WS(I,J) * PHIREC(NJ,1) 1400 CONTINUE R(I,K) = SOM 1410 CONTINUE 1420 CONTINUE C COMPUTATION OF MB.Q(N) STORED IN WS DO 1520 K = 1 , N CALL DMATVC(MB,N,N,N,N,Q(1,K,NRTI),WS(1,K)) 1520 CONTINUE IF (HOM1) THEN C COMPUTATION OF MB.Q(N).Z(N) AND BCV-MB.Q(N).Z(N) STORED IN BB DO 1610 I = 1 , N SOM = 0.D0 DO 1600 J = 1 , N SOM = SOM + WS(I,J) * Z(J,NRTI) 1600 CONTINUE BB(I) = BB(I) - SOM 1610 CONTINUE ENDIF C COMPUTATION OF MB.Q(N).PHIREC(N) AND C MA.Q(1).PHIREC(1) + MB.Q(N).PHIREC(N) STORED IN R DO 1820 K = 1 , N NK = KKK(K) - K DO 1810 I = 1 , N SOM = 0.D+0 DO 1800 J = 1 , K NJ = NK + J SOM = SOM + WS(I,J) * PHIREC(NJ,NRTI) 1800 CONTINUE R(I,K) = R(I,K) + SOM 1810 CONTINUE 1820 CONTINUE RETURN C END OF DBCMAV END SUBROUTINE DFQUS(PHI,N,IRM,IRN,DIAG,BET,K,KKK,U,NU,NTI,Q, 1 WPHI0,WSUM) C ---------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION PHI(IRM,IRN),DIAG(IRN),BET(IRN),U(NU,NTI), 1 Q(IRM,IRM,NTI),WPHI0(IRM,IRM),WSUM(IRN) INTEGER KKK(N) C CALL DQEVAL(N,N,PHI,IRM,IRN,BET,WPHI0,WSUM) DO 1000 I = 1 , N DO 1000 J = 1 , N Q(I,J,K) = WPHI0(I,J) 1000 CONTINUE U(1,K) = DIAG(1) IF (N.LT.2) GOTO 2000 DO 1200 I = 2 , N NI = KKK(I) IMIN1 = I - 1 NJ = NI - I U(NI,K) = DIAG(I) DO 1100 J = 1 , IMIN1 NK = NJ + J U(NK,K) = PHI(J,I) 1100 CONTINUE 1200 CONTINUE 2000 DO 2300 I = 1 , N NI = KKK(I) IF (U(NI,K).GE.0.D0) GOTO 2300 DO 2100 J = 1 , N Q(J,I,K) = -Q(J,I,K) 2100 CONTINUE DO 2200 J = I , N NJ = KKK(J) - J + I U(NJ,K) = -U(NJ,K) 2200 CONTINUE 2300 CONTINUE RETURN C END OF DFQUS END SUBROUTINE DFUNRC(N,KPART,NRTI,U,NU,NTI,KKK,PHIREC,W1,W2, 1 IERROR) C ---------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION PHIREC(NU,NTI),U(NU,NTI),W1(N),W2(N) INTEGER KKK(N) C IERROR = 0 IF (KPART.LT.0) KPART=0 IF (KPART.GT.N) KPART=N KPPL1 = KPART + 1 C COMPUTATION OF THE PHI1'S FOR I <= KPART IF (KPART.GT.0) THEN DO 1400 K = 1 , KPART KK = KKK(K) NK = KK - K NK1 = NK + 1 DO 1150 I = 1 , K PHIREC(NK + I,NRTI) = 0.0D+0 1150 CONTINUE PHIREC(KK,NRTI) = 1.0D+0 DO 1350 I = 2 , NRTI IN = NRTI + 1 - I IN1 = IN + 1 DO 1250 IJ = 1 , K W2(IJ) = PHIREC(NK + IJ,IN1) 1250 CONTINUE CALL DSOLUP(K,U(1,IN1),W2,W1,IERROR) IF (IERROR.GT.0) GOTO 3000 DO 1300 IJ = 1 , K PHIREC(NK + IJ,IN) = W1(IJ) 1300 CONTINUE 1350 CONTINUE 1400 CONTINUE ENDIF C THE REQUIRED PART OF PHIREC1 IS NOW COMPUTED. C THE REQUIRED PART OF PHIREC2 WILL BE COMPUTED NOW. IF (KPPL1.LE.N) THEN C COMPUTATION OF THE PHI2'S FOR I > KPART DO 1700 K = KPPL1 , N KK = KKK(K) NK = KK - K NK1 = NK + KPPL1 DO 1500 I = NK1 , KK PHIREC(I,1) = 0.0D+0 1500 CONTINUE PHIREC(KK,1) = 1.0D+0 DO 1650 I = 2 , NRTI IN1 = I - 1 DO 1600 IJ = KPPL1 , K SOM = 0.0D+0 DO 1550 IN = IJ , K MK = KKK(IN) - IN SOM = SOM + PHIREC(NK + IN , IN1) * U(MK + IJ , I) 1550 CONTINUE PHIREC(NK + IJ , I) = SOM 1600 CONTINUE 1650 CONTINUE 1700 CONTINUE C COMPUTATION OF THE PHI1'S FOR I > KPART, IF KPART > 0 IF (KPART.GT.0) THEN DO 2050 K = KPPL1 , N KK = KKK(K) NK = KK - K DO 1750 I = 1, KPART PHIREC(NK + I, NRTI) = 0.0D+0 1750 CONTINUE DO 2000 I = 2 , NRTI IN = NRTI + 1 - I IN1 = IN + 1 DO 1850 IJ = 1 , KPART SOM = 0.0D+0 DO 1800 NK1 = KPPL1 , K MK = KKK(NK1 - 1) SOM = SOM + U(MK + IJ,IN1) * PHIREC(NK + NK1 ,IN) 1800 CONTINUE W2(IJ) = - SOM + PHIREC(NK + IJ , IN1) 1850 CONTINUE CALL DSOLUP(KPART,U(1,IN1),W2,W1,IERROR) IF (IERROR.GT.0) GOTO 3000 DO 1950 IJ = 1 , KPART PHIREC(NK + IJ , IN) = W1(IJ) 1950 CONTINUE 2000 CONTINUE 2050 CONTINUE ENDIF ENDIF RETURN 3000 IERROR = 250 RETURN C END OF DFUNRC END SUBROUTINE DGTUR(N,IHOM,NRTI,U,KU,NTI,Q,D,ER,KKK,KPART,IP1,IP2, 1 WDIA,WBET,WSUM,WPHI,WPHI0,IERROR) C --------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(KU,NTI),Q(N,N,NTI),D(N,NTI),ER(5),WDIA(N), 1 WBET(N),WSUM(N),WPHI(N,N),WPHI0(N,N) INTEGER KKK(N),IP1(N),IP2(N) LOGICAL NOSORT C IERROR = 0 ISORT = IP1(1) NKP = N IF (ISORT.NE.0) NKP = IP1(2) IF (NKP.LE.0) THEN KPART = 0 RETURN ENDIF NTI1 = NTI - 1 NOSORT = .FALSE. C.....CHECK WHETHER TRANSFORMATION IS ALLOWED. TOL = DMAX1(ER(1),ER(2)) TOL1 = TOL / ER(3) DO 1100 I = 1 , N IP1(I) = I 1100 CONTINUE DO 1300 I = 1 , NKP IK = KKK(I) WDIA(I) = 0.0D+0 DO 1200 J = 2 , NRTI SOM = DABS(U(IK,J)) IF (SOM.GE.TOL1) NOSORT = .TRUE. WDIA(I) = WDIA(I) + DLOG10(SOM) 1200 CONTINUE 1300 CONTINUE KPART = 0 DO 1400 I = 1 , NKP IF (WDIA(I).GT.0) KPART = KPART + 1 1400 CONTINUE C PARTITIONING IS DETERMINED IF ((KPART.EQ.0).OR.(KPART.EQ.NKP)) RETURN IF (NOSORT) RETURN IFLAG = 0 K = 0 N1 = NKP - 1 C DETERMINATION OF THE PERMUTATION WHICH GLOBALLY ORDERS THE MODES DO 1700 I = 1 , N1 SORT = WDIA(I) IF (SORT.LE.0.D0) THEN IPL1 = I + 1 DO 1500 J = IPL1 , NKP SOM = WDIA(J) IF ((SOM.GT.0.D0).AND.(SOM-SORT.GT.1.D0)) THEN IFLAG = 1 WDIA(J) = SORT WDIA(I) = SOM JM = IP1(I) IP1(I) = IP1(J) IP1(J) = JM GOTO 1600 ENDIF 1500 CONTINUE 1600 CONTINUE ENDIF 1700 CONTINUE IF (IFLAG.EQ.0) RETURN C.....THE PERMUTATION P IS NOW DETERMINED C.....COMPUTATION OF Q1P,U2P,INV(P)D1,INV(P)W1 IERROR = 1 DO 2100 I = 1 , N DO 2100 J = 1 , N Q(I,J,NTI) = Q(I,J,1) WPHI(I,J) = 0.0D+0 2100 CONTINUE DO 2500 I = 1 , N IK = IP1(I) IP1(I) = IP2(IK) IK1 = KKK(IK) - IK DO 2300 J = 1 , N Q(J,I,1) = Q(J,IK,NTI) 2300 CONTINUE DO 2400 J = 1 , IK WPHI(J,I) = U(IK1+J,2) 2400 CONTINUE 2500 CONTINUE DO 2600 I = 1 , N IP2(I) = IP1(I) 2600 CONTINUE C.....TRANSFORMATION OF THE U'S,Q'S,D'S AND WI'S DO 3400 I = 2 , NRTI CALL DQUDEC(WPHI,N,N,N,N,WDIA,WBET,WSUM) CALL DFQUS(WPHI,N,N,N,WDIA,WBET,NTI,KKK,U,KU,NTI,Q, 1 WPHI0,WSUM) DO 2700 J = 1 , KU U(J,I) = U(J,NTI) 2700 CONTINUE IF (IHOM.EQ.0) GOTO 2950 CALL DTAMVC(Q(1,1,NTI),N,N,N,N,D(1,I),WBET) DO 2900 J = 1 , N D(J,I) = WBET(J) 2900 CONTINUE 2950 DO 3000 J = 1 , N CALL DMATVC(Q(1,1,I),N,N,N,N,Q(1,J,NTI),WPHI0(1,J)) 3000 CONTINUE DO 3100 J = 1 , N DO 3100 L = 1 , N Q(J,L,I) = WPHI0(J,L) 3100 CONTINUE I1 = I + 1 IF (I1.GT.NRTI) GOTO 3400 C.....COMPUTATION OF UI+1.OI DO 3300 J = 1 , N DO 3300 L = 1 , N SOM = 0.0D+0 DO 3200 K = L , N IK = KKK(K) - K SOM = SOM + U(IK+L,I1) * Q(K,J,NTI) 3200 CONTINUE WPHI(L,J) = SOM 3300 CONTINUE 3400 CONTINUE RETURN C END OF DGTUR END SUBROUTINE DINPRO(A,NA,B,NB,N,C) C -------------------------------- C DOUBLE PRECISION A(*),B(*),C C C LOCAL VARIABLES INTEGER NN,JL,JJ C = 0.0D+0 IF (N.LE.0) GOTO 2000 NN = 1 - NB JL = NA * N DO 1000 JJ = 1 , JL , NA NN = NN + NB C = C + A(JJ) * B(NN) 1000 CONTINUE 2000 CONTINUE RETURN C END OF DINPRO C C DINPRO IS THE DOUBLE PRECISION VERSION OF SUBROUTINE C VVIPP FROM THE ACCU LIBRARY C END SUBROUTINE DKPCH(N,U,NU,NTI,NRTI,KPART,KKK,ER,IER) C -------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(NU,NTI),ER(5) INTEGER KKK(N) C IER = 0 PSI = DMAX1(ER(1),ER(2)) / ER(3) ER(5) = 0.0D+0 BLOWUQ = 0.0D+0 KPPL1 = KPART + 1 IF (KPART.GE.N) GOTO 1300 DO 1200 K = KPPL1 , N PROD = 1.0D+0 KK = KKK(K) DO 1100 I = 2 , NRTI PROD = 1.0D+0 + PROD * DABS(U(KK,I)) IF (BLOWUQ.LT.PROD) BLOWUQ = PROD 1100 CONTINUE 1200 CONTINUE 1300 IF (KPART.LT.1) GOTO 2000 DO 1500 K = 1 , KPART PROD = 1.0D+0 KK = KKK(K) DO 1400 I = 2 , NRTI IN = NRTI + 2 - I PROD = 1.0D+0 + PROD / DABS(U(KK,IN)) IF (ER(5).LT.PROD) ER(5) = PROD 1400 CONTINUE 1500 CONTINUE IF (KPART.GE.N) GOTO 2000 BLOW1 = 0.0D+0 BLOW2 = 0.0D+0 DO 1700 K = 1 , KPART KK = KKK(K) AMX = 1.0D+0 / DABS(U(KK,NRTI)) AMX2 = AMX IF (NRTI.LT.3) GOTO 1650 DO 1600 I = 3 , NRTI IN = NRTI + 2 - I PROD = 1.0D+0 / DABS(U(KK,IN)) AMX1 = PROD * AMX2 AMX2 = DMAX1(AMX1,PROD) AMX = DMAX1(AMX,AMX2) 1600 CONTINUE 1650 IF (AMX.GT.BLOW1) BLOW1 = AMX 1700 CONTINUE DO 1900 K = KPPL1 , N KK = KKK(K) AMX = DABS(U(KK,2)) AMX2 = AMX DO 1800 I = 3 , NRTI PROD = DABS(U(KK,I)) AMX1 = PROD * AMX2 AMX2 = DMAX1(PROD,AMX1) AMX = DMAX1(AMX,AMX2) 1800 CONTINUE IF (AMX.GT.BLOW2) BLOW2 = AMX 1900 CONTINUE BLOW = BLOW1 * BLOW2 ER(5) = ER(5) + BLOW 2000 ER(5) = DMAX1(ER(5),BLOWUQ) IF (ER(5).LT.PSI) RETURN IER = 240 RETURN C END OF DKPCH END SUBROUTINE DMATVC(A,M,N,IRM,IRN,Y,X) C ------------------------------------ C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(IRM,IRN),Y(IRN),X(IRM) C IF (M.LE.0.OR.N.LE.0) GOTO 2000 DO 1100 I=1,M S = 0.0D+0 DO 1000 J=1,N S = S + A(I,J) * Y(J) 1000 CONTINUE X(I) = S 1100 CONTINUE 2000 CONTINUE RETURN C END OF DMATVC END SUBROUTINE DPSR(N,KPART,NRTI,U,NU,NTI,D,KKK,V,W1,W2,IERROR) C ----------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(NU,NTI),D(N,NTI),V(N,NTI),W1(N),W2(N) INTEGER KKK(N) C IERROR = 0 KPPL1 = KPART + 1 C COMPUTATION OF VI2 IF (KPART.GE.N) GOTO 1500 DO 1300 I = KPPL1 , N V(I,1) = 0.0D+0 1300 CONTINUE DO 1450 I = 2 , NRTI IN = I - 1 DO 1400 K = KPPL1 , N SOM = 0.0D+0 DO 1350 L = K , N MK = KKK(L) - L SOM = SOM + U(MK+K,I) * V(L,IN) MK1 = MK + K 1350 CONTINUE V(K,I) = SOM + D(K,I) 1400 CONTINUE 1450 CONTINUE C COMPUTATION OF VI1 1500 IF (KPART.LT.1) RETURN DO 1600 I = 1 , KPART V(I,NRTI) = 0.0D+0 1600 CONTINUE DO 1850 I = 2 , NRTI IN = NRTI + 1 - I IN1 = IN + 1 DO 1700 K = 1 , KPART SOM = 0.0D+0 IF (KPPL1.GT.N) GOTO 1660 DO 1650 L = KPPL1 , N MK = KKK(L) - L SOM = SOM + U(MK+K,IN1) * V(L,IN) 1650 CONTINUE 1660 W2(K) = - SOM - D(K,IN1) + V(K,IN1) 1700 CONTINUE CALL DSOLUP(KPART,U(1,IN1),W2,W1,IERROR) IF (IERROR.NE.0) GOTO 3000 DO 1800 K = 1 , KPART V(K,IN) = W1(K) 1800 CONTINUE 1850 CONTINUE RETURN 3000 IERROR = 250 RETURN C END OF DPSR END SUBROUTINE DQEVAL(M,N,QR,IRM,IRN,BET,QREV,WD) C -------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION QR(IRM,IRN),QREV(IRM,IRM),BET(IRN),WD(IRM) C DO 1100 I=N,M DO 1050 J=I,M QREV(I,J) = -QR(I,N) * QR(J,N) * BET(N) QREV(J,I) = QREV(I,J) 1050 CONTINUE QREV(I,I) = QREV(I,I) + 1.0D+0 1100 CONTINUE IF (N.LT.2) RETURN DO 1300 IK = 2 , N I = N + 2 - IK C THE NECESARY INNERPRODUCTS OF THE UK'S ARE STORED IN THE C ARRAY WD(1:IRM) C IMIN1 = I - 1 WD(IMIN1) = QR(IMIN1,IMIN1) IM = M - I + 1 DO 1150 J = I , M CALL DINPRO(QR(I,IMIN1),1,QREV(I,J),1,IM,WD(J)) C C COMPUTATION OF (I - BET(I-1).U(I-1)') * PRODUCT C ( ( -BET(K).U(K).U(K)') K= I TO N) C 1150 CONTINUE DO 1200 K = I , M DO 1200 L = I , M QREV(K,L) = QREV(K,L) - BET(IMIN1) * QR(K,IMIN1) * WD(L) 1200 CONTINUE DO 1250 J = I , M QREV(J,IMIN1) = - BET(IMIN1) * QR(J,IMIN1) * WD(IMIN1) QREV(IMIN1,J) = - BET(IMIN1) * QR(IMIN1,IMIN1) * WD(J) 1250 CONTINUE QREV(IMIN1,IMIN1) = 1.0D+0 - BET(IMIN1)*WD(IMIN1)*WD(IMIN1) 1300 CONTINUE RETURN C END OF DQEVAL END SUBROUTINE DQUDEC(A,M,N,IRM,IRN,DIAG,BET,WY) C ---------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(IRM,IRN),DIAG(IRN),BET(IRN),WY(IRN) C IF (M.GT.IRM) GOTO 3000 IF (N.GT.IRN) GOTO 3010 C DO 2000 K = 1 , N K1 = K + 1 IK = M - K + 1 CALL DINPRO(A(K,K),1,A(K,K),1,IK,SIGMA) IF (SIGMA.EQ.0.D0) GOTO 1400 QRKK = A(K,K) ALPHAK = DSQRT(SIGMA) IF (QRKK.GE.0.0D+0) ALPHAK = - ALPHAK DIAG(K) = ALPHAK BETA = 1.0D+0 / (SIGMA - QRKK * ALPHAK) BET(K) = BETA A(K,K) = QRKK - ALPHAK IF (K1.GT.N) GOTO 2000 DO 1200 J = K1 , N CALL DINPRO(A(K,K),1,A(K,J),1,IK,WY(J)) WY(J) = WY(J) * BETA 1200 CONTINUE DO 1300 J = K1 , N DO 1300 I = K , M A(I,J) = A(I,J) - A(I,K) * WY(J) 1300 CONTINUE GOTO 2000 1400 BET(K) = 0.D0 DIAG(K) = 0.D0 2000 CONTINUE RETURN 3000 WRITE(*,100) RETURN 3010 WRITE(*,110) RETURN C 100 FORMAT(' ERROR: M > IRM ; NO DECOMPOSITION DONE') 110 FORMAT(' ERROR: N > IRN ; NO DECOMPOSITION DONE') C C END OF DQUDEC END SUBROUTINE DSBVP(N,IHOM,M0,MN,NRTI,Q,NTI,BCV,PHIREC,NU,V, 1 WI,EPS,COND,KKK,IP,WQ,WS,W1,W2,IER) C ------------------------------------------------------------------ C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION M0(N,N),MN(N,N) DIMENSION Q(N,N,NTI),BCV(N),PHIREC(NU,NTI),V(N,NTI), 1 WI(N,NTI),WQ(N,N),WS(N,N),W1(N),W2(N) INTEGER KKK(N),IP(N) C CALL DBCMAV(N,M0,MN,BCV,IHOM,NRTI,Q,NTI,PHIREC,NU,WI,KKK,WQ, 1 W1,WS) CALL DCROUT(WQ,N,N,W1,EPS,IP,W2,IER) IF (IER.NE.0) THEN IER = 260 RETURN ENDIF C COMPUTATION OF THE CONDITION NUMBER COND = 0.0D+0 DO 1300 I = 1 , N DO 1100 J = 1 , N W2(J) = 0.0D+0 1100 CONTINUE W2(I) = 1.0D+0 CALL DSOLDE(WQ,N,N,IP,W2) SOMIN = 0.0D+0 DO 1200 J = 1 , N SOMIN = SOMIN + DABS(W2(J)) 1200 CONTINUE IF (SOMIN.GT.COND) COND = SOMIN 1300 CONTINUE C COMPUTATION OF EI'=PHIREC(I)Y' AND EI=Q(I)EI' DO 1500 K=1,NRTI DO 1500 I=1,N SOMIN=0.0D+0 DO 1400 J=I,N NJ=KKK(J) - J NI= NJ + I SOMIN = SOMIN + PHIREC(NI,K) * W1(J) 1400 CONTINUE IF (IHOM.EQ.0) THEN WI(I,K) = SOMIN ELSE WI(I,K) = WI(I,K) + SOMIN ENDIF 1500 CONTINUE C COMPUTATION OF Q(I)XI' DO 2000 K=1,NRTI CALL DMATVC(Q(1,1,K),N,N,N,N,WI(1,K),W2) DO 1800 J=1,N V(J,K) = W2(J) 1800 CONTINUE 2000 CONTINUE RETURN C END OF DSBVP END SUBROUTINE DSOLUP(M,U,B,X,IER) C ------------------------------ C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(*),X(*),B(*) C IER = 0 IF (M.LE.0) RETURN M1 = M * (M+1) / 2 IF (U(M1).EQ.0.0D+0) GOTO 1200 X(M) = B(M) / U(M1) IF (M.EQ.1) RETURN MMIN1 = M - 1 DO 1100 I = 1 , MMIN1 MI = M - I SOMIN = 0.0D+0 MJ = (MI + 1) * MI / 2 MIJ = MJ DO 1000 J = 1 , I MIJ = MIJ + MI - 1 + J SOMIN = SOMIN + U(MIJ) * X(MI+J) 1000 CONTINUE IF (U(MJ).EQ.0.0D+0) GOTO 1200 X(MI) = (B(MI) - SOMIN) / U(MJ) 1100 CONTINUE RETURN 1200 IER = 1 RETURN C END OF DSOLUP END SUBROUTINE DUPUP(N,U1,U2,U3) C ---------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U1(*),U2(*),U3(*) C DO 2000 I = 1 , N IK = I * (I+1) / 2 - I DO 2000 J = 1 , I SOM = 0.0D+0 DO 1000 L = J , I IL = L * (L+1) / 2 - L SOM = SOM + U1(IL+J) * U2(IK+L) 1000 CONTINUE U3(IK+J) = SOM 2000 CONTINUE RETURN C END OF DUPUP END C SUBROUTINE DAMTES(AMP,WDIA,N,IFLAG) C ----------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WDIA(*) C AMAXU = 0.0D+0 DO 1000 I = 1 , N ABSDIA = DABS(WDIA(I)) IF (ABSDIA.GT.AMAXU) AMAXU = ABSDIA 1000 CONTINUE IF (AMAXU.GT.AMP) GOTO 1020 IFLAG = 0 RETURN 1020 IFLAG = 1 RETURN END SUBROUTINE DTAMVC(A,M,N,IM,IN,Y,X) C ---------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(IM,IN),Y(IM),X(IN) C IF (M.LE.0 .OR. N.LE.0) RETURN DO 1100 I = 1 , N SOM = 0.D0 DO 1000 J = 1 , M SOM = SOM + A(J,I) * Y(J) 1000 CONTINUE X(I) = SOM 1100 CONTINUE RETURN C END OF DTAMVC END SUBROUTINE DCDI(D,N,NTI,Q,K,WSUM) C --------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION D(N,NTI),Q(N,N,NTI),WSUM(N) C CALL DTAMVC(Q(1,1,K),N,N,N,N,D(1,K),WSUM) DO 1000 I = 1 , N D(I,K) = WSUM(I) 1000 CONTINUE RETURN END SUBROUTINE DCNRHS(DIAG,N,NRHS) C ------------------------------ C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DIAG(N) C A = DABS(DIAG(1)) NRHS = 1 DO 1000 I = 2 , N B = DABS(DIAG(I)) IF (B.LT.1.D0) B = 1.D0 / B IF (B.LE.A) GOTO 1000 A = B NRHS = I 1000 CONTINUE RETURN C END OF DCNRHS END SUBROUTINE DCPHIS(FLIN,FDIF,N,A,B,ER,K,IHOM,NRHS,WY,WI,Q,NTI,W, 1 WPHI,WPHI0,IER) C --------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ER(5),WY(N),WI(N,NTI),Q(N,N,NTI), 1 W(N,7),WPHI(N,N),WPHI0(N,N) EXTERNAL FLIN,FDIF C LOCAL ARRAYS DIMENSION HI(5),HI2(5) C C COMPUTATION OF THE FUNDAMENTAL SOLUTION KMIN1 = K - 1 NTI1 = NTI - 1 IF ((NRHS.NE.0).OR.(IHOM.NE.0)) GOTO 1600 C NRHS=0 AND IHOM=0 , DETERMONATION OF NRHS ON THE FIRST MINOR C SHOOTING INTERVAL DO 1500 L = 1 , N T1 = A T2 = B DO 1100 I = 1 , N WY(I) = Q(I,L,1) 1100 CONTINUE IER = 1 CALL DRKFSM(FLIN,N,WY,T1,T2,ER(1),ER(2),ER(3),HI2,NHI2,W,IER) IF (IER.GT.3) RETURN DEL2 = (T2 - A) / NHI2 IF (L.EQ.1) GOTO 1200 IF (DABS(DEL2).GT.DABS(DEL1)) GOTO 1500 1200 DEL1 = DEL2 NRHS = L NHI = NHI2 DO 1300 I = 1 , NHI HI(I) = HI2(I) 1300 CONTINUE T3 = T2 DO 1400 I = 1 , N WPHI(I,NRHS) = WY(I) WPHI0(I,NRHS)= WY(I) 1400 CONTINUE 1500 CONTINUE T1 = A T2 = B B = T3 GOTO 2400 1600 IF (IHOM.EQ.0) GOTO 1900 C INHOMOGENEOUS SYSTEM, COMPUTATION OF THE PARTICULAR SOLUTION DO 1700 I = 1 , N WY(I) = 0.D0 1700 CONTINUE IER = 1 T1 = A CALL DRKFSM(FDIF,N,WY,T1,B,ER(1),ER(2),ER(3),HI,NHI,W,IER) IF (IER.GT.3) RETURN DO 1800 I = 1 , N WI(I,K) = WY(I) 1800 CONTINUE NRHS = 0 GOTO 2400 1900 DO 2000 I = 1 , N WY(I) = Q(I,NRHS,KMIN1) 2000 CONTINUE IER = 1 T1 = A CALL DRKFSM(FLIN,N,WY,T1,B,ER(1),ER(2),ER(3),HI,NHI,W,IER) IF (IER.GT.3) RETURN IF (KMIN1.GT.1) GOTO 2200 DO 2100 I = 1 , N WPHI(I,NRHS) = WY(I) WPHI0(I,NRHS)= WY(I) 2100 CONTINUE GOTO 2400 2200 DO 2300 I = 1 , N WPHI(I,NRHS) = WY(I) 2300 CONTINUE GOTO 2400 2400 CONTINUE C COMPUTATION OF THE BASIC HOMOGENEOUS SOLUTIONS, WHICH HAVE NOT C BEEN COMPUTED YET DO 2900 I = 1 , N IF (I.EQ.NRHS) GOTO 2900 DO 2500 J = 1 , N WY(J) = Q(J,I,KMIN1) 2500 CONTINUE T1 = A CALL DRKFGS(FLIN,N,WY,T1,HI,NHI,W) IF (KMIN1.GT.1) GOTO 2700 DO 2600 J = 1 , N WPHI(J,I) = WY(J) WPHI0(J,I) = WY(J) 2600 CONTINUE GOTO 2900 2700 DO 2800 J = 1 , N WPHI(J,I) = WY(J) 2800 CONTINUE 2900 CONTINUE CALL DQUDEC(WPHI,N,N,N,N,W(1,1),W(1,2),WY) RETURN C END OF DCPHIS END SUBROUTINE DCROUT(A,N,M,B,EPS,P,V,IER) C -------------------------------------- C DOUBLE PRECISION A(M,*),B(*),V(*),EPS INTEGER P(*) C IER = 0 CALL DLUDEC(A,N,M,EPS,P,V,IER) IF (IER.NE.0) RETURN CALL DSOLDE(A,N,M,P,B) RETURN C END OF DCROUT C C DCROUT IS AN DOUBLE PRECISION FORTRAN VERSION OF C T.DEKKER'S SOLDEC END SUBROUTINE DINTCH(A,NA,B,NB,N) C ------------------------------ C DOUBLE PRECISION A(*),B(*),C C IF (N.LE.0) RETURN NN = 1 - NB JL = NA * N DO 1000 JJ = 1 , JL , NA NN = NN + NB C = B(NN) B(NN) = A(JJ) A(JJ) = C 1000 CONTINUE RETURN C END OF DINTCH C C DINTCH IS A DOUBLE PRECISION VERSION OF C J.M. VAN KATS SUBROUTINE INTCHA END SUBROUTINE DLUDEC(A,N,M,EPS,P,V,IER) C ------------------------------------ C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(M,*),V(*) INTEGER P(*),PK C IER = 0 R = -1.0D+0 Z = 1.0D+0 DO 1000 I = 1 , N CALL DINPRO(A(I,1),M,A(I,1),M,N,S) S = DSQRT(S) IF (S.GT.R) R = S V(I) = 1.0D+0 / S 1000 CONTINUE E = EPS * R DO 5000 K = 1 , N R = -1.0D+0 K1 = K - 1 DO 2000 I = K , N CALL DINPRO(A(I,1),M,A(1,K),1,K1,S) A(I,K) = A(I,K) - S S = DABS(A(I,K)) * V(I) IF (S.LE.R) GOTO 2000 R = S PK = I 2000 CONTINUE P(K) = PK V(PK) = V(K) Z = A(PK,K) IF (DABS(Z).LT.E) GOTO 7000 IF (PK.EQ.K) GOTO 3000 CALL DINTCH(A(K,1),M,A(PK,1),M,N) 3000 KK = K + 1 IF (K.EQ.N) RETURN DO 4000 I = KK , N CALL DINPRO(A(K,1),M,A(1,I),1,K1,S) A(K,I) = (A(K,I) - S) / Z 4000 CONTINUE 5000 CONTINUE 7000 IER = 1 RETURN C END OF DLUDEC C C DLUDEC IS AN FORTRAN VERSION OF T. DEKKER'S DEC END SUBROUTINE DSOLDE(A,N,M,P,B) C ---------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(M,*),B(*) INTEGER P(*),PK C DO 1000 K = 1 , N R = B(K) PK = P(K) CALL DINPRO(A(K,1),M,B,1,K-1,C) B(K) = (B(PK) - C) / A(K,K) IF (PK.NE.K) B(PK) = R 1000 CONTINUE IF (N.EQ.1) RETURN K = N 2000 K = K - 1 CALL DINPRO(A(K,K+1),M,B(K+1),1,N-K,C) B(K) = B(K) - C IF (K.GT.1) GOTO 2000 RETURN C END OF DSOLDE C C DSOLDE IS DOUBLE PRECISION VERSION OF H.ZWEERUS-VINK'S C AND J.M. VAN KAT'S SUBROUTINE SOLDEC END