C
C The file hrcu.for -- revised January 2001
C
C *** ESMTH ***
C (NEW VERSION - 5/10/91)
C
C THIS SUBROUTINE FITS A SPLINE CURVE G(T) IN R**3 TO THE UNWRAPPED
C DATA POINTS WSTAR(I), USING THE UNWRAPPED ERROR MATRICES SIGINV(I).
C IN REGARD TO THE SMOOTHING PARAMETER ALAMBDA, THE SUBROUTINE
C OPERATES IN ONE OF FOUR MODES (SEE EXPLANATION OF THE INPUT
C PARAMETER ISMTH BELOW).
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - TIMES
C WSTAR(3,NDATAT) - UNWRAPPED DATA POINTS
C SIGINV(3,3,NDATAT) - UNWRAPPED ERROR MATRICES
C ISMTH - CODE FOR CHOICE OF SMOOTHING PARAMETER ALAMBDA:
C 0 - USE SPECIFIED VALUE OF ALAMBDA
C 1 - DETERMINE ALAMBDA BY CROSS-VALIDATION
C 2 - DETERMINE ALAMBDA BY REQUIRING THE FITTED
C SPLINE POINTS TO LIE IN THE KNOWN CONFIDENCE
C REGIONS
C 3 - DETERMINE ALAMBDA BY THE DISCREPANCY METHOD
C ILAMB - CODE FOR VALUE OF ALAMBDA (IN CASE ISMTH = 0):
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ALAMB - VALUE OF ALAMBDA (IN CASE ISMTH = 0 AND ILAMB = 1)
C IPRINT - CODE FOR PRINTING OPTION:
C 0 - STANDARD
C 261 - THE ARRAYS GHAT AND WFIT
C 262 - THE OBJECTIVE FUNCTION AND TERMS MAKING IT UP
C ITER - ITERATION NO.
C
C OUTPUTS:
C ILAMB - CODE FOR VALUE OF ALAMBDA (FINITE OR INFINITE),
C USED AS OUTPUT IN CASE ISMTH = 1, 2, OR 3
C ALAMB - VALUE OF ALAMBDA, USED AS OUTPUT IN CASE ISMTH = 1,
C 2, OR 3 AND ILAMB = 1
C XVFVAL - CROSS-VALID. SCORE CORRESP. TO ALAMBDA
C GHAT(3,NDATAT) - THE OPTIMAL COEFFICIENTS GAMMA-HAT(K,I)
C WFIT(3,NDATAT) - THE FITTED POINTS
C FOBJ - VALUE OF THE OBJECTIVE FUNCTION AT OPT. GAMMA
C
C NOTE: THE OBJECTIVE FUNCTION IS
C
C F(GAMMA) = ALAMBDA*TRANSPOSE(GAMMA)*OMEGA*GAMMA
C + TRANSPOSE(DIFF)*SIGINV*DIFF ,
C
C WHERE DIFF = WFIT - WSTAR .
C
SUBROUTINE ESMTH(NDATA,TAU,WSTAR,SIGINV,ISMTH,ILAMB,
&ALAMB,IPRINT,ITER,XVFVAL,GHAT,WFIT,FOBJ)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),WSTAR(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT)
DIMENSION GHAT(3,NDATAT),WFIT(3,NDATAT)
DIMENSION DIFF(3,NDATAT),GVAL(3)
DIMENSION OM(0:3,NDATAT)
C
C BRANCH ON THE VALUE OF ISMTH:
C IF ISMTH = 0, THE VALUE OF ALAMBDA IS SPECIFIED. WE ONLY
C CALCULATE THE CORRESP. CROSS-VALID. SCORE BEFORE DOING THE FIT.
C IF ISMTH = 1, WE USE CROSS-VALIDATION TO DETERMINE THE OPTIMAL
C VALUE OF ALAMBDA, BEFORE DOING THE FIT.
C IF ISMTH = 2, WE READ THE VALUE OF KMISS FROM THE TERMINAL AND
C THEN DETERMINE ALAMBDA AS THE LARGEST VALUE OF THE SMOOTHING
C PARAMETER SUCH THAT ALL BUT KMISS FITTED SPLINE POINTS LIE IN
C THE CORRESP. CONFIDENCE REGIONS. WE ALSO CALCULATE THE CROSS-
C VALIDATION SCORE.
C IF ISMTH = 3, WE DETERMINE ALAMBDA AS THE LARGEST VALUE OF THE
C SMOOTHING PARAMETER SUCH THAT
C
C (SUM OF (D(ALAMBDA,I)**T)*SIGINV(I)*D(ALAMBDA,I))/(3*NDATA)
C
C IS LESS OR EQUAL 1.0. HERE D(ALAMBDA,I) IS THE DISCREPANCY AT
C THE I-TH DATA POINT, I.E.:
C D(ALAMBDA,I) = WSTAR(I) - G(ALAMBDA,TAU(I)) ,
C WHERE G(ALAMBDA,T) IS THE FITTED SPLINE.
C
IF (ISMTH.EQ.0) THEN
CALL EVXV(NDATA,TAU,WSTAR,SIGINV,ILAMB,ALAMB,IPRINT,XVFVAL,
& RCMIN,RCMAX)
ELSE IF (ISMTH.EQ.1) THEN
CALL FOLAM(NDATA,TAU,WSTAR,SIGINV,IPRINT,ITER,ILAMB,ALAMB,
& XVFVAL)
ELSE IF (ISMTH.EQ.2) THEN
PRINT *,' '
PRINT *,'SPECIFY THE NO. K OF FITTED SPLINE POINTS YOU ARE'
PRINT *,'WILLING TO ALLOW TO FALL OUTSIDE THE CORRESP.'
PRINT *,'CONFIDENCE REGIONS. (K MUST BE BETWEEN 0 AND'
PRINT *,'NDATA-1, INCLUSIVE.)'
READ *,KMISS
CALL FOLAM2(NDATA,TAU,WSTAR,SIGINV,KMISS,IPRINT,ILAMB,
& ALAMB,PHIVAL)
CALL EVXV(NDATA,TAU,WSTAR,SIGINV,ILAMB,ALAMB,IPRINT,XVFVAL,
& RCMIN,RCMAX)
ELSE IF (ISMTH.EQ.3) THEN
CALL FOLAM3(NDATA,TAU,WSTAR,SIGINV,IPRINT,ILAMB,ALAMB,
& PHI2V)
CALL EVXV(NDATA,TAU,WSTAR,SIGINV,ILAMB,ALAMB,IPRINT,XVFVAL,
& RCMIN,RCMAX)
END IF
C
C THE FITTING PROCESS.
C
TVAL=TAU(1)
CALL SEFIT(NDATA,TAU,WSTAR,SIGINV,ILAMB,ALAMB,TVAL,IPRINT,
&GVAL,GHAT,WFIT,RCOND)
C
IF ((IPRINT.EQ.261).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'THE ARRAY GHAT: GHAT(K,I) IS THE OPTIMAL COEFF. OF'
PRINT *,'THE I-TH B-SPLINE, CORRESP. TO THE K-TH COMPONENT.'
PRINT 950
DO 400 I=1,NDATA
400 PRINT 951,I,(GHAT(K,I),K=1,3)
950 FORMAT(' ',2X,'I',5X,'GHAT(K,I), K = 1, 2, 3')
951 FORMAT(' ',I3,3(2X,G12.5))
PRINT *,' '
PRINT *,'THE ARRAY WFIT'
PRINT 952
DO 410 I=1,NDATA
410 PRINT 951,I,(WFIT(K,I),K=1,3)
952 FORMAT(' ',2X,'I',5X,'WFIT(K,I), K = 1, 2, 3')
END IF
C
C CALCULATE DIFF = WFIT - WSTAR .
C
DO 310 I=1,NDATA
DO 310 K=1,3
310 DIFF(K,I)=WFIT(K,I) - WSTAR(K,I)
C
C CALCULATE FTERM2 = TRANSPOSE(DIFF)*SIGINV*DIFF .
C
FTERM2=0.0D0
DO 320 M=1,NDATA
DO 320 K=1,3
DO 320 L=1,3
320 FTERM2=FTERM2 + DIFF(K,M)*DIFF(L,M)*SIGINV(K,L,M)
C
C CALCULATE FTERM1 = ALAMBDA*TRANSPOSE(GHAT)*OMEGA*GHAT .
C
CALL CALCOM(NDATA,TAU,OM)
C
SUM1=0.0D0
DO 330 I=1,NDATA
DO 330 K=1,3
330 SUM1=SUM1 + GHAT(K,I)*GHAT(K,I)*OM(0,I)
C
SUM2=0.0D0
DO 340 I=1,NDATA-1
JMAX=MIN(NDATA,I+3)
DO 340 J=I+1,JMAX
DO 340 K=1,3
340 SUM2=SUM2 + GHAT(K,I)*GHAT(K,J)*OM(J-I,I)
C
IF (ILAMB.EQ.1) THEN
FTERM1=ALAMB*(SUM1 + 2.0D0*SUM2)
ELSE
FTERM1=0.0D0
END IF
C
C CALCULATE FOBJ.
C
FOBJ=FTERM1 + FTERM2
C
IF ((IPRINT.EQ.262).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'THE OBJECTIVE FUNCTION: FOBJ'
PRINT *,'NOTE: FTERM1 = ALAMBDA*TRANSPOSE(GHAT)*OMEGA*GHAT'
PRINT *,' FTERM2 = TRANSPOSE(DIFF)*SIGINV*DIFF,'
PRINT *,' WHERE DIFF = WFIT - WSTAR'
PRINT *,' FOBJ = FTERM1 + FTERM2'
PRINT 960
PRINT 961,FTERM1,FTERM2,FOBJ
960 FORMAT(' ',8X,'FTERM1',17X,'FTERM2',18X,'FOBJ')
961 FORMAT(' ',G22.15,2(1X,G22.15))
END IF
C
RETURN
END
C
C *** FOLAM ***
C
C THIS SUBROUTINE FINDS THE OPTIMAL VALUE OF THE SMOOTHING
C PARAMETER ALAMBDA, FOR A SET OF DATA IN R**3, USING THE
C METHOD OF CROSS-VALIDATION.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - TIMES
C Y(3,NDATAT) - DATA POINTS
C SIGINV(3,3,NDATAT) - ERROR MATRICES
C IPRINT - CODE FOR PRINTING OPTION
C 0 - STANDARD
C 251 - RHO
C 252 - TABLE OF CROSS-VALIDATION SCORES
C 253 - TABLE OF EST. RECIPS. OF COND. NOS.
C ITER - ITERATION NO.
C
C OUTPUTS:
C ILAMB - CODE FOR OPTIMAL VALUE OF ALAMBDA:
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ALAMB - VALUE OF ALAMBDA, IN CASE ILAMB = 1
C XVFVAL - CROSS-VALID. SCORE CORRESP. TO THE OPT. VALUE
C OF ALAMBDA
C
C THE SUBROUTINE FOLAM CALLS THE SUBROUTINES CLCEIG AND EVXV.
C
SUBROUTINE FOLAM(NDATA,TAU,Y,SIGINV,IPRINT,ITER,ILAMB,
&ALAMB,XVFVAL)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),Y(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT)
DIMENSION XVTAB(0:100)
DIMENSION RCMIN(0:100),RCMAX(0:100)
C
C CALCULATE RHO.
C
CALL CLCEIG(NDATA,TAU,SIGINV,IPRINT,RHO)
IF ((IPRINT.EQ.251).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT 905,RHO
905 FORMAT(' ','RHO = ',G12.5)
END IF
C
C CALCULATE TABLE OF CROSS-VALIDATION SCORES, CORRESP. TO
C DIFFERENT VALUES OF ALAMBDA.
C
PRINT *,'TYPE MMAX = GRID SIZE FOR DETERMINING OPTIMAL VALUE'
PRINT *,'OF ALAMBDA (AT MOST 100; MMAX = 10 IS USUALLY A'
PRINT *,'REASONABLE VALUE)'
910 READ *,MMAX
IF (MMAX.LT.1) THEN
PRINT *,'MMAX MUST BE AT LEAST 1. PLEASE RETYPE.'
GO TO 910
END IF
IF (MMAX.GT.100) THEN
PRINT *,'MMAX MUST BE AT MOST 100. PLEASE RETYPE.'
GO TO 910
END IF
C
GNUZ=0.0D0
GNUINC=1.0D0/DBLE(MMAX)
C
C ARRANGE PRINTING OF CROSS-VALIDATION SCORES.
C
PRINT *,' '
PRINT *,'TABLE OF CROSS-VALIDATION SCORES'
PRINT 920
920 FORMAT(' ',2X,'M',7X,'GNU',10X,'ALAMB',12X,'XVTAB(M)')
C
DO 110 M=0,MMAX
GNU=GNUZ + DBLE(M)*GNUINC
IF (ABS(1.0D0 - GNU).LT.1.0D-12) THEN
ILAMB=2
ELSE
ILAMB=1
ALAMB=(GNU/(1.0D0 - GNU))*(1.0D0/RHO)
END IF
CALL EVXV(NDATA,TAU,Y,SIGINV,ILAMB,ALAMB,IPRINT,XVTAB(M),
&RCMIN(M),RCMAX(M))
C
IF (ILAMB.EQ.1) THEN
PRINT 930,M,GNU,ALAMB,XVTAB(M)
ELSE
PRINT 931,M,GNU,XVTAB(M)
END IF
930 FORMAT(' ',I3,2X,G12.5,2X,G12.5,2X,G22.15)
931 FORMAT(' ',I3,2X,G12.5,2X,' INFINITY ',2X,G22.15)
C
110 CONTINUE
C
IF ((IPRINT.EQ.253).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'TABLE OF EST. RECIPROCALS OF COND. NOS.'
PRINT *,'(MIN. AND MAX. OVER IDEL)'
PRINT 940
DO 120 M=0,MMAX
120 PRINT 941,M,RCMIN(M),RCMAX(M)
940 FORMAT(' ',2X,'M',6X,'RCMIN',9X,'RCMAX')
941 FORMAT(' ',I3,2X,G12.5,2X,G12.5)
END IF
C
C FIND OPTIMAL VALUE OF ALAMBDA, CORRESP. TO MIN. CROSS-
C VALIDATION SCORE.
C
MSTAR=0
XVFVAL=XVTAB(0)
DO 130 M=1,MMAX
IF (XVTAB(M).LT.XVFVAL) THEN
MSTAR=M
XVFVAL=XVTAB(M)
END IF
130 CONTINUE
C
GNUSTR=GNUZ + DBLE(MSTAR)*GNUINC
IF (ABS(1.0D0 - GNUSTR).LT.1.0D-12) THEN
ILAMB=2
ELSE
ILAMB=1
ALAMB=(GNUSTR/(1.0D0 - GNUSTR))*(1.0D0/RHO)
END IF
C
C PRINT THE OPTIMAL VALUE OF ALAMBDA.
C
PRINT *,' '
IF (ILAMB.EQ.1) THEN
PRINT 950,ALAMB
950 FORMAT(' ','OPT. VALUE OF ALAMBDA = ',G12.5)
ELSE
PRINT *,'OPT. VALUE OF ALAMBDA = INFINITY.'
END IF
C
RETURN
END
C
C *** EVXV ***
C
C SUBROUTINE TO EVALUATE THE CROSS-VALIDATION SCORE.
C
C INPUTS:
C NDATA - N0. OF DATA POINTS
C TAU(NDATAT) - TIMES
C Y(3,NDATAT) - DATA POINTS
C SIGINV(3,3,NDATAT) - ERROR MATRICES
C ILAMB - CODE FOR VALUE OF THE SMOOTHING PARAMETER ALAMBDA
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ALAMB - VALUE OF ALAMBDA, IN THE CASE ILAMB = 1
C IPRINT - PRINTING PARAMETER
C 221 - DETAILS OF CALC. OF XVSC
C
C OUTPUTS:
C XVSC - CROSS-VALIDATION SCORE
C RCMIN - MIN. OF RCOND (EST. OF RECIP. OF COND. NO.
C OF COEFF. MATRIX), OVER ALL VALUES OF IDEL
C RCMAX - MAX. OF RCOND, OVER ALL VALUES OF IDEL
C
SUBROUTINE EVXV(NDATA,TAU,Y,SIGINV,ILAMB,ALAMB,IPRINT,XVSC,
& RCMIN,RCMAX)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),Y(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT),GVAL(3,NDATAT)
DIMENSION GHAT(3,NDATAT),YFIT(3,NDATAT)
DIMENSION RTAU(NDATAT),RY(3,NDATAT)
DIMENSION RSIGIV(3,3,NDATAT),R(3,NDATAT)
DIMENSION TEMP(3),TERM(NDATAT),RCOND(NDATAT)
C
XVSC=0.0D0
C
C MAIN LOOP.
C
DO 100 IDEL=1,NDATA
TVAL=TAU(IDEL)
C
ITARG=0
DO 50 I=1,NDATA
IF (I.NE.IDEL) THEN
ITARG=ITARG + 1
RTAU(ITARG)=TAU(I)
DO 55 K=1,3
55 RY(K,ITARG)=Y(K,I)
DO 57 K=1,3
DO 57 L=1,3
57 RSIGIV(K,L,ITARG)=SIGINV(K,L,I)
END IF
50 CONTINUE
C
CALL SEFIT(NDATA-1,RTAU,RY,RSIGIV,ILAMB,ALAMB,TVAL,IPRINT,
& GVAL(1,IDEL),GHAT,YFIT,RCOND(IDEL))
C
DO 60 K=1,3
60 R(K,IDEL)=GVAL(K,IDEL) - Y(K,IDEL)
C
DO 70 K=1,3
TEMP(K)=0.0D0
DO 65 L=1,3
65 TEMP(K)=TEMP(K) + SIGINV(K,L,IDEL)*R(L,IDEL)
70 CONTINUE
C
TERM(IDEL)=0.0D0
DO 75 K=1,3
75 TERM(IDEL)=TERM(IDEL) + R(K,IDEL)*TEMP(K)
C
XVSC=XVSC + TERM(IDEL)
IF (IDEL.EQ.1) THEN
RCMIN=RCOND(IDEL)
RCMAX=RCOND(IDEL)
ELSE
RCMIN=MIN(RCMIN,RCOND(IDEL))
RCMAX=MAX(RCMAX,RCOND(IDEL))
END IF
100 CONTINUE
C
C END OF MAIN LOOP.
C
IF (IPRINT.EQ.221) THEN
PRINT *,' '
IF (ILAMB.EQ.1) THEN
PRINT 905,ALAMB
ELSE
PRINT *,'ILAMB = 2; ALAMBDA = INFINITY'
END IF
PRINT 910,XVSC
905 FORMAT(' ','ILAMB = 1; ALAMB = ',G12.5)
910 FORMAT(' ','CROSS-VALID. SCORE = ',G22.15)
C
PRINT *,' '
PRINT *,'DETAILS OF THE CALC. OF THE CROSS-VALIDATION'
PRINT *,'SCORE. FOR EACH DATA POINT, WE GIVE:'
PRINT *,' DATA VECTOR'
PRINT *,' FITTED DATA VECTOR'
PRINT *,' COMPONENT OF XVSC DUE TO THIS DATA POINT;'
PRINT *,' EST. OF RECIP. OF COND. NO. OF COEFF.'
PRINT *,' MATRIX'
C
DO 120 I=1,NDATA
PRINT *,' '
PRINT *,'DATA POINT NO. I = ',I
PRINT 911,(Y(K,I),K=1,3)
PRINT 911,(GVAL(K,I),K=1,3)
PRINT 912,TERM(I),RCOND(I)
120 CONTINUE
911 FORMAT(' ',3(2X,G12.5))
912 FORMAT(' ','COMP. OF XVSC = ',G22.15,
& ' ;RCOND = ',G12.5)
END IF
C
RETURN
END
C
C *** SEFIT ***
C
C THIS SUBROUTINE FITS A SMOOTHING SPLINE IN R**3 TO A SET OF DATA
C POINTS HAVING ELLIPTICAL ERROR DISTRIBUTIONS.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - TIMES
C Y(3,NDATAT) - DATA POINTS
C SIGINV(3,3,NDATAT) - ERROR MATRICES
C ILAMB - CODE FOR VALUE OF THE SMOOTHING PARAMETER ALAMBDA
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ALAMB - VALUE OF ALAMBDA, IN THE CASE ILAMB = 1
C TVAL - VALUE OF T AT WHICH TO EVALUATE THE FITTED
C SPLINE
C IPRINT - PRINTING PARAMETER
C 0 - NO PRINTING
C 201 - THE MATRIX OMEGA
C 202 - THE ARRAYS BMAT, ETAU
C 203 - THE ARRAY GMAT
C 204 - THE ARRAYS Z, H
C 205 - THE ARRAY SMAT
C 206 - THE EIGENVALUES AND EIGENVECTORS
C OF THE MATRIX OMEGA-S
C 207 - THE VECTOR DELTA
C 208 - THE MATRIX ALPHA
C 209 - THE SOLUTION VECTOR EPS
C 210 - THE ARRAY GHAT
C 211 - THE ARRAY YFIT
C
C OUTPUTS:
C GVAL(3) - VALUE OF THE FITTED SPLINE AT T = TVAL
C GHAT(3,NDATAT) - THE OPTIMAL COEFFS. GAMMA-HAT(K,I)
C YFIT(3,NDATAT) - THE FITTED POINTS
C RCOND - EST. OF THE RECIPROCAL OF THE CONDITION
C NUMBER OF THE COEFF. MATRIX
C
SUBROUTINE SEFIT(NDATA,TAU,Y,SIGINV,ILAMB,ALAMB,TVAL,IPRINT,
& GVAL,GHAT,YFIT,RCOND)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),Y(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT),GVAL(3)
DIMENSION OM(0:3,NDATAT),BMAT(-1:1,NDATAT)
DIMENSION GMAT(3,3,-2:2,NDATAT)
DIMENSION Z(3,NDATAT),H(3,NDATAT)
DIMENSION GHAT(3,NDATAT),YFIT(3,NDATAT)
DIMENSION BETA(NDATAT)
DIMENSION ETAU(-2:NDATAT+3)
DIMENSION SMAT(10,3*NDATAT),C(3*NDATAT)
DIMENSION ZWORK(3*NDATAT)
DIMENSION ALPHA(6,6),DELTA(6),EPS(6)
DIMENSION AMAT(6,6),ZAWORK(6)
DIMENSION OMSMAT(NDATAT,4),OMEVAL(NDATAT)
DIMENSION OMEVEC(NDATAT,NDATAT)
DIMENSION EOM(NDATAT),EOM2(NDATAT)
DIMENSION XI(NDATAT),ETA(NDATAT)
LOGICAL MATZ
C
C CALCULATE THE MATRIX OMEGA.
C
CALL CALCOM(NDATA,TAU,OM)
IF (IPRINT.EQ.201) THEN
PRINT *,' '
PRINT *,'THE DIAGONALS OF THE MATRIX OMEGA'
PRINT 900
DO 400 I=1,NDATA-3
400 PRINT 901,I,OM(0,I),OM(1,I),OM(2,I),OM(3,I)
DO 410 I=NDATA-2,NDATA
410 PRINT 901,I,(OM(K,I),K=0,NDATA-I)
900 FORMAT(' ',2X,'I',4X,'OMEGA(I,I)',2X,'OMEGA(I,I+1)',
& 2X,'OMEGA(I,I+2)',2X,'OMEGA(I,I+3)')
901 FORMAT(' ',I3,4(2X,G12.5))
END IF
C
C CALCULATE THE MATRIX B AND THE ARRAY ETAU.
C
CALL CBMAT(NDATA,TAU,BMAT,ETAU)
IF (IPRINT.EQ.202) THEN
PRINT *,' '
PRINT *,'THE ARRAY BMAT:'
PRINT *,' BMAT(-1,I) CONTAINS B(I-1,I) = BETA(I,TAU(I-1))'
PRINT *,' BMAT(0,I) CONTAINS B(I,I) = BETA(I,TAU(I))'
PRINT *,' BMAT(1,I) CONTAINS B(I+1,I) = BETA(I,TAU(I+1))'
PRINT 922
PRINT 923,1,BMAT(0,1),BMAT(1,1)
DO 440 I=2,NDATA-1
440 PRINT 924,I,(BMAT(J,I),J= -1,1)
PRINT 925,NDATA,BMAT(-1,NDATA),BMAT(0,NDATA)
922 FORMAT(' ',2X,'I',3X,'BMAT(-1,I)',5X,'BMAT(0,I)',
& 5X,'BMAT(1,I)')
923 FORMAT(' ',I3,14X,2(2X,G12.5))
924 FORMAT(' ',I3,3(2X,G12.5))
925 FORMAT(' ',I3,2(2X,G12.5))
PRINT *,' '
PRINT *,'THE ARRAY ETAU:'
PRINT 930
DO 450 I=-2,NDATA+3
450 PRINT 931,I,ETAU(I)
930 FORMAT(' ',3X,'I',6X,'ETAU(I)')
931 FORMAT(' ',I4,2X,G12.5)
END IF
C
C FOR CONVENIENCE, WE SET N = 3*NDATA.
C
N=3*NDATA
C
C CALCULATE THE ARRAY GMAT WHICH CONTAINS THE MATRIX G:
C
C G = TRANSPOSE(B)*SIGINV*B .
C
C SUPPOSE WE CONSIDER G TO BE A BLOCK-MATRIX, WITH 3X3 SUB-BLOCKS
C G(I,J). WE DENOTE THE (K,L)-TH ELEMENT OF G(I,J) BY G(K,L,I,J).
C THEN:
C
C GMAT(K,L,0,I) CONTAINS G(K,L,I,I) ,
C GMAT(K,L,1,I) CONTAINS G(K,L,I,I+1) ,
C GMAT(K,L,2,I) CONTAINS G(K,L,I,I+2) .
C
C BRIEFLY STATED, G(I,J) IS THE SUM OVER M OF THE TERMS:
C
C B(M,I)*B(M,J)*SIGINV(M) .
C
C FIRST WE CALCULATE THE SUB-BLOCKS G(I,I).
C
DO 20 I=1,NDATA
MSTART=MAX(1,I-1)
MSTOP=MIN(NDATA,I+1)
DO 15 K=1,3
DO 15 L=1,3
GMAT(K,L,0,I)=0.0D0
DO 10 M=MSTART,MSTOP
TERM=BMAT(M-I,I)*BMAT(M-I,I)*SIGINV(K,L,M)
10 GMAT(K,L,0,I)=GMAT(K,L,0,I) + TERM
15 CONTINUE
20 CONTINUE
C
C NEXT WE CALCULATE THE SUB-BLOCKS G(I,I+1).
C
DO 40 I=1,NDATA-1
DO 35 K=1,3
DO 35 L=1,3
GMAT(K,L,1,I)=0.0D0
DO 30 M=I,I+1
TERM=BMAT(M-I,I)*BMAT(M-I-1,I+1)*SIGINV(K,L,M)
30 GMAT(K,L,1,I)=GMAT(K,L,1,I) + TERM
35 CONTINUE
40 CONTINUE
C
C FINALLY WE CALCULATE THE SUB-BLOCKS G(I,I+2).
C
DO 60 I=1,NDATA-2
DO 55 K=1,3
DO 55 L=1,3
GMAT(K,L,2,I)=BMAT(1,I)*BMAT(-1,I+2)*SIGINV(K,L,I+1)
55 CONTINUE
60 CONTINUE
C
C NOW WE FILL OUT THE PART OF THE ARRAY GMAT WHICH CONTAINS THE
C SUB-BLOCKS G(I,J) WITH I GREATER THAN J:
C
C GMAT(K,L,-2,I) CONTAINS G(K,L,I,I-2) ,
C GMAT(K,L,-1,I) CONTAINS G(K,L,I,I-1) .
C
DO 64 I=2,NDATA
DO 62 K=1,3
DO 62 L=1,3
GMAT(K,L,-1,I)=GMAT(L,K,1,I-1)
62 CONTINUE
64 CONTINUE
C
DO 68 I=3,NDATA
DO 66 K=1,3
DO 66 L=1,3
GMAT(K,L,-2,I)=GMAT(L,K,2,I-2)
66 CONTINUE
68 CONTINUE
C
IF (IPRINT.EQ.203) THEN
PRINT *,' '
PRINT *,'THE ARRAY GMAT: GMAT CONTAINS THE (3*NDATA)X(3*NDATA)'
PRINT *,'MATRIX G = TRANSPOSE(B)*SIGINV*B .'
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I)'
DO 460 I=1,NDATA
PRINT *,'I = ',I
DO 465 K=1,3
465 PRINT 940,(GMAT(K,L,0,I),L=1,3)
460 CONTINUE
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I+1)'
DO 470 I=1,NDATA-1
PRINT *,'I = ',I
DO 475 K=1,3
475 PRINT 940,(GMAT(K,L,1,I),L=1,3)
470 CONTINUE
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I+2)'
DO 480 I=1,NDATA-2
PRINT *,'I = ',I
DO 485 K=1,3
485 PRINT 940,(GMAT(K,L,2,I),L=1,3)
480 CONTINUE
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I-1)'
DO 488 I=2,NDATA
PRINT *,'I = ',I
DO 486 K=1,3
486 PRINT 940,(GMAT(K,L,-1,I),L=1,3)
488 CONTINUE
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I-2)'
DO 492 I=3,NDATA
PRINT *,'I = ',I
DO 490 K=1,3
490 PRINT 940,(GMAT(K,L,-2,I),L=1,3)
492 CONTINUE
940 FORMAT(' ',3(2X,G12.5))
END IF
C
C CALCULATE THE ARRAY Z(3,NDATAT). THE VECTOR Z(M) IS GIVEN BY:
C
C Z(M) = SIGINV(M)*Y(M) .
C
DO 80 M=1,NDATA
DO 75 K=1,3
Z(K,M)=0.0D0
DO 70 L=1,3
70 Z(K,M)=Z(K,M) + SIGINV(K,L,M)*Y(L,M)
75 CONTINUE
80 CONTINUE
C
C CALCULATE THE ARRAY H(3,NDATAT). THE VECTOR H(I) IS GIVEN BY:
C
C H(I) = SUM OVER M OF THE TERMS B(M,I)*Z(M) .
C
DO 100 I=1,NDATA
MSTART=MAX(1,I-1)
MSTOP=MIN(NDATA,I+1)
DO 95 K=1,3
H(K,I)=0.0D0
DO 90 M=MSTART,MSTOP
90 H(K,I)=H(K,I) + BMAT(M-I,I)*Z(K,M)
95 CONTINUE
100 CONTINUE
C
IF (IPRINT.EQ.204) THEN
PRINT *,' '
PRINT *,'THE ARRAY Z: THE VECTOR Z(M) = SIGINV(M)*Y(M)'
PRINT 950
DO 495 M=1,NDATA
495 PRINT 951,M,(Z(K,M),K=1,3)
950 FORMAT(' ',2X,'M',6X,'Z(K,M), K = 1, 2, 3')
951 FORMAT(' ',I3,1X,3(2X,G12.5))
PRINT *,' '
PRINT *,'THE ARRAY H: THE VECTOR H(I) = THE SUM OVER M'
PRINT *,'OF THE TERMS B(M,I)*Z(M) .'
PRINT 952
DO 500 I=1,NDATA
500 PRINT 951,I,(H(K,I),K=1,3)
952 FORMAT(' ',2X,'I',6X,'H(K,I), K = 1, 2, 3')
END IF
C
C BRANCH DEPENDING ON THE VALUE OF ILAMB:
C IF ILAMB = 1, WE SET UP THE MATRIX S = ALAMB*OMEGA + G , AND
C DETERMINE GHAT AS THE SOLUTION OF THE LINEAR SYSTEM
C S*GHAT = H .
C IF ILAMB = 2, WE CALCULATE THE EIGENVECTORS OF THE MATRIX
C OMEGA-S. THESE EIGENVECTORS, TOGETHER WITH THE MATRIX G AND
C THE VECTOR H, DETERMINE A LINEAR SYSTEM: ALPHA*EPS = DELTA.
C THE VECTOR GHAT CAN BE DETERMINED FROM EPS AND THE EIGENVECTORS
C OF OMEGA-S. (NOTE: IN THE CASE ILAMB = 2, THE FITTED SPLINE
C IS A STRAIGHT LINE.)
C
IF (ILAMB.EQ.2) GO TO 105
C
C THE CASE ILAMB = 1. IN THIS CASE THE SMOOTHING PARAMETER ALAMBDA
C IS FINITE.
C
C SET UP THE ARRAY SMAT, SO THAT IT CONTAINS THE MATRIX S IN BAND
C SYMMETRIC STORAGE MODE (LINPACK VERSION). NOTE THAT BAND SYMMETRIC
C STORAGE MODE IN LINPACK IS ESSENTIALLY THE TRANSPOSE OF BAND
C SYMMETRIC STORAGE MODE IN EISPACK. SPECIFICALLY, THE MAIN DIAGONAL
C OF S GOES INTO ROW 10 OF SMAT. THE 1-ST, 2-ND, ... , 9-TH SUPER-
C DIAGONALS OF S GO INTO ROWS 9, 8, ... , 1 RESP. OF SMAT, FILLED
C OUT WITH ZEROS ON THE LEFT.
C
C WE FIRST INITIALIZE SMAT TO ZERO.
C
DO 705 MROW=1,10
DO 705 MCOL=1,N
705 SMAT(MROW,MCOL)=0.0D0
C
C NEXT WE FILL IN THE NON-ZERO ENTRIES IN SMAT. THIS IS DONE BY
C CONSIDERING THE APPROPRIATE 3X3 SUB-BLOCKS S(I,J) IN TURN, AND
C SETTING UP THAT PART OF SMAT WHICH CORRESPONDS TO S(I,J). THE
C APPROPRIATE SUB-BLOCKS S(I,J) ARE THOSE SUCH THAT J RANGES FROM I
C TO MIN(I+3,NDATA), WHILE I RANGES FROM 1 TO NDATA.
C
DO 740 I=1,NDATA
JMAX=MIN(I+3,NDATA)
DO 735 J=I,JMAX
C
C DENOTE THE (K,L)-TH ENTRY OF S(I,J) BY S(K,L,I,J). THEN THE ROW AND
C COLUMN INDICES OF S(K,L,I,J) IN THE BIG MATRIX S (I.E., THE (3*NDATA)X
C (3*NDATA) MATRIX S, CONSIDERED AS ONE BIG MATRIX) ARE:
C IROW = 3*(I - 1) + K ,
C ICOL = 3*(J - 1) + L .
C WE ARE INTERESTED ONLY IN ELEMENTS ON THE DIAGONAL OR ON ONE OF THE
C FIRST 9 SUPER-DIAGONALS OF S. THEREFORE, WE REQUIRE ICOL TO LIE
C BETWEEN IROW AND IROW + 9.
C
DO 730 K=1,3
DO 730 L=1,3
IROW=3*(I-1) + K
ICOL=3*(J-1) + L
IF ((ICOL.LT.IROW).OR.(ICOL.GT.IROW+9)) GO TO 720
C
C THE ROW AND COLUMN INDICES OF THE POSITION IN THE ARRAY SMAT WHICH
C IS TO RECEIVE S(K,L,I,J) ARE:
C MROW = 10 - (ICOL - IROW) ,
C MCOL = ICOL .
C
C TO OBTAIN S(K,L,I,J), WE USE THE FORMULA:
C
C S(K,L,I,J) = ALAMBDA*OMEGA(K,L,I,J) + G(K,L,I,J) .
C
C IF K.NE.L, OMEGA(K,L,I,J) = 0. OTHERWISE, OMEGA(K,L,I,J) IS CONTAINED
C IN OM(J-I,I). IF J.EQ.I+3, G(K,L,I,J) = 0. OTHERWISE, G(K,L,I,J) IS
C CONTAINED IN GMAT(K,L,J-I,I).
C
MROW=10 - (ICOL - IROW)
MCOL=ICOL
IF (K.EQ.L) THEN
SMAT(MROW,MCOL)=ALAMB*OM(J-I,I)
END IF
IF (J.LE.I+2) THEN
SMAT(MROW,MCOL)=SMAT(MROW,MCOL) + GMAT(K,L,J-I,I)
END IF
720 CONTINUE
730 CONTINUE
C
C END OF LOOP ON K AND L.
C
735 CONTINUE
740 CONTINUE
C
C END OF LOOP ON I AND J.
C
IF (IPRINT.EQ.205) THEN
PRINT *,' '
PRINT *,'THE ARRAY SMAT: SMAT CONTAINS THE MATRIX S,'
PRINT *,'IN BAND SYMMETRIC STORAGE MODE (LINPACK VERS.).'
PRINT *,'ROW 10 OF SMAT CONTAINS THE DIAGONAL OF S,'
PRINT *,'ROW 9 CONTAINS THE FIRST SUPER-DIAGONAL, ETC.'
PRINT *,'(THE ROWS FOR THE SUPER-DIAGONALS ARE FILLED OUT'
PRINT *,'WITH ZEROS AT THE FRONT.)'
PRINT *,'ROWS 1 THROUGH 5, PRINTED VERTICALLY'
DO 800 J=1,N
800 PRINT 1050,J,(SMAT(I,J),I=1,5)
1050 FORMAT(' ',I3,1X,5(1X,G12.5))
PRINT *,'ROWS 6 THROUGH 10, PRINTED VERTICALLY'
DO 805 J=1,N
805 PRINT 1050,J,(SMAT(I,J),I=6,10)
END IF
C
C SET UP THE ARRAY C.
C
DO 750 I=1,NDATA
DO 750 K=1,3
MROW=3*(I-1) + K
C(MROW)=H(K,I)
750 CONTINUE
C
C SOLVE THE LINEAR SYSTEM: S*X = C. TO DO THIS, WE USE THE ROUTINES
C DPBCO AND DPBSL FROM LINPACK.
C LDS = ROW DIMENSION OF THE ARRAY SMAT = 10
C N = ORDER OF S = 3*NDATA
C MDIAG = NO. OF SUPERDIAGONALS = 9
C
LDS=10
MDIAG=9
CALL DPBCO(SMAT,LDS,N,MDIAG,RCOND,ZWORK,INFO)
IF (INFO.NE.0) THEN
PRINT *,' '
PRINT *,'ERROR IN CALCULATING THE FACTORIZATION OF THE'
PRINT *,'MATRIX S IN SEFIT, USING DPBCO'
PRINT *,'INFO = ',INFO
STOP
END IF
CALL DPBSL(SMAT,LDS,N,MDIAG,C)
C
C SET UP THE ARRAY GHAT(3,NDATAT), BASED ON THE SOLUTION VECTOR
C WHICH IS STORED IN C.
C
DO 760 I=1,NDATA
DO 760 K=1,3
MROW=3*(I-1) + K
GHAT(K,I)=C(MROW)
760 CONTINUE
C
C END OF THE CASE ILAMB = 1. WE NOW TRANSFER BACK TO THE MAINLINE
C FLOW.
C
GO TO 220
C
105 CONTINUE
C
C THE CASE ILAMB = 2. IN THIS CASE THE SMOOTHING PARAMETER ALAMBDA
C IS INFINITE. THE FITTED SPLINE IS A STRAIGHT LINE.
C
C CALCULATE THE EIGENVALUES AND EIGENVECTORS OF THE MATRIX
C OMEGA-S. THIS IS DONE IN ORDER TO GET A BASIS FOR THE
C NULLSPACE OF THE MATRIX OMEGA.
C
C SET UP THE ARRAY OMSMAT, WHICH CONTAINS THE MATRIX OMEGA-S IN
C BAND SYMMETRIC STORAGE MODE (EISPACK VERSION). THE 4-TH COL. OF
C OMSMAT CONTAINS THE MAIN DIAGONAL OF OMEGA-S. COLS. 3, 2, 1 OF
C OMSMAT CONTAIN THE 1-ST, 2-ND, 3-RD SUPERDIAGONALS OF OMEGA-S,
C RESPECTIVELY, FILLED IN BY ZEROS ON THE TOP.
C
DO 110 I=1,NDATA
110 OMSMAT(I,4)=OM(0,I)
C
DO 115 I=2,NDATA
115 OMSMAT(I,3)=OM(1,I-1)
OMSMAT(1,3)=0.0D0
C
DO 120 I=3,NDATA
120 OMSMAT(I,2)=OM(2,I-2)
OMSMAT(1,2)=0.0D0
OMSMAT(2,2)=0.0D0
C
DO 125 I=4,NDATA
125 OMSMAT(I,1)=OM(3,I-3)
OMSMAT(1,1)=0.0D0
OMSMAT(2,1)=0.0D0
OMSMAT(3,1)=0.0D0
C
C CALCULATE THE EIGENVALUES AND EIGENVECTORS OF THE MATRIX
C OMEGA-S: TO DO THIS, WE USE THE ROUTINES BANDR AND IMTQL2
C FROM EISPACK.
C NM = ROW DIMENSION OF OMSMAT = NDATAT
C NEIG = ORDER OF OMEGA-S = NDATA
C MB = NO. OF SUPERDIAGONALS + 1 = 4
C MATZ = .TRUE. (CALCULATE THE ORTH. TRANSF. MATRIX)
C
NM=NDATAT
NEIG=NDATA
MB=4
MATZ=.TRUE.
CALL BANDR(NM,NEIG,MB,OMSMAT,OMEVAL,EOM,EOM2,MATZ,OMEVEC)
CALL IMTQL2(NM,NEIG,OMEVAL,EOM,OMEVEC,IERR)
IF (IERR.NE.0) THEN
PRINT *,' '
PRINT *,'ERROR IN FINDING THE EIGENVALUES AND EIGENVECTORS'
PRINT *,'OF OMEGA-S IN SEFIT, USING IMTQL2'
PRINT *,'IERR = ',IERR
STOP
END IF
C
C THE VECTOR OMEVAL CONTAINS THE EIGENVALUES OF OMEGA-S IN ASCENDING
C ORDER. THE MATRIX OMEVEC CONTAINS THE CORRESP. ORTHONORMAL EIGEN-
C VECTORS AS ITS COLUMNS.
C
DO 130 I=1,NDATA
XI(I)=OMEVEC(I,1)
130 ETA(I)=OMEVEC(I,2)
C
IF (IPRINT.EQ.206) THEN
PRINT *,' '
PRINT *,'THE FIRST FOUR EIGENVALUES OF OMEGA-S,'
PRINT *,'TOGETHER WITH THEIR EIGENVECTORS'
PRINT *,' '
PRINT *,'EIGENVALUES:'
PRINT *,' '
PRINT 955,(OMEVAL(I),I=1,4)
955 FORMAT(' ',4(2X,G12.5))
PRINT *,' '
PRINT *,'EIGENVECTORS (COLUMNS):'
PRINT *,' '
DO 510 K=1,NDATA
510 PRINT 955,(OMEVEC(K,I),I=1,4)
PRINT *,' '
PRINT *,'NOTE THAT THE FIRST TWO EIGENVECTORS OF'
PRINT *,'OMEGA-S BECOME THE VECTORS XI AND ETA.'
END IF
C
C CALCULATE THE VECTOR DELTA.
C
DO 140 L=1,3
SUM=0.0D0
DO 145 K=1,NDATA
145 SUM=SUM + XI(K)*H(L,K)
140 DELTA(L)=SUM
C
DO 150 L=4,6
SUM=0.0D0
DO 155 K=1,NDATA
155 SUM=SUM + ETA(K)*H(L-3,K)
150 DELTA(L)=SUM
C
IF (IPRINT.EQ.207) THEN
PRINT *,' '
PRINT *,'THE VECTOR DELTA'
DO 520 L=1,6
520 PRINT 960,DELTA(L)
960 FORMAT(' ',2X,G12.5)
END IF
C
C CALCULATE THE MATRIX ALPHA.
C
DO 160 L=1,3
DO 160 K=L,3
SUM=0.0D0
DO 165 I=1,NDATA
JSTART=MAX(1,I-2)
JSTOP=MIN(NDATA,I+2)
DO 165 J=JSTART,JSTOP
165 SUM=SUM + XI(I)*XI(J)*GMAT(L,K,J-I,I)
160 ALPHA(L,K)=SUM
C
DO 170 L=4,6
DO 170 K=L,6
SUM=0.0D0
DO 175 I=1,NDATA
JSTART=MAX(1,I-2)
JSTOP=MIN(NDATA,I+2)
DO 175 J=JSTART,JSTOP
175 SUM=SUM + ETA(I)*ETA(J)*GMAT(L-3,K-3,J-I,I)
170 ALPHA(L,K)=SUM
C
DO 180 L=1,3
DO 180 K=4,6
SUM=0.0D0
DO 185 I=1,NDATA
JSTART=MAX(1,I-2)
JSTOP=MIN(NDATA,I+2)
DO 185 J=JSTART,JSTOP
185 SUM=SUM + XI(I)*ETA(J)*GMAT(L,K-3,J-I,I)
180 ALPHA(L,K)=SUM
C
DO 190 L=2,6
DO 190 K=1,L-1
190 ALPHA(L,K)=ALPHA(K,L)
C
IF (IPRINT.EQ.208) THEN
PRINT *,' '
PRINT *,'THE MATRIX ALPHA'
DO 530 L=1,6
PRINT 965,(ALPHA(L,K),K=1,3)
530 PRINT 966,(ALPHA(L,K),K=4,6)
965 FORMAT(' ',3(2X,G12.5))
966 FORMAT(' ',2X,3(2X,G12.5))
END IF
C
C SET UP THE ARRAY AMAT, SO THAT IT CONTAINS THE MATRIX ALPHA IN
C BAND SYMMETRIC STORAGE MODE (LINPACK VERSION). THE MAIN DIAGONAL
C OF ALPHA GOES INTO ROW 6 OF AMAT. THE 1-ST, 2-ND, ... , 5-TH
C SUPERDIAGONALS OF ALPHA GO INTO ROWS 5, 4, ... , 1 RESP. OF AMAT,
C FILLED OUT WITH ZEROS ON THE LEFT.
C
C FIRST INITIALIZE AMAT TO ZERO.
C
DO 192 K=1,6
DO 192 I=1,6
192 AMAT(K,I)=0.0D0
C
C NEXT TRANSFER THE DIAGONALS OF ALPHA.
C
DO 196 K=0,5
DO 194 I=K+1,6
194 AMAT(6-K,I)=ALPHA(I-K,I)
196 CONTINUE
C
C SOLVE THE SYSTEM: ALPHA*EPS = DELTA. TO DO THIS, WE USE THE
C ROUTINES DPBCO AND DPBSL FROM LINPACK.
C LDA = ROW DIMENSION OF AMAT = 6
C NEQU = ORDER OF ALPHA = 6
C MDIAG = NO. OF SUPERDIAGONALS = 5
C
LDA=6
NEQU=6
MDIAG=5
CALL DPBCO(AMAT,LDA,NEQU,MDIAG,RCOND,ZAWORK,INFO)
IF (INFO.NE.0) THEN
PRINT *,' '
PRINT *,'ERROR IN CALCULATING THE FACTORIZATION OF THE'
PRINT *,'MATRIX ALPHA IN SEFIT, USING DPBCO'
PRINT *,'INFO = ',INFO
STOP
END IF
CALL DPBSL(AMAT,LDA,NEQU,MDIAG,DELTA)
DO 198 I=1,6
198 EPS(I)=DELTA(I)
C
IF (IPRINT.EQ.209) THEN
PRINT *,' '
PRINT *,'THE SOLUTION VECTOR EPS'
DO 540 K=1,6
540 PRINT 970,EPS(K)
970 FORMAT(' ',2X,G12.5)
END IF
C
C CALCULATE THE ARRAY GHAT.
C
DO 200 I=1,NDATA
DO 200 K=1,3
200 GHAT(K,I)=XI(I)*EPS(K) + ETA(I)*EPS(K+3)
C
C END OF THE CASE ILAMB = 2. NOW RETURN TO THE MAINLINE FLOW.
C
220 CONTINUE
IF (IPRINT.EQ.210) THEN
PRINT *,' '
PRINT *,'THE ARRAY GHAT: GHAT(K,I) IS THE OPTIMAL COEFF. OF'
PRINT *,'THE I-TH B-SPLINE, CORRESP. TO THE K-TH COMPONENT.'
PRINT 980
DO 550 I=1,NDATA
550 PRINT 981,I,(GHAT(K,I),K=1,3)
980 FORMAT(' ',2X,'I',5X,'GHAT(K,I), K = 1, 2, 3')
981 FORMAT(' ',I3,3(2X,G12.5))
END IF
C
C CALCULATE THE ARRAY YFIT = B* GHAT.
C
DO 230 I=1,NDATA
JSTART=MAX(1,I-1)
JSTOP=MIN(NDATA,I+1)
DO 230 K=1,3
YFIT(K,I)=0.0D0
DO 235 J=JSTART,JSTOP
235 YFIT(K,I)=YFIT(K,I) + BMAT(I-J,J)*GHAT(K,J)
230 CONTINUE
C
IF (IPRINT.EQ.211) THEN
PRINT *,' '
PRINT *,'THE ARRAY YFIT = B*GHAT:'
PRINT 985
DO 560 I=1,NDATA
560 PRINT 986,I,(YFIT(K,I),K=1,3)
985 FORMAT(' ',2X,'I',6X,'YFIT(K,I), K = 1, 2, 3')
986 FORMAT(' ',I3,1X,3(2X,G12.5))
END IF
C
C CALCULATE THE VALUE OF THE SMOOTHING SPLINE G(T) AT THE GIVEN
C VALUE T = TVAL.
C
C WE FIRST DETERMINE L SUCH THAT TVAL IS IN THE L-TH INTERVAL
C [TAU(L),TAU(L+1)]. NOTE: L = 0 MEANS THAT TVAL.LT.TAU(1).
C L = NDATA MEANS THAT TAU(NDATA).LE.TVAL .
C
DO 320 I=1,NDATA
IF (TVAL.LT.TAU(I)) THEN
L=I - 1
GO TO 325
END IF
320 CONTINUE
L=NDATA
325 CONTINUE
C
C CASE I: L = 0
C
IF (L.EQ.0) THEN
TEMP=2.0D0*(TAU(2) - TAU(1)) + (TAU(3) - TAU(2))
BPRIM1=( -1.5D0)/TEMP
BPRIM2= -BPRIM1
BETA(1)=BMAT(0,1) - BPRIM1*(TAU(1) - TVAL)
BETA(2)=BMAT(-1,2) - BPRIM2*(TAU(1) - TVAL)
DO 330 M=1,3
330 GVAL(M)=GHAT(M,1)*BETA(1) + GHAT(M,2)*BETA(2)
END IF
C
C CASE II: L = 1, 2, ... , NDATA-1 .
C
IF ((L.GE.1).AND.(L.LE.NDATA-1)) THEN
ISTART=MAX(1,L-1)
IEND=MIN(NDATA,L+2)
DO 340 I=ISTART,IEND
340 CALL CBETA(NDATA,ETAU,I,TVAL,L,BETA(I))
DO 350 M=1,3
GVAL(M)=0.0D0
DO 355 I=ISTART,IEND
355 GVAL(M)=GVAL(M) + GHAT(M,I)*BETA(I)
350 CONTINUE
END IF
C
C CASE III: L = NDATA
C
IF (L.EQ.NDATA) THEN
TEMP=2.0D0*(TAU(NDATA) - TAU(NDATA-1))
& + (TAU(NDATA-1) - TAU(NDATA-2))
BPRIM1=(1.5D0)/TEMP
BPRIM2= -BPRIM1
BETA(NDATA)=BMAT(0,NDATA) + BPRIM1*(TVAL - TAU(NDATA))
BETA(NDATA-1)=BMAT(1,NDATA-1) + BPRIM2*(TVAL - TAU(NDATA))
DO 360 M=1,3
360 GVAL(M)=GHAT(M,NDATA-1)*BETA(NDATA-1)
& + GHAT(M,NDATA)*BETA(NDATA)
END IF
C
RETURN
END
C
C *** QMULT ***
C
C SUBROUTINE TO CALCULATE THE PRODUCT OF TWO QUATERNIONS: C = A*B .
C
SUBROUTINE QMULT(A,B,C)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(0:3),B(0:3),C(0:3)
C
C(0)=A(0)*B(0) - A(1)*B(1) - A(2)*B(2) - A(3)*B(3)
C(1)=A(0)*B(1) + A(1)*B(0) + A(2)*B(3) - A(3)*B(2)
C(2)=A(0)*B(2) + A(2)*B(0) - A(1)*B(3) + A(3)*B(1)
C(3)=A(0)*B(3) + A(3)*B(0) + A(1)*B(2) - A(2)*B(1)
C
RETURN
END
C
C *** CALCOM ***
C
C THIS SUBROUTINE CALCULATES THE MATRIX OMEGA, BASED ON THE SEQUENCE
C OF TIMES: TAU(1), TAU(2), ... , TAU(N). NOTE: IN THIS SUBROUTINE,
C N IS USED AS SHORT FOR NDATA.
C
C OMEGA IS THE NXN MATRIX WHOSE ELEMENTS ARE:
C OMEGA(I,J) = THE INTEGRAL OF THE PRODUCT OF THE SECOND DERIVATIVES
C OF BETA(I,X) AND BETA(J,X), OVER THE INTERVAL (TAU(1), TAU(N)).
C HERE BETA(I,X) IS THE I-TH B-SPLINE ASSOCIATED WITH THE SEQUENCE
C TAU(1), TAU(2), ... , TAU(N).
C
C (REMARK: FOR I = 3, 4, ... , N-2, BETA(I,X) IS THE FUNCTION B(I-2,4,X)
C DEFINED ON P. 108 OF DE BOOR'S BOOK: A PRACTICAL GUIDE TO SPLINES.
C FOR I = 1, 2, BETA(I,X) IS DEFINED BY:
C BETA(1,X) = B(-2,4,X) + 0.5*B(-1,4,X) ,
C BETA(2,X) = B(0,4,X) + 0.5*B(-1,4,X) .
C FOR I = N-1, N:
C BETA(N-1,X) = B(N-3,4,X) + 0.5*B(N-2,4,X) ,
C BETA(N,X) = B(N-1,4,X) + 0.5*B(N-2,4,X) .
C THE SEQUENCE TAU(I) IS EXTENDED BELOW BY SYMMETRY ABOUT TAU(1), AND
C ABOVE BY SYMMETRY ABOUT TAU(N). THESE DEFINITIONS CAUSE BETA(1,X)
C AND BETA(2,X) TO HAVE ZERO SECOND DERIVATIVE AT X = TAU(1), AND
C CAUSE BETA(N-1,X) AND BETA(N,X) TO HAVE ZERO SECOND DERIVATIVE AT
C X = TAU(N).)
C
C INPUTS:
C N - NO. OF DATA POINTS (SHORT FOR NDATA). MUST BE AT LEAST 4.
C TAU - ARRAY CONTAINING THE TIMES TAU(I)
C
C OUTPUTS:
C OM - ARRAY CONATINING THE NON-ZERO DIAGONALS OF THE MATRIX OMEGA
C OM(0,I) = OMEGA(I,I)
C OM(1,I) = OMEGA(I,I+1)
C OM(2,I) = OMEGA(I,I+2)
C OM(3,I) = OMEGA(I,I+3)
C (NOTE: OMEGA IS SYMMETRIC.)
C
SUBROUTINE CALCOM(N,TAU,OM)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT)
DIMENSION OM(0:3,NDATAT)
DIMENSION AA(-1:1,NDATAT),QQ(0:1,NDATAT)
DIMENSION CC(-2:2,NDATAT),H(NDATAT-1)
C
C THE MATRIX OMEGA IS CALCULATED USING THE EQUATION:
C
C OMEGA = A*Q*(A**T) .
C
C A IS THE NXN MATRIX WHOSE ELEMENTS A(I,K) ARE THE COEFFICIENTS
C USED IN EXPRESSING THE SECOND DERIVATIVE OF BETA(I,X) AS A LINEAR
C COMBINATION OF THE FUNCTIONS B(K-1,2,X), DEFINED ON P. 108 OF
C DE BOOR'S BOOK. THAT IS, THE SECOND DERIVATIVE OF BETA(I,X) IS
C EQUAL TO THE SUM OVER K, K = 1, 2, ... , N, OF THE TERMS:
C
C A(I,K)*B(K-1,2,X) .
C
C (THE FORMULAS FOR THE A(I,K) ARE DERIVED FROM THE FORMULAS ON
C P. 138 IN DE BOOR'S BOOK.)
C
C THE MATRIX Q IS THE NXN MATRIX WITH THE ELEMENTS:
C Q(K,L) = THE INTEGRAL OF THE PRODUCT OF B(K-1,2,X) AND
C B(L-1,2,X) OVER THE INTERVAL (TAU(1), TAU(N)).
C
C NOTE:
C THE ARRAY AA CONTAINS THE NON-ZERO DIAGONALS OF THE MATRIX A:
C AA(-1,I) = A(I,I-1)
C AA(0,I) = A(I,I)
C AA(1,I) = A(I,I+1)
C THE ARRAY QQ CONTAINS THE NON-ZERO DIAGONALS OF THE MATRXI Q:
C QQ(0,I) = Q(I,I)
C QQ(1,I) = Q(I,I+1)
C (NOTE: THE MATRIX Q IS SYMMETRIC.)
C THE ARRAY CC CONTAINS THE NON-ZERO DIAGONALS OF THE MATRIX C = A*Q:
C CC(-2,I)=C(I,I-2)
C CC(-1,I)=C(I,I-1)
C CC(0,I)=C(I,I)
C CC(1,I)=C(I,I+1)
C CC(2,I)=C(I,I+2)
C
C CALCULATION OF THE MATRIX A.
C
IF (N.GE.5) THEN
DO 10 I=3,N-2
U=1.0D0/(TAU(I+1) - TAU(I-2))
V=1.0D0/(TAU(I+2) - TAU(I-1))
X=1.0D0/(TAU(I) - TAU(I-2))
Y=1.0D0/(TAU(I+1) - TAU(I-1))
Z=1.0D0/(TAU(I+2) - TAU(I))
AA(-1,I)=6.0D0*U*X
AA(0,I)= -6.0D0*(U + V)*Y
AA(1,I)=6.0D0*V*Z
10 CONTINUE
END IF
C
AA(0,1)=0.0D0
U=1.0D0/(TAU(2) + TAU(3) - 2.0D0*TAU(1))
X=1.0D0/(TAU(3) - TAU(1))
AA(1,1)=3.0D0*U*X
C
AA(-1,2)=0.0D0
V=1.0D0/(TAU(4) - TAU(1))
Y=1.0D0/(TAU(4) - TAU(2))
AA(0,2)= -3.0D0*X*(2.0D0*V + U)
AA(1,2)=6.0D0*V*Y
C
AA(0,N)=0.0D0
U=1.0D0/(2.0D0*TAU(N) - TAU(N-1) - TAU(N-2))
X=1.0D0/(TAU(N) - TAU(N-2))
AA(-1,N)=3.0D0*U*X
C
AA(1,N-1)=0.0D0
V=1.0D0/(TAU(N) - TAU(N-3))
Y=1.0D0/(TAU(N-1) - TAU(N-3))
AA(0,N-1)= -3.0D0*X*(U + 2.0D0*V)
AA(-1,N-1)=6.0D0*V*Y
C
C CALCULATION OF THE MATRIX Q.
C
DO 20 I=1,N-1
20 H(I)=TAU(I+1) - TAU(I)
C
DO 25 K=2,N-1
25 QQ(0,K)=(H(K-1) + H(K))/3.0D0
C
QQ(0,1)=H(1)/3.0D0
QQ(0,N)=H(N-1)/3.0D0
C
DO 30 K=1,N-1
30 QQ(1,K)=H(K)/6.0D0
C
C CALCULATION OF THE PRODUCT C = A*Q.
C
DO 40 I=2,N-1
40 CC(0,I)=AA(-1,I)*QQ(1,I-1) + AA(0,I)*QQ(0,I)
& + AA(1,I)*QQ(1,I)
C
CC(0,1)=AA(0,1)*QQ(0,1) + AA(1,1)*QQ(1,1)
CC(0,N)=AA(-1,N)*QQ(1,N-1) + AA(0,N)*QQ(0,N)
C
DO 45 I=1,N-1
45 CC(1,I)=AA(0,I)*QQ(1,I) + AA(1,I)*QQ(0,I+1)
C
DO 50 I=1,N-2
50 CC(2,I)=AA(1,I)*QQ(1,I+1)
C
DO 55 I=2,N
55 CC(-1,I)=AA(-1,I)*QQ(0,I-1) + AA(0,I)*QQ(1,I-1)
C
DO 60 I=3,N
60 CC(-2,I)=AA(-1,I)*QQ(1,I-2)
C
C CALCULATION OF OMEGA = C*(A**T).
C
DO 70 I=2,N-1
70 OM(0,I)=CC(-1,I)*AA(-1,I) + CC(0,I)*AA(0,I)
& + CC(1,I)*AA(1,I)
C
OM(0,1)=CC(0,1)*AA(0,1) + CC(1,1)*AA(1,1)
OM(0,N)=CC(-1,N)*AA(-1,N) + CC(0,N)*AA(0,N)
C
DO 75 I=1,N-2
75 OM(1,I)=CC(0,I)*AA(-1,I+1) + CC(1,I)*AA(0,I+1)
& + CC(2,I)*AA(1,I+1)
C
OM(1,N-1)=CC(0,N-1)*AA(-1,N) + CC(1,N-1)*AA(0,N)
C
DO 80 I=1,N-2
80 OM(2,I)=CC(1,I)*AA(-1,I+2) + CC(2,I)*AA(0,I+2)
C
DO 85 I=1,N-3
85 OM(3,I)=CC(2,I)*AA(-1,I+3)
C
RETURN
END
C
C *** CBETA ***
C
C THIS SUBROUTINE EVALUATES THE B-SPLINE BETA(I,X), I = 1, 2, ... , NDATA.
C
C INPUTS:
C N - NO. OF DATA POINTS (N IS SHORT FOR NDATA)
C TAU(-2:NDATAT+3)
C - ARRAY CONTAINING THE LOCATIONS OF THE ACTUAL, AS WELL AS
C FICTITIOUS, DATA POINTS
C I - INDEX OF THE B-SPLINE BETA(I,X)
C X - VALUE OF THE INDEPENDENT VARIABLE. (X MUST LIE IN THE
C CLOSED INTERVAL FROM TAU(1) TO TAU(N).)
C L - INDEX OF THE INTERVAL IN WHICH X LIES. (THAT IS, X LIES IN
C THE CLOSED INTERVAL FROM TAU(L) TO TAU(L+1). L CAN TAKE THE
C VALUES: L = 1, 2, ... N-1.)
C
C OUTPUTS:
C BETAV - VALUE OF BETA(I,X)
C
C NOTE:
C FOR I = 3, 4, ... , N-2, BETA(I,X) IS THE FUNCTION B(I-2,4,X) DEFINED
C ON P. 108 OF DE BOOR'S BOOK: A PRACTICAL GUIDE TO SPLINES.
C FOR I = 1, 2:
C BETA(1,X) = B(-2,4,X) + 0.5*B(-1,4,X)
C BETA(2,X) = B(0,4,X) + 0.5*B(-1,4,X)
C
C FOR I = N-1, N:
C BETA(N-1,X) = B(N-3,4,X) + 0.5*B(N-2,4,X)
C BETA(N,X) = B(N-1,4,X) + 0.5*B(N-2,4,X)
C
SUBROUTINE CBETA(N,TAU,I,X,L,BETAV)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(-2:NDATAT+3)
C
IF (I.EQ.1) THEN
CALL BSPL(N,TAU,-1,X,L,BVAL1)
CALL BSPL(N,TAU,-2,X,L,BVAL2)
BETAV=BVAL2 + 0.5D0*BVAL1
ELSE IF (I.EQ.2) THEN
CALL BSPL(N,TAU,-1,X,L,BVAL1)
CALL BSPL(N,TAU,0,X,L,BVAL2)
BETAV=BVAL2 + 0.5D0*BVAL1
ELSE IF (I.EQ.N-1) THEN
CALL BSPL(N,TAU,N-2,X,L,BVAL1)
CALL BSPL(N,TAU,N-3,X,L,BVAL2)
BETAV=BVAL2 + 0.5D0*BVAL1
ELSE IF (I.EQ.N) THEN
CALL BSPL(N,TAU,N-2,X,L,BVAL1)
CALL BSPL(N,TAU,N-1,X,L,BVAL2)
BETAV=BVAL2 + 0.5D0*BVAL1
ELSE
CALL BSPL(N,TAU,I-2,X,L,BETAV)
END IF
C
RETURN
END
C
C *** BSPL ***
C
C THIS SUBROUTINE EVALUATES THE B-SPLINE B(I,4,X), DEFINED ON P. 108
C OF DE BOOR'S BOOK: A PRACTICAL GUIDE TO SPLINES.
C
C INPUTS:
C N - NO. OF DATA POINTS (N IS SHORT FOR NDATA)
C TAU(-2:NDATAT+3)
C - ARRAY CONTAINING THE LOCATIONS OF THE ACTUAL, AS WELL AS
C FICTITIOUS, DATA POINTS
C I - INDEX OF THE B-SPLINE TO BE EVALUATED. (I CAN TAKE THE
C VALUES: I = -2, -1, 0, 1, ... , N-1.)
C X - VALUE OF THE INDEPENDENT VARIABLE AT WHICH TO EVALUATE
C THE B-SPLINE. (X MUST LIE IN THE CLOSED INTERVAL FROM
C TAU(1) TO TAU(N).)
C L - INDEX OF THE INTERVAL IN WHICH X LIES. (THAT IS, X LIES
C IN THE CLOSED INTERVAL FROM TAU(L) TO TAU(L+1). L CAN
C TAKE THE VALUES: L = 1, 2, ... , N-1.)
C
C OUTPUTS:
C BVALUE- VALUE OF B(I,4,X)
C
C NOTE: THE LOCATIONS OF THE FICTITIOUS DATA POINTS ARE:
C TAU(-2), TAU(-1), TAU(0)
C AND:
C TAU(N+1), TAU(N+2), TAU(N+3) .
C
SUBROUTINE BSPL(N,TAU,I,X,L,BVALUE)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(-2:NDATAT+3)
C
IF (L.EQ.I) THEN
U=1.0D0/(TAU(I+3) - TAU(I))
V=1.0D0/(TAU(I+2) - TAU(I))
W=1.0D0/(TAU(I+1) - TAU(I))
D=X - TAU(I)
BVALUE=U*V*W*D*D*D
ELSE IF (L.EQ.I+1) THEN
U=1.0D0/(TAU(I+3) - TAU(I))
V=1.0D0/(TAU(I+2) - TAU(I))
R=1.0D0/(TAU(I+4) - TAU(I+1))
S=1.0D0/(TAU(I+3) - TAU(I+1))
T=1.0D0/(TAU(I+2) - TAU(I+1))
D=X - TAU(I)
E=X - TAU(I+1)
F=TAU(I+2) - X
G=TAU(I+3) - X
H=TAU(I+4) - X
B03=V*T*D*F + S*T*E*G
B13=S*T*E*E
BVALUE=U*D*B03 + R*H*B13
ELSE IF (L.EQ.I+2) THEN
U=1.0D0/(TAU(I+3) - TAU(I))
R=1.0D0/(TAU(I+4) - TAU(I+1))
S=1.0D0/(TAU(I+3) - TAU(I+1))
P=1.0D0/(TAU(I+4) - TAU(I+2))
Q=1.0D0/(TAU(I+3) - TAU(I+2))
D=X - TAU(I)
E=X - TAU(I+1)
F1=X - TAU(I+2)
G=TAU(I+3) - X
H=TAU(I+4) - X
B03=S*Q*G*G
B13=S*Q*E*G + P*Q*H*F1
BVALUE=U*D*B03 + R*H*B13
ELSE IF (L.EQ.I+3) THEN
R=1.0D0/(TAU(I+4) - TAU(I+1))
P=1.0D0/(TAU(I+4) - TAU(I+2))
Z=1.0D0/(TAU(I+4) - TAU(I+3))
H=TAU(I+4) - X
BVALUE=R*P*Z*H*H*H
ELSE
BVALUE=0.0D0
END IF
C
RETURN
END
C
C *** ANORM ***
C
DOUBLE PRECISION FUNCTION ANORM(X,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION X(N)
C CALCULATES THE NORM OF A VECTOR X(NX1)
C
ANORM = 0.0D0
DO 1 I=1,N
ANORM=ANORM+X(I)*X(I)
1 CONTINUE
ANORM=DSQRT(ANORM)
RETURN
END
C
C *** CLONG ***
C
C CONVERTS LONGITUDE VALUES TO THE RANGE [SLONG,SLONG + 360).
C
DOUBLE PRECISION FUNCTION CLONG(PHI,SLONG)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
RATIO=(PHI - SLONG)/360.0D0
K=LARGIN(RATIO)
CLONG=PHI - DBLE(K)*360.0D0
C
RETURN
END
C
C *** LARGIN ***
C
C LARGIN(X) = THE LARGEST INTEGER WHICH IS LESS THAN OR EQUAL TO X.
C
INTEGER FUNCTION LARGIN(X)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
LARGIN=INT(X)
IF (X.LT.0.0D0) THEN
R=X - DBLE(LARGIN)
IF (R.LT.0.0D0) LARGIN=LARGIN - 1
END IF
C
RETURN
END
C
C *** CBMAT ***
C
C THIS SUBROUTINE CALCULATES THE ARRAY BMAT. BMAT IS A 3XNDATA ARRAY
C CONTAINING THE QUANTITIES:
C
C B(I,J) = BETA(J,TAU(I)) .
C
C IN PARTICULAR:
C
C BMAT(-1,J) CONTAINS B(J-1,J) ,
C BMAT(0,J) CONTAINS B(J,J) ,
C BMAT(1,J) CONTAINS B(J+1,J) .
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU - ARRAY CONTAINING THE TIMES TAU(I), I = 1, 2, ... , NDATA
C
C OUTPUTS:
C BMAT - THE 3XNDATA ARRAY BMAT
C ETAU - EXTENDED ARRAY OF TIMES
C
SUBROUTINE CBMAT(NDATA,TAU,BMAT,ETAU)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),BMAT(-1:1,NDATAT)
DIMENSION ETAU(-2:NDATAT+3)
C
C EXTEND THE ARRAY TAU TO AN EXTENDED ARRAY ETAU.
C
CALL CETAU(NDATA,TAU,ETAU)
C
C CALCULATE THE ARRAY BMAT.
C NOTE: I = INDEX OF THE ARGUMENT TAU(I)
C L = INDEX OF THE INTERVAL IN WHICH TAU(I) LIES
C J = INDEX OF THE B-SPLINE
C
DO 90 I=1,NDATA
L=MIN(NDATA-1,I)
JSTART=MAX(1,I-1)
JEND=MIN(NDATA,I+1)
DO 85 J=JSTART,JEND
CALL CBETA(NDATA,ETAU,J,TAU(I),L,BMAT(I-J,J))
85 CONTINUE
90 CONTINUE
C
RETURN
END
C
C *** CLCEIG ***
C
C SUBROUTINE TO CALCULATE THE EIGENVALUES OF THE MATRICES OMEGA
C AND
C
C G = TRANSPOSE(B)*SIGINV*B ,
C
C AND TO CALCULATE THE RATIO:
C
C RHO = (MAX. EIGENVALUE OF OMEGA)/(MAX. EIGENVALUE OF G) .
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - TIMES
C SIGINV(3,3,NDATAT) - ERROR MATRICES
C IPRINT - PRINTING PARAMETER
C 0 - NO PRINTING
C 231 - THE MATRIX OMEGA
C 232 - THE MATRIX B
C 233 - THE MATRIX G
C 234 - THE EIGENVALUES OF OMEGA
C 235 - THE EIGENVALUES OF G
C
C OUTPUTS:
C RHO - RATIO OF MAX EIGENVALUE OF OMEGA TO MAX. EIGENVALUE OF G
C
SUBROUTINE CLCEIG(NDATA,TAU,SIGINV,IPRINT,RHO)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),SIGINV(3,3,NDATAT)
DIMENSION OM(0:3,NDATAT)
DIMENSION BMAT(-1:1,NDATAT),ETAU(-2:NDATAT+3)
DIMENSION GMAT(3,3,0:2,NDATAT)
DIMENSION OMSMAT(NDATAT,4),OMEVAL(NDATAT)
DIMENSION EOM(NDATAT),EOM2(NDATAT),ZB(1,1)
DIMENSION GBSMAT(3*NDATAT,9),GEVAL(3*NDATAT)
DIMENSION EG(3*NDATAT),EG2(3*NDATAT)
LOGICAL MATZ
C
C CALCULATE THE MATRIX OMEGA.
C
CALL CALCOM(NDATA,TAU,OM)
IF (IPRINT.EQ.231) THEN
PRINT *,' '
PRINT *,'THE DIAGONALS OF THE MATRIX OMEGA'
PRINT 900
DO 400 I=1,NDATA-3
400 PRINT 901,I,OM(0,I),OM(1,I),OM(2,I),OM(3,I)
DO 410 I=NDATA-2,NDATA
410 PRINT 901,I,(OM(K,I),K=0,NDATA-I)
900 FORMAT(' ',2X,'I',4X,'OMEGA(I,I)',2X,'OMEGA(I,I+1)',
& 2X,'OMEGA(I,I+2)',2X,'OMEGA(I,I+3)')
901 FORMAT(' ',I3,4(2X,G12.5))
END IF
C
C CALCULATE THE MATRIX B AND THE ARRAY ETAU.
C
CALL CBMAT(NDATA,TAU,BMAT,ETAU)
IF (IPRINT.EQ.232) THEN
PRINT *,' '
PRINT *,'THE ARRAY BMAT:'
PRINT *,' BMAT(-1,I) CONTAINS B(I-1,I) = BETA(I,TAU(I-1))'
PRINT *,' BMAT(0,I) CONTAINS B(I,I) = BETA(I,TAU(I))'
PRINT *,' BMAT(1,I) CONTAINS B(I+1,I) = BETA(I,TAU(I+1))'
PRINT 922
PRINT 923,1,BMAT(0,1),BMAT(1,1)
DO 440 I=2,NDATA-1
440 PRINT 924,I,(BMAT(J,I),J= -1,1)
PRINT 925,NDATA,BMAT(-1,NDATA),BMAT(0,NDATA)
922 FORMAT(' ',2X,'I',3X,'BMAT(-1,I)',5X,'BMAT(0,I)',
& 5X,'BMAT(1,I)')
923 FORMAT(' ',I3,14X,2(2X,G12.5))
924 FORMAT(' ',I3,3(2X,G12.5))
925 FORMAT(' ',I3,2(2X,G12.5))
PRINT *,' '
PRINT *,'THE ARRAY ETAU:'
PRINT 930
DO 450 I=-2,NDATA+3
450 PRINT 931,I,ETAU(I)
930 FORMAT(' ',3X,'I',6X,'ETAU(I)')
931 FORMAT(' ',I4,2X,G12.5)
END IF
C
C FOR CONVENIENCE, WE SET N = 3*NDATA.
C
N=3*NDATA
C
C CALCULATE THE ARRAY GMAT WHICH CONTAINS THE MATRIX G:
C
C G = TRANSPOSE(B)*SIGINV*B .
C
C SUPPOSE WE CONSIDER G TO BE A BLOCK-MATRIX, WITH 3X3 SUB-BLOCKS
C G(I,J). WE DENOTE THE (K,L)-TH ELEMENT OF G(I,J) BY G(K,L,I,J).
C THEN:
C
C GMAT(K,L,0,I) CONTAINS G(K,L,I,I) ,
C GMAT(K,L,1,I) CONTAINS G(K,L,I,I+1) ,
C GMAT(K,L,2,I) CONTAINS G(K,L,I,I+2) .
C
C BRIEFLY STATED, G(I,J) IS THE SUM OVER M OF THE TERMS:
C
C B(M,I)*B(M,J)*SIGINV(M) .
C
C FIRST WE CALCULATE THE SUB-BLOCKS G(I,I).
C
DO 20 I=1,NDATA
MSTART=MAX(1,I-1)
MSTOP=MIN(NDATA,I+1)
DO 15 K=1,3
DO 15 L=1,3
GMAT(K,L,0,I)=0.0D0
DO 10 M=MSTART,MSTOP
TERM=BMAT(M-I,I)*BMAT(M-I,I)*SIGINV(K,L,M)
10 GMAT(K,L,0,I)=GMAT(K,L,0,I) + TERM
15 CONTINUE
20 CONTINUE
C
C NEXT WE CALCULATE THE SUB-BLOCKS G(I,I+1).
C
DO 40 I=1,NDATA-1
DO 35 K=1,3
DO 35 L=1,3
GMAT(K,L,1,I)=0.0D0
DO 30 M=I,I+1
TERM=BMAT(M-I,I)*BMAT(M-I-1,I+1)*SIGINV(K,L,M)
30 GMAT(K,L,1,I)=GMAT(K,L,1,I) + TERM
35 CONTINUE
40 CONTINUE
C
C FINALLY WE CALCULATE THE SUB-BLOCKS G(I,I+2).
C
DO 60 I=1,NDATA-2
DO 55 K=1,3
DO 55 L=1,3
GMAT(K,L,2,I)=BMAT(1,I)*BMAT(-1,I+2)*SIGINV(K,L,I+1)
55 CONTINUE
60 CONTINUE
C
IF (IPRINT.EQ.233) THEN
PRINT *,' '
PRINT *,'THE ARRAY GMAT: GMAT CONTAINS THE (3*NDATA)X(3*NDATA)'
PRINT *,'MATRIX G = TRANSPOSE(B)*SIGINV*B .'
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I)'
DO 460 I=1,NDATA
PRINT *,'I = ',I
DO 465 K=1,3
465 PRINT 940,(GMAT(K,L,0,I),L=1,3)
460 CONTINUE
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I+1)'
DO 470 I=1,NDATA-1
PRINT *,'I = ',I
DO 475 K=1,3
475 PRINT 940,(GMAT(K,L,1,I),L=1,3)
470 CONTINUE
PRINT *,' '
PRINT *,'THE SUB-BLOCKS G(I,I+2)'
DO 480 I=1,NDATA-2
PRINT *,'I = ',I
DO 485 K=1,3
485 PRINT 940,(GMAT(K,L,2,I),L=1,3)
480 CONTINUE
940 FORMAT(' ',3(2X,G12.5))
END IF
C
C SET UP THE ARRAY OMSMAT, WHICH CONTAINS THE MATRIX OMEGA-S IN
C BAND SYMMETRIC STORAGE MODE (EISPACK VERSION). THE 4-TH COL. OF
C OMSMAT CONTAINS THE MAIN DIAGONAL OF OMEGA-S. COLS. 3, 2, 1 OF
C OMSMAT CONTAIN THE 1-ST, 2-ND, 3-RD SUPERDIAGONALS OF OMEGA-S,
C RESPECTIVELY, FILLED IN BY ZEROS ON THE TOP.
C
DO 80 I=1,NDATA
80 OMSMAT(I,4)=OM(0,I)
C
DO 85 I=2,NDATA
85 OMSMAT(I,3)=OM(1,I-1)
OMSMAT(1,3)=0.0D0
C
DO 90 I=3,NDATA
90 OMSMAT(I,2)=OM(2,I-2)
OMSMAT(1,2)=0.0D0
OMSMAT(2,2)=0.0D0
C
DO 95 I=4,NDATA
95 OMSMAT(I,1)=OM(3,I-3)
OMSMAT(1,1)=0.0D0
OMSMAT(2,1)=0.0D0
OMSMAT(3,1)=0.0D0
C
C CALCULATE THE EIGENVALUES OF OMEGA-S: TO DO THIS, WE USE THE
C ROUTINES BANDR AND IMTQL1 FROM EISPACK.
C NM = ROW DIMENSION OF OMSMAT = NDATAT
C NEIG = ORDER OF OMEGA-S = NDATA
C MB = NO. OF SUPERDIAGONALS + 1 = 4
C MATZ = .FALSE. (DO NOT CALCULATE ORTH. TRANSF. MATRIX ZB)
C
NM=NDATAT
NEIG=NDATA
MB=4
MATZ=.FALSE.
C
C NOTE THAT THE ARRAY ZB(1,1) IS A DUMMY ARRAY AND IS NOT USED.
C
CALL BANDR(NM,NEIG,MB,OMSMAT,OMEVAL,EOM,EOM2,MATZ,ZB)
CALL IMTQL1(NEIG,OMEVAL,EOM,IERR)
IF (IERR.NE.0) THEN
PRINT *,' '
PRINT *,'ERROR IN FINDING THE EIGENVALUES OF OMEGA-S IN'
PRINT *,'CLCEIG, USING IMTQL1'
PRINT *,'IERR = ',IERR
STOP
END IF
C
C THE VECTOR OMEVAL CONTAINS THE EIGENVALUES OF OMEGA-S IN ASCENDING
C ORDER.
C
IF (IPRINT.EQ.234) THEN
PRINT *,' '
PRINT *,'THE EIGENVALUES OF OMEGA-S'
PRINT 950
DO 490 I=1,NDATA
490 PRINT 951,I,OMEVAL(I)
950 FORMAT(' ',2X,'I',3X,'OMEVAL(I)')
951 FORMAT(' ',I3,2X,G12.5)
END IF
C
C SET UP THE ARRAY GBSMAT, WHICH CONTAINS THE MATRIX G IN BAND
C SYMMETRIC STORAGE MODE (EISPACK VERSION). THE 9-TH COL. OF
C GBSMAT CONTAINS THE MAIN DIAGONAL OF G. COLS. 8, 7, ... , 1 OF
C GBSMAT CONTAIN THE 1-ST, 2-ND, ... , 8-TH SUPERDIAGONALS OF G,
C RESPECTIVELY, FILLED IN BY ZEROS ON THE TOP.
C
C FIRST INITIALIZE GBSMAT TO ZERO.
C
DO 105 MROW=1,N
DO 105 MCOL=1,9
105 GBSMAT(MROW,MCOL)=0.0D0
C
C NEXT WE FILL IN THE NON-ZERO ENTRIES IN GBSMAT. THIS IS DONE
C BY CONSIDERING THE APPROPRIATE 3X3 SUB-BLOCKS G(I,J) IN TURN,
C AND SETTING UP THAT PART OF GBSMAT WHICH CORRESPONDS TO G(I,J).
C THE APPROPRIATE SUB-BLOCKS G(I,J) ARE THOSE SUCH THAT J RANGES
C FROM I TO MIN(I+2,NDATA), WHILE I RANGES FROM 1 TO NDATA.
C
DO 140 I=1,NDATA
JMAX=MIN(I+2,NDATA)
DO 135 J=I,JMAX
C
C DENOTE THE (K,L)-TH ENTRY OF G(I,J) BY G(K,L,I,J). THEN THE
C ROW AND COLUMN INDICES IN THE BIG MATRIX G ARE:
C IROW = 3*(I-1) + K ,
C ICOL = 3*(J-1) + L .
C WE ARE INTERESTED ONLY IN ELEMENTS ON THE DIAGONAL OR ON ONE OF
C THE FIRST 8 SUPERDIAGONALS OF G. THEREFORE, WE REQUIRE ICOL TO
C BE BETWEEN IROW AND IROW+8.
C
DO 130 K=1,3
DO 130 L=1,3
IROW=3*(I-1) + K
ICOL=3*(J-1) + L
IF ((ICOL.LT.IROW).OR.(ICOL.GT.IROW+8)) GO TO 120
C
C THE ROW AND COLUMN INDICES OF THE POSITION IN THE ARRAY GBSMAT
C WHICH IS TO RECEIVE G(K,L,I,J) ARE:
C MROW = ICOL
C MCOL = 9 - (ICOL - IROW)
C NOTE: G(K,L,I,J) IS CONTAINED IN GMAT(K,L,J-I,I).
C
MROW=ICOL
MCOL=9 - (ICOL - IROW)
GBSMAT(MROW,MCOL)=GMAT(K,L,J-I,I)
120 CONTINUE
130 CONTINUE
C
C END OF LOOP ON K AND L.
C
135 CONTINUE
140 CONTINUE
C
C END OF LOOP ON I AND J.
C
C CALCULATE THE EIGENVALUES OF THE MATRIX G: AGAIN, WE USE THE
C ROUTINES BANDR AND IMTQL1 FROM EISPACK.
C NM = ROW DIMENSION OF GBSMAT = 3*NDATAT
C NEIG = ORDER OF G = 3*NDATA
C MB = NO. OF SUPERDIAGONALS + 1 = 9
C MATZ = .FALSE. (DO NOT CALCULATE ORTH. TRANSF. MATRIX ZB)
C
NM=3*NDATAT
NEIG=3*NDATA
MB=9
MATZ=.FALSE.
C
C AGAIN, THE ARRAY ZB(1,1) IS A DUMMY ARRAY AND IS NOT USED.
C
CALL BANDR(NM,NEIG,MB,GBSMAT,GEVAL,EG,EG2,MATZ,ZB)
CALL IMTQL1(NEIG,GEVAL,EG,IERR)
IF (IERR.NE.0) THEN
PRINT *,' '
PRINT *,'ERROR IN FINDING THE EIGENVALUES OF G IN'
PRINT *,'CLCEIG, USING IMTQL1'
PRINT *,'IERR = ',IERR
STOP
END IF
C
C THE VECTOR GEVAL CONTAINS THE EIGENVALUES OF G IN ASCENDING ORDER.
C
IF (IPRINT.EQ.235) THEN
PRINT *,' '
PRINT *,'THE EIGENVALUES OF G'
PRINT 960
DO 500 I=1,N
500 PRINT 961,I,GEVAL(I)
960 FORMAT(' ',2X,'I',4X,'GEVAL(I)')
961 FORMAT(' ',I3,2X,G12.5)
END IF
C
C CALCULATE RHO.
C
RHO=OMEVAL(NDATA)/GEVAL(N)
C
RETURN
END
C
C *** CALCR1 ***
C
C SUBROUTINE TO CALCULATE THE ROTATION R OF S**3 WHICH CORRESPONDS TO
C PARALLEL TRANSPORT ALONG THE GEODESIC JOINING Y TO X. HERE X AND Y
C ARE POINTS ON S**3, NOT ANTI-PODAL.
C
C NOTE THAT R MAPS Y TO X, AND THE TANGENT SPACE TO S**3 AT Y ONTO THE
C TANGENT SPACE TO S**3 AT X.
C
C INPUTS:
C Y(4) - STARTING POINT, ON S**3
C X(4) - ENDING POINT, ON S**3
C
C OUTPUTS:
C R(4,4) - ROTATION OF S**3 WHICH CORRESPONDS TO PARALLEL TRANSPORT
C ALONG THE GEODESIC JOINING Y TO X
C IERR - ERROR INDICATOR:
C 0 - ALL IS WELL
C 1 - ERROR, SINCE Y IS ANTI-PODAL TO X
C
C CALCR1 USES THE SUBPROGRAM ANORM.
C
SUBROUTINE CALCR1(Y,X,R,IERR)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION Y(4),X(4),A(4),U(4),R(4,4)
DIMENSION A1(4),A2(4),P(4)
C
DO 10 I=1,4
IF (Y(I).NE. -X(I)) GO TO 15
10 CONTINUE
C
C ERROR RETURN, SINCE Y = -X.
C
IERR=1
RETURN
C
15 CONTINUE
IERR=0
C
C CALCULATE D = THE ANGULAR DISTANCE (IN RADIANS ON S**3) BETWEEN Y
C AND X.
C
PROD=0.0D0
DO 30 J=1,4
30 PROD=PROD + Y(J)*X(J)
PROD=MIN(PROD,1.0D0)
D=DACOS(PROD)
C
C CALCULATE A = THE COMPONENT OF Y WHICH IS ORTHOGONAL TO X.
C
DO 40 J=1,4
40 A(J)=Y(J) - PROD*X(J)
C
C IF A = 0, THEN R = I. OTHERWISE, WE CALCULATE U = A/(NORM OF A).
C
AN=ANORM(A,4)
C
IF (AN.EQ.0.0D0) THEN
DO 42 I=1,4
DO 42 J=1,4
IF (I.EQ.J) THEN
R(I,J)=1.0D0
ELSE
R(I,J)=0.0D0
END IF
42 CONTINUE
RETURN
ELSE
DO 45 J=1,4
45 U(J)=A(J)/AN
END IF
C
C LET A1 = X AND A2 = U. EXTEND THE PAIR A1, A2 TO AN ORTHONORMAL BASIS
C FOR R**4: A1, A2, A3, A4. THE ROTATION R SATISFIES:
C
C R*A1 = COS(D)*A1 - SIN(D)*A2
C R*A2 = SIN(D)*A1 + COS(D)*A2
C R*A3 = A3
C R*A4 = A4 .
C
C WE CAN CALCULATE R AS FOLLOWS. LET E(1), E(2), E(3), E(4) BE THE
C STANDARD BASIS FOR R**4. WE KNOW:
C
C E(L) = A1(L)*A1 + A2(L)*A2 + P(L) ,
C
C WHERE P(L) IS THE ORTHOGONAL PROJECTION OF E(L) ONTO THE SUBSPACE OF
C R**4 SPANNED BY A3 AND A4. HENCE:
C
C R*E(L) = A1(L)*R*A1 + A2(L)*R*A2 + P(L) .
C
C THIS GIVES THE EXPRESSION:
C
C R*E(L) = (A1(L)*COS(D) + A2(L)*SIN(D))*A1
C
C + ( -A1(L)*SIN(D) + A2(L)*COS(D))*A2 + P(L) .
C
DO 48 J=1,4
A1(J)=X(J)
48 A2(J)=U(J)
C
CD=DCOS(D)
SD=DSIN(D)
C
DO 70 L=1,4
C
C CALCULATE P = P(L).
C
DO 50 J=1,4
IF (J.EQ.L) THEN
P(J)=1.0D0
ELSE
P(J)=0.0D0
END IF
50 CONTINUE
DO 55 J=1,4
55 P(J)=P(J) - A1(L)*A1(J) - A2(L)*A2(J)
C
C CALCULATE R*E(L).
C
C1=A1(L)*CD + A2(L)*SD
C2= -A1(L)*SD + A2(L)*CD
DO 60 J=1,4
60 R(J,L)=C1*A1(J) + C2*A2(J) + P(J)
C
70 CONTINUE
C
RETURN
END
C
C *** IVEXP2 ***
C
C SUBROUTINE TO CALCULATE THE INVERSE EXPONENTIAL OF A UNIT QUATERNION
C WHICH IS NEAR 1.
C
C INPUTS:
C U(0:3) - A UNIT QUATERNION NEAR 1
C
C OUTPUTS:
C X(3) - THE PURE QUATERNION SUCH THAT EXP(X) = U. (X(0) IS NOT
C STORED SINCE IT IS ZERO.)
C IERR - ERROR INDICATOR:
C 0 - ALL IS WELL
C 1 - ERROR, SINCE U = -1
C
SUBROUTINE IVEXP2(U,X,IERR)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION U(0:3),X(3)
C
TEMP=U(1)*U(1) + U(2)*U(2) + U(3)*U(3)
YNORM=DSQRT(TEMP)
IF (YNORM.EQ.0.0D0) THEN
IF (U(0).LT.0.0D0) THEN
IERR=1
RETURN
ELSE
IERR=0
DO 10 I=1,3
10 X(I)=0.0D0
RETURN
END IF
END IF
C
IERR=0
IF (U(0).GT.0.5D0) THEN
FACT=DASIN(YNORM)/YNORM
ELSE IF (U(0).LT.-0.5D0) THEN
PI=4.0D0*DATAN(1.0D0)
FACT=(PI - DASIN(YNORM))/YNORM
ELSE
FACT=DACOS(U(0))/YNORM
END IF
DO 20 I=1,3
20 X(I)=FACT*U(I)
C
RETURN
END
C
C *** CLCEXP ***
C
C SUBROUTINE TO CALCULATE THE EXPONENTIAL OF A PURE QUATERNION.
C
C INPUTS:
C X(3) - A PURE QUATERNION. (X(0) IS NOT STORED SINCE IT IS ZERO.)
C
C OUTPUTS:
C U(0:3) - THE UNIT QUATERNION = EXP(X).
C
C NOTE THAT THE MAPPING U = EXP(X) IS NOT ONE-TO-ONE UNLESS X IS
C RESTRICTED, FOR EX., BY THE CONDITION THAT IT HAS NORM LESS THAN PI.
C
SUBROUTINE CLCEXP(X,U)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION X(3),U(0:3)
C
XNORM=ANORM(X,3)
IF (XNORM.EQ.0.0D0) THEN
U(0)=1.0D0
DO 10 I=1,3
10 U(I)=0.0D0
ELSE
U(0)=COS(XNORM)
FACT=SIN(XNORM)/XNORM
DO 20 I=1,3
20 U(I)=FACT*X(I)
END IF
C
RETURN
END
C
C *** CLCRZ2 ***
C
C GIVEN A VECTOR A ON S**3, THIS SUBROUTINE FINDS A ROTATION RZ OF R**4
C WHICH MAPS A INTO E(1) = (1,0,0,0)**T.
C
C SPECIFICALLY, IF A IS NOT ANTI-PODAL TO E(1), RZ IS TAKEN TO BE THE
C ROTATION OF R**4 WHICH CORRESPONDS TO PARALLEL TRANSPORT ALONG THE
C GEODESIC ON S**3 JOINING A TO E(1). THE CASE WHERE A IS ANTI-PODAL
C TO E(1) IS TREATED SEPARATELY.
C
SUBROUTINE CLCRZ2(A,RZ)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(4),RZ(4,4),E1(4)
C
E1(1)=1.0D0
DO 10 J=2,4
10 E1(J)=0.0D0
C
CALL CALCR1(A,E1,RZ,IERR)
IF (IERR.EQ.1) THEN
DO 20 J=1,4
DO 20 L=1,4
20 RZ(J,L)=0.0D0
RZ(1,1)= -1.0D0
RZ(2,2)= -1.0D0
RZ(3,3)=1.0D0
RZ(4,4)=1.0D0
END IF
C
RETURN
END
C
C *** CONV ***
C
C SUBROUTINE TO CONVERT A POINT ON S**3 TO AXIS-LATITUDE, AXIS-LONGITUDE,
C ANGLE-OF-ROTATION FORM.
C
C INPUTS:
C QVECT(4) - POINT ON S**3
C SLONG - STARTING VALUE FOR REPORTING AXIS LONGITUDES (SO THAT
C THESE LONGITUDES LIE IN THE RANGE [SLONG,SLONG + 360))
C PI - PI (IN DOUBLE PRECISION)
C
C OUTPUTS:
C ALAT - AXIS LATITUDE (IN DEGREES)
C ALONG - AXIS LONGITUDE (IN DEGREES)
C RHO - ANGLE OF ROTATION (IN DEGREES)
C
C THE SUBROUTINE CONV CALLS THE ROUTINES ANORM AND CLONG.
C
SUBROUTINE CONV(QVECT,SLONG,PI,ALAT,ALONG,RHO)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION QVECT(4),UVECT(3)
C
FACT=180.0D0/PI
RHO=2.0D0*DACOS(QVECT(1))*FACT
DO 10 J=1,3
10 UVECT(J)=QVECT(J+1)
UVNORM=ANORM(UVECT(1),3)
IF (UVNORM.EQ.0.0D0) THEN
ALAT=0.0D0
ALONG=0.0D0
ELSE
DO 20 J=1,3
20 UVECT(J)=UVECT(J)/UVNORM
ALAT=DASIN(UVECT(3))*FACT
U12=DSQRT(UVECT(1)*UVECT(1) + UVECT(2)*UVECT(2))
IF (U12.EQ.0.0D0) THEN
ALONG=0.0D0
ELSE
ALONG=DATAN2(UVECT(2)/U12,UVECT(1)/U12)*FACT
END IF
END IF
ALONG=CLONG(ALONG,SLONG)
C
RETURN
END
C
C *** SORT ***
C
C SUBROUTINE TO SORT AN ARRAY OF NDATA NOS. INTO ASCENDING ORDER.
C
C INPUTS:
C NDATA - NO. OF VALUES IN ARRAY TO BE SORTED (MUST BE AT LEAST 2;
C AT MOST 50)
C VAL(NDATAT) - ARRAY TO BE SORTED
C
C OUTPUTS:
C SVAL(NDATAT) - SORTED ARRAY (IN ASCENDING ORDER)
C INDEX(NDATAT)- ARRAY OF INDICES, GIVING THE ORIGINAL INDEX VALUES
C FOR ITEMS IN THE ARRAY SVAL
C
SUBROUTINE SORT(NDATA,VAL,SVAL,INDEX)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION VAL(NDATAT),SVAL(NDATAT)
DIMENSION INDEX(NDATAT)
C
C INITIALIZE THE ARRAYS INDEX AND SVAL.
C
DO 10 I=1,NDATA
INDEX(I)=I
10 SVAL(I)=VAL(I)
C
C THE SORT OPERATES BY MAKING SUCCESSIVE PASSES THROUGH THE ARRAY SVAL,
C INTERCHANGING ANY ITEMS WHICH ARE IN REVERSE ORDER. THE SORT STOPS
C WHEN A PASS IS COMPLETED WITH NO INTERCHANGES.
C
20 ICNT=0
C
C NOTE: ICNT COUNTS THE NO. OF INTERCHANGES IN A GIVEN PASS.
C
DO 30 I=1,NDATA - 1
IF (SVAL(I).GT.SVAL(I+1)) THEN
TEMP=SVAL(I)
SVAL(I)=SVAL(I+1)
SVAL(I+1)=TEMP
ITEMP=INDEX(I)
INDEX(I)=INDEX(I+1)
INDEX(I+1)=ITEMP
ICNT=ICNT + 1
END IF
30 CONTINUE
C
IF (ICNT.GT.0) GO TO 20
C
RETURN
END
C
C *** CPSIF ***
C
C SUBROUTINE TO CALCULATE THE PSI FUNCTION. TO DEFINE THIS FUNCTION,
C LET GNU BE A NO. IN THE INTERVAL [0,1]. LET ALAMBDA BE THE VALUE
C CORRESP. TO GNU, ALAMBDA LYING IN THE INTERVAL [0,INFINITY]. IF
C G(ALAMBDA,T) IS THE SPLINE FIT TO THE DATA, USING THE PARAMETER
C ALAMBDA, DEFINE D(ALAMBDA,I) TO BE THE RESIDUAL:
C D(ALAMBDA,I) = WSTAR(I) - G(ALAMBDA,TAU(I)) .
C LET PHI(ALAMBDA,I) DENOTE THE SIZE OF THIS RESIDUAL RELATIVE TO 7.815:
C
C PHI(ALAMBDA,I) =
C (D(ALAMBDA,I)**T)*SIGINV(I)*D(ALAMBDA,I)/7.815 .
C
C PHIF(ALAMBDA) IS THE (NDATA-K)-TH LARGEST VALUE PHI(ALAMBDA,I), I.E.,
C THE SMALLEST VALUE PHI SUCH THAT NDATA-K VALUES PHI(ALAMBDA,I) ARE
C SMALLER THAN IT. PSIF(GNU) = PHIF(ALAMBDA).
C
C NOTE: ALAMBDA IS DETERMINED FROM GNU ACCORDING TO THE RULE:
C
C ALAMBDA = (GNU/(1.0 - GNU))*(1.0/RHO) FOR GNU LESS THAN 1.O ,
C ALAMBDA = INFINITY FOR GNU EQUAL TO 1.O .
C
C HERE RHO IS THE RATIO OF THE MAX. EIGENVALUE OF THE MATRIX OMEGA TO
C THE MAX. EIGENVALUE OF THE MATRIX G.
C
C REMARK: THE FITTED VALUE G(ALAMBDA,TAU(I)) LIES IN THE 95% CONFIDENCE
C REGION AT WSTAR(I) IF AND ONLY IF PHI(ALAMBDA,I) IS LESS OR EQUAL 1.0.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - DATA TIMES
C WSTAR(3,NDATAT) - UNROLLED DATA POINTS
C SIGINV(3,3,NDATAT) - INVERSES OF THE COVARIANCE MATRICES AT THE DATA
C POINTS
C RHO - RATIO OF MAX. EIGENVALUE OF OMEGA TO MAX. EIGENVALUE OF G
C GNU - SMOOTHING PARAMETER
C K - NO. OF PHI VALUES LARGER THAN PSIVAL. (K IS BETWEEN ZERO
C AND NDATA-1, INCLUSIVE.)
C IPRINT - PRINTING PARAMETER
C
C OUTPUTS:
C PSIVAL - VALUE OF PSIF(GNU)
C
SUBROUTINE CPSIF(NDATA,TAU,WSTAR,SIGINV,RHO,GNU,K,IPRINT,PSIVAL)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),WSTAR(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT)
DIMENSION GVAL(3),GHAT(3,NDATAT)
DIMENSION YFIT(3,NDATAT)
DIMENSION D(3,NDATAT),PHI(NDATAT)
DIMENSION SPHI(NDATAT)
DIMENSION INDEX(NDATAT)
C
C CALCULATE ALAMBDA:
C ILAMB = 1: ALAMBDA IS FINITE
C ILAMB = 2: ALAMBDA IS INFINITE
C
IF (ABS(1.0D0-GNU).LT.1.0D-12) THEN
ILAMB=2
ELSE
ILAMB=1
ALAMB=(GNU/(1.0D0-GNU))*(1.0D0/RHO)
END IF
C
C CALCULATE G(ALAMBDA,TAU(I)), I = 1, 2, ... , NDATA.
C
TVAL=(TAU(1) + TAU(2))/2.0D0
CALL SEFIT(NDATA,TAU,WSTAR,SIGINV,ILAMB,ALAMB,TVAL,IPRINT,GVAL,
&GHAT,YFIT,RCOND)
C
C CALCULATE D(ALAMBDA,I) AND PHI(ALAMBDA,I).
C
DO 20 I=1,NDATA
DO 22 J=1,3
22 D(J,I)=WSTAR(J,I) - YFIT(J,I)
SUM=0.0D0
DO 25 J=1,3
DO 25 L=1,3
25 SUM=SUM + D(J,I)*SIGINV(J,L,I)*D(L,I)
PHI(I)=SUM/7.815D0
20 CONTINUE
C
C SORT THE ARRAY PHI.
C
CALL SORT(NDATA,PHI,SPHI,INDEX)
C
PSIVAL=SPHI(NDATA-K)
C
RETURN
END
C
C *** FOLAM2 ***
C
C SUBROUTINE TO DETERMINE ALAMBDA BY REQUIRING THE FITTED SPLINE POINTS
C TO LIE IN THE KNOWN CONFIDENCE REGIONS.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - DATA TIMES
C WSTAR(3,NDATAT) - UNROLLED DATA POINTS
C SIGINV(3,3,NDATAT) - INVERSES OF THE COVARIANCE MATRICES AT THE DATA
C POINTS
C K - NO. OF DATA POINTS AT WHICH WE DO NOT REQUIRE THE FITTED
C POINT TO LIE IN THE CONFIDENCE REGION (MUST BE BETWEEN
C 0 AND NDATA-1, INCLUSIVE)
C IPRINT - PRINTING PARAMETER
C
C OUTPUTS:
C ILAMB - CODE FOR OPTIMAL VALUE OF ALAMBDA:
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ALAMB - VALUE OF ALAMBDA, IN CASE ILAMB = 1
C PHIVAL - VALUE OF THE PHI FUNCTION AT THE OPTIMAL VALUE OF ALAMBDA
C
C THE SUBROUTINE FOLAM2 CALLS THE ROUTINE CLCEIG.
C
C THE SUBROUTINE FOLAM2 FINDS THE LARGEST VALUE OF ALAMBDA SUCH THAT
C PHI(ALAMBDA) IS LESS OR EQUAL TO 1.0. NOTE:
C (A) IF PHI(ALAMBDA=INFINITY) IS LESS OR EQUAL TO 1.0, THEN THE OPTIMAL
C VALUE OF ALAMBDA IS ALAMBDA = INFINITY.
C (B) IF PHI(ALAMBDA=INFINITY) IS GREATER THAN 1.0, THEN THE OPTIMAL
C VALUE OF ALAMBDA IS GOTTEN BY SOLVING THE EQUATION:
C PHI(ALAMBDA) = 1.0 .
C
C FOR CONVENIENCE, THE ROUTINE FOLAM2 OPERATES NOT WITH THE FUNCTION
C PHI(ALAMBDA) BUT RATHER WITH THE FUNCTION PSI(GNU) = PHI(ALAMBDA)
C WHERE ALAMBDA IS DETERMINED FROM GNU ACCORDING TO THE RULE:
C ALAMBDA = (GNU/(1.0 - GNU))*(1.0/RHO) FOR GNU LESS THAN 1.0 ,
C ALAMBDA = INFINITY FOR GNU EQUAL TO 1.0 .
C THE VARIABLE GNU RANGES OVER THE INTERVAL [0,1]. HERE RHO IS THE RATIO
C OF THE MAX. EIGENVALUE OF THE MATRIX OMEGA TO THE MAX. EIGENVALUE OF
C THE MATRIX G.
C
C NOTE: WHEN SOLVING THE EQUATION:
C PSI(GNU = 1.0 ,
C IN THE CASE WHERE PSI(1.0) IS GREATER THAN 1.0, THE ROUTINE FOLAM2
C USES THE MODIFIED FALSE POSITION METHOD (SEE CONTE AND BE BOOR, ELEM.
C NUM. ANAL., 3RD ED., P. 77).
C
SUBROUTINE FOLAM2(NDATA,TAU,WSTAR,SIGINV,K,IPRINT,ILAMB,ALAMB,
&PHIVAL)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),WSTAR(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT)
DIMENSION A(0:50),B(0:50),W(0:50),F(0:50),G(0:50),P(0:50)
C
C CALCULATE RHO.
C
CALL CLCEIG(NDATA,TAU,SIGINV,IPRINT,RHO)
C
C CALCULATE THE VALUE OF THE PSI FUNCTION WHEN GNU = 1.0 .
C
CALL CPSIF(NDATA,TAU,WSTAR,SIGINV,RHO,1.0D0,K,IPRINT,PSIVAL)
IF (PSIVAL.LE.1.0D0) THEN
ILAMB=2
PHIVAL=PSIVAL
RETURN
ELSE
ILAMB=1
END IF
C
C AT THIS POINT WE KNOW THAT THE VALUE OF GNU SUCH THAT PSI(GNU) = 1.0
C IS LESS THAN 1.0 . WE DETERMINE THIS VALUE BY MEANS OF THE MODIFIED
C FALSE POSITION METHOD.
C
G(0)=PSIVAL - 1.0D0
CALL CPSIF(NDATA,TAU,WSTAR,SIGINV,RHO,0.0D0,K,IPRINT,PSIVAL)
F(0)=PSIVAL - 1.0D0
IF (F(0).GT.0.0D0) THEN
PRINT *,'ERROR IN THE ROUTINE FOLAM2. THE VALUE OF THE PSI'
PRINT *,'FUNCTION IS GREATER THAN 1.0 WHEN GNU = 0.0 .'
STOP
END IF
A(0)=0.0D0
B(0)=1.0D0
W(0)=A(0)
P(0)=PSIVAL
C
IF (IPRINT.EQ.270) THEN
PRINT *,' '
PRINT *,'ITERATION TO DETERMINE VALUE OF GNU SUCH THAT'
PRINT *,'PSI(GNU) = 1.0 :'
PRINT *,' '
PRINT 90,0
PRINT 91,A(0),B(0),W(0)
PRINT 92,F(0),G(0),P(0)
90 FORMAT(' ','ITERATION NO. N = ',I3)
91 FORMAT(' ',3X,'A = ',G12.5,' B = ',G12.5,' W = ',G12.5)
92 FORMAT(' ',3X,'F = ',G12.5,' G = ',G12.5,' P = ',G12.5)
END IF
C
DO 50 N=0,49
DEL=G(N) - F(N)
IF (DEL.LT.0.0D0) THEN
PRINT *,'ERROR IN THE ROUTINE FOLAM2. THE QUANTITY DEL ='
PRINT *,'G(N) - F(N) IS NEGATIVE.'
PRINT *,'N = ',N
PRINT *,'DEL = ',DEL
STOP
ELSE IF (DEL.EQ.0.0D0) THEN
PRINT *,'THE QUANTITY DEL = G(N) - F(N) IS ZERO.'
PRINT *,'N = ',N
PRINT *,'WE TAKE GNU TO BE THE AVERAGE OF A(N) AND B(N).'
ASTAR=(A(N) + B(N))/2.0D0
BSTAR=ASTAR
PSTAR=P(N)
NSTAR=N
GO TO 60
END IF
W(N+1)=(G(N)*A(N) - F(N)*B(N))/DEL
C
CALL CPSIF(NDATA,TAU,WSTAR,SIGINV,RHO,W(N+1),K,IPRINT,P(N+1))
C
IF (P(N+1).GE.1.0D0) THEN
A(N+1)=A(N)
B(N+1)=W(N+1)
G(N+1)=P(N+1) - 1.0D0
IF (P(N).GE.1.0D0) THEN
F(N+1)=F(N)/2.0D0
ELSE
F(N+1)=F(N)
END IF
ELSE
A(N+1)=W(N+1)
B(N+1)=B(N)
F(N+1)=P(N+1) - 1.0D0
IF (P(N).LT.1.0D0) THEN
G(N+1)=G(N)/2.0D0
ELSE
G(N+1)=G(N)
END IF
END IF
C
IF (IPRINT.EQ.270) THEN
PRINT *,' '
PRINT 90,N+1
PRINT 91,A(N+1),B(N+1),W(N+1)
PRINT 92,F(N+1),G(N+1),P(N+1)
END IF
C
EPS=0.05D0
IF ((B(N+1) - A(N+1)).LT.(EPS*A(N+1))) THEN
ASTAR=A(N+1)
BSTAR=B(N+1)
PSTAR=P(N+1)
NSTAR=N+1
GO TO 60
END IF
50 CONTINUE
C
PRINT *,' '
PRINT *,'NO CONVERGENCE IN FOLAM2 AFTER 50 ITERATIONS.'
STOP
C
60 CONTINUE
IF (IPRINT.EQ.270) THEN
PRINT *,' '
PRINT 93,NSTAR
PRINT 94,ASTAR,BSTAR,PSTAR
93 FORMAT(' ','CONVERGENCE AT ITERATION NO. N = ',I3)
94 FORMAT(' ',3X,'ASTAR = ',G12.5,' BSTAR = ',G12.5,
&' PSTAR = ',G12.5)
END IF
GNU=ASTAR
ALAMB=(GNU/(1.0D0 - GNU))*(1.0D0/RHO)
PHIVAL=PSTAR
C
RETURN
END
C
C *** CETAU ***
C
C SUBROUTINE TO CALCULATE EXTENDED ARRAY OF TIMES.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - ARRAY CONTAINING THE TIMES TAU(I), I = 1, ... , NDATA
C
C OUTPUTS:
C ETAU(-2:NDATAT+3) - EXTENDED ARRAY OF TIMES
C
SUBROUTINE CETAU(NDATA,TAU,ETAU)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),ETAU(-2:NDATAT+3)
C
DO 10 I=1,NDATA
10 ETAU(I)=TAU(I)
C
DO 15 I= -2,0
15 ETAU(I)=2.0D0*TAU(1) - TAU(2-I)
C
DO 20 I=NDATA+1,NDATA+3
20 ETAU(I)=2.0D0*TAU(NDATA) - TAU(2*NDATA-I)
C
RETURN
END
C
C *** CP2F ***
C
C SUBROUTINE TO CALCULATE THE PSI2 FUNCTION. TO DEFINE THIS FUNCTION,
C LET GNU BE A NO. IN THE INTERVAL [0,1]. LET ALAMBDA BE THE VALUE
C CORRESP. TO GNU, ALAMBDA LYING IN THE INTERVAL [0,INFINITY]. IF
C G(ALAMBDA,T) IS THE SPLINE FIT TO THE DATA, USING THE PARAMETER
C ALAMBDA, DEFINE D(ALAMBDA,I) TO BE THE RESIDUAL:
C D(ALAMBDA,I) = WSTAR(I) - G(ALAMBDA,TAU(I)) .
C LET PHI(ALAMBDA,I) DENOTE THE SIZE OF THIS RESIDUAL:
C
C PHI(ALAMBDA,I) =
C (D(ALAMBDA,I)**T)*SIGINV(I)*D(ALAMBDA,I) .
C
C PHI2F(ALAMBDA) IS DEFINED BY:
C
C PHI2F(ALAMBDA) =
C (SUM OF PHI(ALAMBDA,I) OVER I = 1, 2, ... , NDATA)/(3*NDATA) .
C
C PSI2F(GNU) = PHI2F(ALAMBDA).
C
C NOTE: ALAMBDA IS DETERMINED FROM GNU ACCORDING TO THE RULE:
C
C ALAMBDA = (GNU/(1.0 - GNU))*(1.0/RHO) FOR GNU LESS THAN 1.O ,
C ALAMBDA = INFINITY FOR GNU EQUAL TO 1.O .
C
C HERE RHO IS THE RATIO OF THE MAX. EIGENVALUE OF THE MATRIX OMEGA TO
C THE MAX. EIGENVALUE OF THE MATRIX G.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - DATA TIMES
C WSTAR(3,NDATAT) - UNROLLED DATA POINTS
C SIGINV(3,3,NDATAT) - INVERSES OF THE COVARIANCE MATRICES AT THE DATA
C POINTS
C RHO - RATIO OF MAX. EIGENVALUE OF OMEGA TO MAX. EIGENVALUE OF G
C GNU - SMOOTHING PARAMETER
C IPRINT - PRINTING PARAMETER
C
C OUTPUTS:
C PSI2V - VALUE OF PSI2F(GNU)
C
SUBROUTINE CP2F(NDATA,TAU,WSTAR,SIGINV,RHO,GNU,IPRINT,PSI2V)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),WSTAR(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT)
DIMENSION GVAL(3),GHAT(3,NDATAT)
DIMENSION YFIT(3,NDATAT)
DIMENSION D(3,NDATAT),PHI(NDATAT)
C
C CALCULATE ALAMBDA:
C ILAMB = 1: ALAMBDA IS FINITE
C ILAMB = 2: ALAMBDA IS INFINITE
C
IF (ABS(1.0D0-GNU).LT.1.0D-12) THEN
ILAMB=2
ELSE
ILAMB=1
ALAMB=(GNU/(1.0D0-GNU))*(1.0D0/RHO)
END IF
C
C CALCULATE G(ALAMBDA,TAU(I)), I = 1, 2, ... , NDATA.
C
TVAL=(TAU(1) + TAU(2))/2.0D0
CALL SEFIT(NDATA,TAU,WSTAR,SIGINV,ILAMB,ALAMB,TVAL,IPRINT,GVAL,
&GHAT,YFIT,RCOND)
C
C CALCULATE D(ALAMBDA,I) AND PHI(ALAMBDA,I).
C
DO 20 I=1,NDATA
DO 22 J=1,3
22 D(J,I)=WSTAR(J,I) - YFIT(J,I)
SUM=0.0D0
DO 25 J=1,3
DO 25 L=1,3
25 SUM=SUM + D(J,I)*SIGINV(J,L,I)*D(L,I)
PHI(I)=SUM
20 CONTINUE
C
C CALCULATE PSI2V.
C
SUM=0.0D0
DO 30 I=1,NDATA
30 SUM=SUM + PHI(I)
PSI2V=SUM/DBLE(3*NDATA)
C
RETURN
END
C
C *** FOLAM3 ***
C
C SUBROUTINE TO DETERMINE ALAMBDA BY THE DISCREPANCY METHOD.
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - DATA TIMES
C WSTAR(3,NDATAT) - UNROLLED DATA POINTS
C SIGINV(3,3,NDATAT) - INVERSES OF THE COVARIANCE MATRICES AT THE DATA
C POINTS
C IPRINT - PRINTING PARAMETER
C
C OUTPUTS:
C ILAMB - CODE FOR OPTIMAL VALUE OF ALAMBDA:
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ALAMB - VALUE OF ALAMBDA, IN CASE ILAMB = 1
C PHI2V - VALUE OF THE PHI2 FUNCTION AT THE OPTIMAL VALUE OF ALAMBDA
C
C THE SUBROUTINE FOLAM3 CALLS THE ROUTINE CLCEIG.
C
C THE SUBROUTINE FOLAM3 FINDS THE LARGEST VALUE OF ALAMBDA SUCH THAT
C PHI2(ALAMBDA) IS LESS OR EQUAL TO 1.0. NOTE:
C (A) IF PHI2(ALAMBDA=INFINITY) IS LESS OR EQUAL TO 1.0, THEN THE
C OPTIMAL VALUE OF ALAMBDA IS ALAMBDA = INFINITY.
C (B) IF PHI2(ALAMBDA=INFINITY) IS GREATER THAN 1.0, THEN THE OPTIMAL
C VALUE OF ALAMBDA IS GOTTEN BY SOLVING THE EQUATION:
C PHI2(ALAMBDA) = 1.0 .
C
C FOR CONVENIENCE, THE ROUTINE FOLAM3 OPERATES NOT WITH THE FUNCTION
C PHI2(ALAMBDA) BUT RATHER WITH THE FUNCTION PSI2(GNU) = PHI2(ALAMBDA)
C WHERE ALAMBDA IS DETERMINED FROM GNU ACCORDING TO THE RULE:
C ALAMBDA = (GNU/(1.0 - GNU))*(1.0/RHO) FOR GNU LESS THAN 1.0 ,
C ALAMBDA = INFINITY FOR GNU EQUAL TO 1.0 .
C THE VARIABLE GNU RANGES OVER THE INTERVAL [0,1]. HERE RHO IS THE RATIO
C OF THE MAX. EIGENVALUE OF THE MATRIX OMEGA TO THE MAX. EIGENVALUE OF
C THE MATRIX G.
C
C NOTE: WHEN SOLVING THE EQUATION:
C PSI2(GNU) = 1.0 ,
C IN THE CASE WHERE PSI2(1.0) IS GREATER THAN 1.0, THE ROUTINE FOLAM3
C USES THE MODIFIED FALSE POSITION METHOD (SEE CONTE AND BE BOOR, ELEM.
C NUM. ANAL., 3RD ED., P. 77).
C
SUBROUTINE FOLAM3(NDATA,TAU,WSTAR,SIGINV,IPRINT,ILAMB,ALAMB,
&PHI2V)
PARAMETER (NDATAT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),WSTAR(3,NDATAT)
DIMENSION SIGINV(3,3,NDATAT)
DIMENSION A(0:50),B(0:50),W(0:50),F(0:50),G(0:50),P(0:50)
C
C CALCULATE RHO.
C
CALL CLCEIG(NDATA,TAU,SIGINV,IPRINT,RHO)
C
C CALCULATE THE VALUE OF THE PSI2 FUNCTION WHEN GNU = 1.0 .
C
CALL CP2F(NDATA,TAU,WSTAR,SIGINV,RHO,1.0D0,IPRINT,PSI2V)
IF (PSI2V.LE.1.0D0) THEN
ILAMB=2
PHI2V=PSI2V
RETURN
ELSE
ILAMB=1
END IF
C
C AT THIS POINT WE KNOW THAT THE VALUE OF GNU SUCH THAT PSI2(GNU) = 1.0
C IS LESS THAN 1.0 . WE DETERMINE THIS VALUE BY MEANS OF THE MODIFIED
C FALSE POSITION METHOD.
C
G(0)=PSI2V - 1.0D0
CALL CP2F(NDATA,TAU,WSTAR,SIGINV,RHO,0.0D0,IPRINT,PSI2V)
F(0)=PSI2V - 1.0D0
IF (F(0).GT.0.0D0) THEN
PRINT *,'ERROR IN THE ROUTINE FOLAM3. THE VALUE OF THE PSI2'
PRINT *,'FUNCTION IS GREATER THAN 1.0 WHEN GNU = 0.0 .'
STOP
END IF
A(0)=0.0D0
B(0)=1.0D0
W(0)=A(0)
P(0)=PSI2V
C
IF (IPRINT.EQ.275) THEN
PRINT *,' '
PRINT *,'ITERATION TO DETERMINE VALUE OF GNU SUCH THAT'
PRINT *,'PSI2(GNU) = 1.0 :'
PRINT *,' '
PRINT 90,0
PRINT 91,A(0),B(0),W(0)
PRINT 92,F(0),G(0),P(0)
90 FORMAT(' ','ITERATION NO. N = ',I3)
91 FORMAT(' ',3X,'A = ',G12.5,' B = ',G12.5,' W = ',G12.5)
92 FORMAT(' ',3X,'F = ',G12.5,' G = ',G12.5,' P = ',G12.5)
END IF
C
DO 50 N=0,49
DEL=G(N) - F(N)
IF (DEL.LT.0.0D0) THEN
PRINT *,'ERROR IN THE ROUTINE FOLAM3. THE QUANTITY DEL ='
PRINT *,'G(N) - F(N) IS NEGATIVE.'
PRINT *,'N = ',N
PRINT *,'DEL = ',DEL
STOP
ELSE IF (DEL.EQ.0.0D0) THEN
PRINT *,'THE QUANTITY DEL = G(N) - F(N) IS ZERO.'
PRINT *,'N = ',N
PRINT *,'WE TAKE GNU TO BE THE AVERAGE OF A(N) AND B(N).'
ASTAR=(A(N) + B(N))/2.0D0
BSTAR=ASTAR
PSTAR=P(N)
NSTAR=N
GO TO 60
END IF
W(N+1)=(G(N)*A(N) - F(N)*B(N))/DEL
C
CALL CP2F(NDATA,TAU,WSTAR,SIGINV,RHO,W(N+1),IPRINT,P(N+1))
C
IF (P(N+1).GE.1.0D0) THEN
A(N+1)=A(N)
B(N+1)=W(N+1)
G(N+1)=P(N+1) - 1.0D0
IF (P(N).GE.1.0D0) THEN
F(N+1)=F(N)/2.0D0
ELSE
F(N+1)=F(N)
END IF
ELSE
A(N+1)=W(N+1)
B(N+1)=B(N)
F(N+1)=P(N+1) - 1.0D0
IF (P(N).LT.1.0D0) THEN
G(N+1)=G(N)/2.0D0
ELSE
G(N+1)=G(N)
END IF
END IF
C
IF (IPRINT.EQ.275) THEN
PRINT *,' '
PRINT 90,N+1
PRINT 91,A(N+1),B(N+1),W(N+1)
PRINT 92,F(N+1),G(N+1),P(N+1)
END IF
C
EPS=0.05D0
IF ((B(N+1) - A(N+1)).LT.(EPS*A(N+1))) THEN
ASTAR=A(N+1)
BSTAR=B(N+1)
PSTAR=P(N+1)
NSTAR=N+1
GO TO 60
END IF
50 CONTINUE
C
PRINT *,' '
PRINT *,'NO CONVERGENCE IN FOLAM3 AFTER 50 ITERATIONS.'
STOP
C
60 CONTINUE
IF (IPRINT.EQ.275) THEN
PRINT *,' '
PRINT 93,NSTAR
PRINT 94,ASTAR,BSTAR,PSTAR
93 FORMAT(' ','CONVERGENCE AT ITERATION NO. N = ',I3)
94 FORMAT(' ',3X,'ASTAR = ',G12.5,' BSTAR = ',G12.5,
&' PSTAR = ',G12.5)
END IF
GNU=ASTAR
ALAMB=(GNU/(1.0D0 - GNU))*(1.0D0/RHO)
PHI2V=PSTAR
C
RETURN
END
C
C We include several subroutines which we use from the package LINPACK,
C namely: dpbco, dpbfa, dpbsl, ddot, daxpy, dscal, dasum. Reference:
C LINPACK User's Guide, by J. J. Dongarra et al., SIAM, 1979.
C
cat > dpbco.f <<'CUT HERE............'
subroutine dpbco(abd,lda,n,m,rcond,z,info)
integer lda,n,m,info
double precision abd(lda,1),z(1)
double precision rcond
c
c dpbco factors a double precision symmetric positive definite
c matrix stored in band form and estimates the condition of the
c matrix.
c
c if rcond is not needed, dpbfa is slightly faster.
c to solve a*x = b , follow dpbco by dpbsl.
c to compute inverse(a)*c , follow dpbco by dpbsl.
c to compute determinant(a) , follow dpbco by dpbdi.
c
c on entry
c
c abd double precision(lda, n)
c the matrix to be factored. the columns of the upper
c triangle are stored in the columns of abd and the
c diagonals of the upper triangle are stored in the
c rows of abd . see the comments below for details.
c
c lda integer
c the leading dimension of the array abd .
c lda must be .ge. m + 1 .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c 0 .le. m .lt. n .
c
c on return
c
c abd an upper triangular matrix r , stored in band
c form, so that a = trans(r)*r .
c if info .ne. 0 , the factorization is not complete.
c
c rcond double precision
c an estimate of the reciprocal condition of a .
c for the system a*x = b , relative perturbations
c in a and b of size epsilon may cause
c relative perturbations in x of size epsilon/rcond .
c if rcond is so small that the logical expression
c 1.0 + rcond .eq. 1.0
c is true, then a may be singular to working
c precision. in particular, rcond is zero if
c exact singularity is detected or the estimate
c underflows. if info .ne. 0 , rcond is unchanged.
c
c z double precision(n)
c a work vector whose contents are usually unimportant.
c if a is singular to working precision, then z is
c an approximate null vector in the sense that
c norm(a*z) = rcond*norm(a)*norm(z) .
c if info .ne. 0 , z is unchanged.
c
c info integer
c = 0 for normal return.
c = k signals an error condition. the leading minor
c of order k is not positive definite.
c
c band storage
c
c if a is a symmetric positive definite band matrix,
c the following program segment will set up the input.
c
c m = (band width above diagonal)
c do 20 j = 1, n
c i1 = max0(1, j-m)
c do 10 i = i1, j
c k = i-j+m+1
c abd(k,j) = a(i,j)
c 10 continue
c 20 continue
c
c this uses m + 1 rows of a , except for the m by m
c upper left triangle, which is ignored.
c
c example.. if the original matrix is
c
c 11 12 13 0 0 0
c 12 22 23 24 0 0
c 13 23 33 34 35 0
c 0 24 34 44 45 46
c 0 0 35 45 55 56
c 0 0 0 46 56 66
c
c then n = 6 , m = 2 and abd should contain
c
c * * 13 24 35 46
c * 12 23 34 45 56
c 11 22 33 44 55 66
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c linpack dpbfa
c blas daxpy,ddot,dscal,dasum
c fortran dabs,dmax1,max0,min0,dreal,dsign
c
c internal variables
c
double precision ddot,ek,t,wk,wkm
double precision anorm,s,dasum,sm,ynorm
integer i,j,j2,k,kb,kp1,l,la,lb,lm,mu
c
c
c find norm of a
c
do 30 j = 1, n
l = min0(j,m+1)
mu = max0(m+2-j,1)
z(j) = dasum(l,abd(mu,j),1)
k = j - l
if (m .lt. mu) go to 20
do 10 i = mu, m
k = k + 1
z(k) = z(k) + dabs(abd(i,j))
10 continue
20 continue
30 continue
anorm = 0.0d0
do 40 j = 1, n
anorm = dmax1(anorm,z(j))
40 continue
c
c factor
c
call dpbfa(abd,lda,n,m,info)
if (info .ne. 0) go to 180
c
c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
c estimate = norm(z)/norm(y) where a*z = y and a*y = e .
c the components of e are chosen to cause maximum local
c growth in the elements of w where trans(r)*w = e .
c the vectors are frequently rescaled to avoid overflow.
c
c solve trans(r)*w = e
c
ek = 1.0d0
do 50 j = 1, n
z(j) = 0.0d0
50 continue
do 110 k = 1, n
if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
if (dabs(ek-z(k)) .le. abd(m+1,k)) go to 60
s = abd(m+1,k)/dabs(ek-z(k))
call dscal(n,s,z,1)
ek = s*ek
60 continue
wk = ek - z(k)
wkm = -ek - z(k)
s = dabs(wk)
sm = dabs(wkm)
wk = wk/abd(m+1,k)
wkm = wkm/abd(m+1,k)
kp1 = k + 1
j2 = min0(k+m,n)
i = m + 1
if (kp1 .gt. j2) go to 100
do 70 j = kp1, j2
i = i - 1
sm = sm + dabs(z(j)+wkm*abd(i,j))
z(j) = z(j) + wk*abd(i,j)
s = s + dabs(z(j))
70 continue
if (s .ge. sm) go to 90
t = wkm - wk
wk = wkm
i = m + 1
do 80 j = kp1, j2
i = i - 1
z(j) = z(j) + t*abd(i,j)
80 continue
90 continue
100 continue
z(k) = wk
110 continue
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
c
c solve r*y = w
c
do 130 kb = 1, n
k = n + 1 - kb
if (dabs(z(k)) .le. abd(m+1,k)) go to 120
s = abd(m+1,k)/dabs(z(k))
call dscal(n,s,z,1)
120 continue
z(k) = z(k)/abd(m+1,k)
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
t = -z(k)
call daxpy(lm,t,abd(la,k),1,z(lb),1)
130 continue
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
c
ynorm = 1.0d0
c
c solve trans(r)*v = y
c
do 150 k = 1, n
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
z(k) = z(k) - ddot(lm,abd(la,k),1,z(lb),1)
if (dabs(z(k)) .le. abd(m+1,k)) go to 140
s = abd(m+1,k)/dabs(z(k))
call dscal(n,s,z,1)
ynorm = s*ynorm
140 continue
z(k) = z(k)/abd(m+1,k)
150 continue
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
ynorm = s*ynorm
c
c solve r*z = w
c
do 170 kb = 1, n
k = n + 1 - kb
if (dabs(z(k)) .le. abd(m+1,k)) go to 160
s = abd(m+1,k)/dabs(z(k))
call dscal(n,s,z,1)
ynorm = s*ynorm
160 continue
z(k) = z(k)/abd(m+1,k)
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
t = -z(k)
call daxpy(lm,t,abd(la,k),1,z(lb),1)
170 continue
c make znorm = 1.0
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
ynorm = s*ynorm
c
if (anorm .ne. 0.0d0) rcond = ynorm/anorm
if (anorm .eq. 0.0d0) rcond = 0.0d0
180 continue
return
end
CUT HERE............
cat > dpbfa.f <<'CUT HERE............'
subroutine dpbfa(abd,lda,n,m,info)
integer lda,n,m,info
double precision abd(lda,1)
c
c dpbfa factors a double precision symmetric positive definite
c matrix stored in band form.
c
c dpbfa is usually called by dpbco, but it can be called
c directly with a saving in time if rcond is not needed.
c
c on entry
c
c abd double precision(lda, n)
c the matrix to be factored. the columns of the upper
c triangle are stored in the columns of abd and the
c diagonals of the upper triangle are stored in the
c rows of abd . see the comments below for details.
c
c lda integer
c the leading dimension of the array abd .
c lda must be .ge. m + 1 .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c 0 .le. m .lt. n .
c
c on return
c
c abd an upper triangular matrix r , stored in band
c form, so that a = trans(r)*r .
c
c info integer
c = 0 for normal return.
c = k if the leading minor of order k is not
c positive definite.
c
c band storage
c
c if a is a symmetric positive definite band matrix,
c the following program segment will set up the input.
c
c m = (band width above diagonal)
c do 20 j = 1, n
c i1 = max0(1, j-m)
c do 10 i = i1, j
c k = i-j+m+1
c abd(k,j) = a(i,j)
c 10 continue
c 20 continue
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas ddot
c fortran max0,dsqrt
c
c internal variables
c
double precision ddot,t
double precision s
integer ik,j,jk,k,mu
c begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0d0
ik = m + 1
jk = max0(j-m,1)
mu = max0(m+2-j,1)
if (m .lt. mu) go to 20
do 10 k = mu, m
t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1)
t = t/abd(m+1,jk)
abd(k,j) = t
s = s + t*t
ik = ik - 1
jk = jk + 1
10 continue
20 continue
s = abd(m+1,j) - s
c ......exit
if (s .le. 0.0d0) go to 40
abd(m+1,j) = dsqrt(s)
30 continue
info = 0
40 continue
return
end
CUT HERE............
cat > ddot.f <<'CUT HERE............'
double precision function ddot(n,dx,incx,dy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end
CUT HERE............
cat > daxpy.f <<'CUT HERE............'
subroutine daxpy(n,da,dx,incx,dy,incy)
c
c constant times a vector plus a vector.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),da
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (da .eq. 0.0d0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i + 1) = dy(i + 1) + da*dx(i + 1)
dy(i + 2) = dy(i + 2) + da*dx(i + 2)
dy(i + 3) = dy(i + 3) + da*dx(i + 3)
50 continue
return
end
CUT HERE............
cat > dscal.f <<'CUT HERE............'
subroutine dscal(n,da,dx,incx)
c
c scales a vector by a constant.
c uses unrolled loops for increment equal to one.
c jack dongarra, linpack, 3/11/78.
c modified to correct problem with negative increment, 8/21/90.
c
double precision da,dx(1)
integer i,incx,ix,m,mp1,n
c
if(n.le.0)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
do 10 i = 1,n
dx(ix) = da*dx(ix)
ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dx(i) = da*dx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 i = mp1,n,5
dx(i) = da*dx(i)
dx(i + 1) = da*dx(i + 1)
dx(i + 2) = da*dx(i + 2)
dx(i + 3) = da*dx(i + 3)
dx(i + 4) = da*dx(i + 4)
50 continue
return
end
CUT HERE............
cat > dasum.f <<'CUT HERE............'
double precision function dasum(n,dx,incx)
c
c takes the sum of the absolute values.
c uses unrolled loops for increment equal to one.
c jack dongarra, linpack, 3/11/78.
c modified to correct problem with negative increment, 8/21/90.
c
double precision dx(1),dtemp
integer i,incx,ix,m,mp1,n
c
dasum = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
do 10 i = 1,n
dtemp = dtemp + dabs(dx(ix))
ix = ix + incx
10 continue
dasum = dtemp
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,6)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dabs(dx(i))
30 continue
if( n .lt. 6 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,6
dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
* + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
50 continue
60 dasum = dtemp
return
end
CUT HERE............
cat > dpbsl.f <<'CUT HERE............'
subroutine dpbsl(abd,lda,n,m,b)
integer lda,n,m
double precision abd(lda,1),b(1)
c
c dpbsl solves the double precision symmetric positive definite
c band system a*x = b
c using the factors computed by dpbco or dpbfa.
c
c on entry
c
c abd double precision(lda, n)
c the output from dpbco or dpbfa.
c
c lda integer
c the leading dimension of the array abd .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c
c b double precision(n)
c the right hand side vector.
c
c on return
c
c b the solution vector x .
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal. technically this indicates
c singularity but it is usually caused by improper subroutine
c arguments. it will not occur if the subroutines are called
c correctly and info .eq. 0 .
c
c to compute inverse(a) * c where c is a matrix
c with p columns
c call dpbco(abd,lda,n,rcond,z,info)
c if (rcond is too small .or. info .ne. 0) go to ...
c do 10 j = 1, p
c call dpbsl(abd,lda,n,c(1,j))
c 10 continue
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,ddot
c fortran min0
c
c internal variables
c
double precision ddot,t
integer k,kb,la,lb,lm
c
c solve trans(r)*y = b
c
do 10 k = 1, n
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
t = ddot(lm,abd(la,k),1,b(lb),1)
b(k) = (b(k) - t)/abd(m+1,k)
10 continue
c
c solve r*x = y
c
do 20 kb = 1, n
k = n + 1 - kb
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
b(k) = b(k)/abd(m+1,k)
t = -b(k)
call daxpy(lm,t,abd(la,k),1,b(lb),1)
20 continue
return
end
CUT HERE............
C
C We also include some subroutines which we use from the package EISPACK,
C namely: bandr, imtql1, pythag, imtql2. References:
C (1) Matrix Eigensystem Routines - EISPACK Guide, by B. T. Smith et al.,
C Second Edition, Springer-Verlag, 1976
C (2) Matrix Eigensystem Routines - EISPACK Guide Extension, by B. S.
C Garbow et al., Springer-Verlag, 1977
C
cat > bandr.f <<'CUT HERE............'
subroutine bandr(nm,n,mb,a,d,e,e2,matz,z)
C REFORMULATED S2 IN LOOP 500 TO AVOID OVERFLOW. (9/29/89 BSG)
c
integer j,k,l,n,r,i1,i2,j1,j2,kr,mb,mr,m1,nm,n2,r1,ugl,maxl,maxr
double precision a(nm,mb),d(n),e(n),e2(n),z(nm,n)
double precision g,u,b1,b2,c2,f1,f2,s2,dmin,dminrt
logical matz
c
c this subroutine is a translation of the algol procedure bandrd,
c num. math. 12, 231-241(1968) by schwarz.
c handbook for auto. comp., vol.ii-linear algebra, 273-283(1971).
c
c this subroutine reduces a real symmetric band matrix
c to a symmetric tridiagonal matrix using and optionally
c accumulating orthogonal similarity transformations.
c
c on input
c
c nm must be set to the row dimension of two-dimensional
c array parameters as declared in the calling program
c dimension statement.
c
c n is the order of the matrix.
c
c mb is the (half) band width of the matrix, defined as the
c number of adjacent diagonals, including the principal
c diagonal, required to specify the non-zero portion of the
c lower triangle of the matrix.
c
c a contains the lower triangle of the symmetric band input
c matrix stored as an n by mb array. its lowest subdiagonal
c is stored in the last n+1-mb positions of the first column,
c its next subdiagonal in the last n+2-mb positions of the
c second column, further subdiagonals similarly, and finally
c its principal diagonal in the n positions of the last column.
c contents of storages not part of the matrix are arbitrary.
c
c matz should be set to .true. if the transformation matrix is
c to be accumulated, and to .false. otherwise.
c
c on output
c
c a has been destroyed, except for its last two columns which
c contain a copy of the tridiagonal matrix.
c
c d contains the diagonal elements of the tridiagonal matrix.
c
c e contains the subdiagonal elements of the tridiagonal
c matrix in its last n-1 positions. e(1) is set to zero.
c
c e2 contains the squares of the corresponding elements of e.
c e2 may coincide with e if the squares are not needed.
c
c z contains the orthogonal transformation matrix produced in
c the reduction if matz has been set to .true. otherwise, z
c is not referenced.
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
c this version dated september 1989.
c
c ------------------------------------------------------------------
c
dmin = 2.0d0**(-64)
dminrt = 2.0d0**(-32)
c .......... initialize diagonal scaling matrix ..........
do 30 j = 1, n
30 d(j) = 1.0d0
c
if (.not. matz) go to 60
c
do 50 j = 1, n
c
do 40 k = 1, n
40 z(j,k) = 0.0d0
c
z(j,j) = 1.0d0
50 continue
c
60 m1 = mb - 1
if (m1 - 1) 900, 800, 70
70 n2 = n - 2
c
do 700 k = 1, n2
maxr = min0(m1,n-k)
c .......... for r=maxr step -1 until 2 do -- ..........
do 600 r1 = 2, maxr
r = maxr + 2 - r1
kr = k + r
mr = mb - r
g = a(kr,mr)
a(kr-1,1) = a(kr-1,mr+1)
ugl = k
c
do 500 j = kr, n, m1
j1 = j - 1
j2 = j1 - 1
if (g .eq. 0.0d0) go to 600
b1 = a(j1,1) / g
b2 = b1 * d(j1) / d(j)
IF (ABS(B1) .GT. 1.0D0) THEN
U = 1.0D0 / B1
S2 = U / (U + B2)
ELSE
S2 = 1.0D0 / (1.0D0 + B1 * B2)
ENDIF
c
if (s2 .ge. 0.5d0 ) go to 450
b1 = g / a(j1,1)
b2 = b1 * d(j) / d(j1)
c2 = 1.0d0 - s2
d(j1) = c2 * d(j1)
d(j) = c2 * d(j)
f1 = 2.0d0 * a(j,m1)
f2 = b1 * a(j1,mb)
a(j,m1) = -b2 * (b1 * a(j,m1) - a(j,mb)) - f2 + a(j,m1)
a(j1,mb) = b2 * (b2 * a(j,mb) + f1) + a(j1,mb)
a(j,mb) = b1 * (f2 - f1) + a(j,mb)
c
do 200 l = ugl, j2
i2 = mb - j + l
u = a(j1,i2+1) + b2 * a(j,i2)
a(j,i2) = -b1 * a(j1,i2+1) + a(j,i2)
a(j1,i2+1) = u
200 continue
c
ugl = j
a(j1,1) = a(j1,1) + b2 * g
if (j .eq. n) go to 350
maxl = min0(m1,n-j1)
c
do 300 l = 2, maxl
i1 = j1 + l
i2 = mb - l
u = a(i1,i2) + b2 * a(i1,i2+1)
a(i1,i2+1) = -b1 * a(i1,i2) + a(i1,i2+1)
a(i1,i2) = u
300 continue
c
i1 = j + m1
if (i1 .gt. n) go to 350
g = b2 * a(i1,1)
350 if (.not. matz) go to 500
c
do 400 l = 1, n
u = z(l,j1) + b2 * z(l,j)
z(l,j) = -b1 * z(l,j1) + z(l,j)
z(l,j1) = u
400 continue
c
go to 500
c
450 u = d(j1)
d(j1) = s2 * d(j)
d(j) = s2 * u
f1 = 2.0d0 * a(j,m1)
f2 = b1 * a(j,mb)
u = b1 * (f2 - f1) + a(j1,mb)
a(j,m1) = b2 * (b1 * a(j,m1) - a(j1,mb)) + f2 - a(j,m1)
a(j1,mb) = b2 * (b2 * a(j1,mb) + f1) + a(j,mb)
a(j,mb) = u
c
do 460 l = ugl, j2
i2 = mb - j + l
u = b2 * a(j1,i2+1) + a(j,i2)
a(j,i2) = -a(j1,i2+1) + b1 * a(j,i2)
a(j1,i2+1) = u
460 continue
c
ugl = j
a(j1,1) = b2 * a(j1,1) + g
if (j .eq. n) go to 480
maxl = min0(m1,n-j1)
c
do 470 l = 2, maxl
i1 = j1 + l
i2 = mb - l
u = b2 * a(i1,i2) + a(i1,i2+1)
a(i1,i2+1) = -a(i1,i2) + b1 * a(i1,i2+1)
a(i1,i2) = u
470 continue
c
i1 = j + m1
if (i1 .gt. n) go to 480
g = a(i1,1)
a(i1,1) = b1 * a(i1,1)
480 if (.not. matz) go to 500
c
do 490 l = 1, n
u = b2 * z(l,j1) + z(l,j)
z(l,j) = -z(l,j1) + b1 * z(l,j)
z(l,j1) = u
490 continue
c
500 continue
c
600 continue
c
if (mod(k,64) .ne. 0) go to 700
c .......... rescale to avoid underflow or overflow ..........
do 650 j = k, n
if (d(j) .ge. dmin) go to 650
maxl = max0(1,mb+1-j)
c
do 610 l = maxl, m1
610 a(j,l) = dminrt * a(j,l)
c
if (j .eq. n) go to 630
maxl = min0(m1,n-j)
c
do 620 l = 1, maxl
i1 = j + l
i2 = mb - l
a(i1,i2) = dminrt * a(i1,i2)
620 continue
c
630 if (.not. matz) go to 645
c
do 640 l = 1, n
640 z(l,j) = dminrt * z(l,j)
c
645 a(j,mb) = dmin * a(j,mb)
d(j) = d(j) / dmin
650 continue
c
700 continue
c .......... form square root of scaling matrix ..........
800 do 810 j = 2, n
810 e(j) = dsqrt(d(j))
c
if (.not. matz) go to 840
c
do 830 j = 1, n
c
do 820 k = 2, n
820 z(j,k) = e(k) * z(j,k)
c
830 continue
c
840 u = 1.0d0
c
do 850 j = 2, n
a(j,m1) = u * e(j) * a(j,m1)
u = e(j)
e2(j) = a(j,m1) ** 2
a(j,mb) = d(j) * a(j,mb)
d(j) = a(j,mb)
e(j) = a(j,m1)
850 continue
c
d(1) = a(1,mb)
e(1) = 0.0d0
e2(1) = 0.0d0
go to 1001
c
900 do 950 j = 1, n
d(j) = a(j,mb)
e(j) = 0.0d0
e2(j) = 0.0d0
950 continue
c
1001 return
end
CUT HERE............
cat > imtql1.f <<'CUT HERE............'
subroutine imtql1(n,d,e,ierr)
c
integer i,j,l,m,n,ii,mml,ierr
double precision d(n),e(n)
double precision b,c,f,g,p,r,s,tst1,tst2,pythag
c
c this subroutine is a translation of the algol procedure imtql1,
c num. math. 12, 377-383(1968) by martin and wilkinson,
c as modified in num. math. 15, 450(1970) by dubrulle.
c handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
c
c this subroutine finds the eigenvalues of a symmetric
c tridiagonal matrix by the implicit ql method.
c
c on input
c
c n is the order of the matrix.
c
c d contains the diagonal elements of the input matrix.
c
c e contains the subdiagonal elements of the input matrix
c in its last n-1 positions. e(1) is arbitrary.
c
c on output
c
c d contains the eigenvalues in ascending order. if an
c error exit is made, the eigenvalues are correct and
c ordered for indices 1,2,...ierr-1, but may not be
c the smallest eigenvalues.
c
c e has been destroyed.
c
c ierr is set to
c zero for normal return,
c j if the j-th eigenvalue has not been
c determined after 30 iterations.
c
c calls pythag for dsqrt(a*a + b*b) .
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
c this version dated august 1983.
c
c ------------------------------------------------------------------
c
ierr = 0
if (n .eq. 1) go to 1001
c
do 100 i = 2, n
100 e(i-1) = e(i)
c
e(n) = 0.0d0
c
do 290 l = 1, n
j = 0
c .......... look for small sub-diagonal element ..........
105 do 110 m = l, n
if (m .eq. n) go to 120
tst1 = dabs(d(m)) + dabs(d(m+1))
tst2 = tst1 + dabs(e(m))
if (tst2 .eq. tst1) go to 120
110 continue
c
120 p = d(l)
if (m .eq. l) go to 215
if (j .eq. 30) go to 1000
j = j + 1
c .......... form shift ..........
g = (d(l+1) - p) / (2.0d0 * e(l))
r = pythag(g,1.0d0)
g = d(m) - p + e(l) / (g + dsign(r,g))
s = 1.0d0
c = 1.0d0
p = 0.0d0
mml = m - l
c .......... for i=m-1 step -1 until l do -- ..........
do 200 ii = 1, mml
i = m - ii
f = s * e(i)
b = c * e(i)
r = pythag(f,g)
e(i+1) = r
if (r .eq. 0.0d0) go to 210
s = f / r
c = g / r
g = d(i+1) - p
r = (d(i) - g) * s + 2.0d0 * c * b
p = s * r
d(i+1) = g + p
g = c * r - b
200 continue
c
d(l) = d(l) - p
e(l) = g
e(m) = 0.0d0
go to 105
c .......... recover from underflow ..........
210 d(i+1) = d(i+1) - p
e(m) = 0.0d0
go to 105
c .......... order eigenvalues ..........
215 if (l .eq. 1) go to 250
c .......... for i=l step -1 until 2 do -- ..........
do 230 ii = 2, l
i = l + 2 - ii
if (p .ge. d(i-1)) go to 270
d(i) = d(i-1)
230 continue
c
250 i = 1
270 d(i) = p
290 continue
c
go to 1001
c .......... set error -- no convergence to an
c eigenvalue after 30 iterations ..........
1000 ierr = l
1001 return
end
CUT HERE............
cat > pythag.f <<'CUT HERE............'
double precision function pythag(a,b)
double precision a,b
c
c finds dsqrt(a**2+b**2) without overflow or destructive underflow
c
double precision p,r,s,t,u
p = dmax1(dabs(a),dabs(b))
if (p .eq. 0.0d0) go to 20
r = (dmin1(dabs(a),dabs(b))/p)**2
10 continue
t = 4.0d0 + r
if (t .eq. 4.0d0) go to 20
s = r/t
u = 1.0d0 + 2.0d0*s
p = u*p
r = (s/u)**2 * r
go to 10
20 pythag = p
return
end
CUT HERE............
cat > imtql2.f <<'CUT HERE............'
subroutine imtql2(nm,n,d,e,z,ierr)
c
integer i,j,k,l,m,n,ii,nm,mml,ierr
double precision d(n),e(n),z(nm,n)
double precision b,c,f,g,p,r,s,tst1,tst2,pythag
c
c this subroutine is a translation of the algol procedure imtql2,
c num. math. 12, 377-383(1968) by martin and wilkinson,
c as modified in num. math. 15, 450(1970) by dubrulle.
c handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
c
c this subroutine finds the eigenvalues and eigenvectors
c of a symmetric tridiagonal matrix by the implicit ql method.
c the eigenvectors of a full symmetric matrix can also
c be found if tred2 has been used to reduce this
c full matrix to tridiagonal form.
c
c on input
c
c nm must be set to the row dimension of two-dimensional
c array parameters as declared in the calling program
c dimension statement.
c
c n is the order of the matrix.
c
c d contains the diagonal elements of the input matrix.
c
c e contains the subdiagonal elements of the input matrix
c in its last n-1 positions. e(1) is arbitrary.
c
c z contains the transformation matrix produced in the
c reduction by tred2, if performed. if the eigenvectors
c of the tridiagonal matrix are desired, z must contain
c the identity matrix.
c
c on output
c
c d contains the eigenvalues in ascending order. if an
c error exit is made, the eigenvalues are correct but
c unordered for indices 1,2,...,ierr-1.
c
c e has been destroyed.
c
c z contains orthonormal eigenvectors of the symmetric
c tridiagonal (or full) matrix. if an error exit is made,
c z contains the eigenvectors associated with the stored
c eigenvalues.
c
c ierr is set to
c zero for normal return,
c j if the j-th eigenvalue has not been
c determined after 30 iterations.
c
c calls pythag for dsqrt(a*a + b*b) .
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
c this version dated august 1983.
c
c ------------------------------------------------------------------
c
ierr = 0
if (n .eq. 1) go to 1001
c
do 100 i = 2, n
100 e(i-1) = e(i)
c
e(n) = 0.0d0
c
do 240 l = 1, n
j = 0
c .......... look for small sub-diagonal element ..........
105 do 110 m = l, n
if (m .eq. n) go to 120
tst1 = dabs(d(m)) + dabs(d(m+1))
tst2 = tst1 + dabs(e(m))
if (tst2 .eq. tst1) go to 120
110 continue
c
120 p = d(l)
if (m .eq. l) go to 240
if (j .eq. 30) go to 1000
j = j + 1
c .......... form shift ..........
g = (d(l+1) - p) / (2.0d0 * e(l))
r = pythag(g,1.0d0)
g = d(m) - p + e(l) / (g + dsign(r,g))
s = 1.0d0
c = 1.0d0
p = 0.0d0
mml = m - l
c .......... for i=m-1 step -1 until l do -- ..........
do 200 ii = 1, mml
i = m - ii
f = s * e(i)
b = c * e(i)
r = pythag(f,g)
e(i+1) = r
if (r .eq. 0.0d0) go to 210
s = f / r
c = g / r
g = d(i+1) - p
r = (d(i) - g) * s + 2.0d0 * c * b
p = s * r
d(i+1) = g + p
g = c * r - b
c .......... form vector ..........
do 180 k = 1, n
f = z(k,i+1)
z(k,i+1) = s * z(k,i) + c * f
z(k,i) = c * z(k,i) - s * f
180 continue
c
200 continue
c
d(l) = d(l) - p
e(l) = g
e(m) = 0.0d0
go to 105
c .......... recover from underflow ..........
210 d(i+1) = d(i+1) - p
e(m) = 0.0d0
go to 105
240 continue
c .......... order eigenvalues and eigenvectors ..........
do 300 ii = 2, n
i = ii - 1
k = i
p = d(i)
c
do 260 j = ii, n
if (d(j) .ge. p) go to 260
k = j
p = d(j)
260 continue
c
if (k .eq. i) go to 300
d(k) = d(i)
d(i) = p
c
do 280 j = 1, n
p = z(j,i)
z(j,i) = z(j,k)
z(j,k) = p
280 continue
c
300 continue
c
go to 1001
c .......... set error -- no convergence to an
c eigenvalue after 30 iterations ..........
1000 ierr = l
1001 return
end
CUT HERE............