C.hr utility.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     utility.f
C
C     Copyright (C) 1990.
C
C     A. Ronald Gallant & George E. Tauchen 
C     P.O. Box 5513, Raleigh NC 27650-5513, USA   
C
C     Permission to use, copy, modify, and distribute this
C     software and its documentation for any purpose and
C     without fee is hereby granted, provided that the
C     above copyright notice appear in all copies and that
C     both that copyright notice and this permission notice
C     appear in supporting documentation.
C
C     This software is provided "as is" without any
C     expressed or implied warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*
C.hr BLOCK DATA
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     IF DGMPNT OR HEADER ARE COPIED FROM THIS SET OF ROUTINES,
C     BE SURE TO INCLUDE THIS BLOCK DATA PROGRAM IN THE COPY.
C     REMOVE THE COMMENTS IF THESE UTILITIES ARE USED WITH A ROUTINE
C     OTHER THAN NLSYS OR GMM.
      BLOCK DATA
      COMMON /ZLNSIZ/ LNSIZE
      DATA LNSIZE/79/
      END
C.hr DGMPNT
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMPNT   7/16/87
C
C     PURPOSE
C     PRINT A MATRIX.
C
C     USAGE
C     CALL DGMPNT(A,M,N)
C
C     ARGUMENTS
C     A - AN M BY N MATRIX STORED COLUMNWISE (STORAGE MODE 0).
C         REAL*8
C     M - NUMBER OF ROWS IN A
C         INTEGER*4
C     N - NUMBER OF COLUMNS IN A
C         INTEGER*4
C
C     COMMENT
C     THE DEFAULT LINESIZE IS 133, 132 PLUS CARRIAGE CONTROL CHARACTER.
C     THE USAGE:
C     COMMON /ZLNSIZ/ LNSIZE
C     LNSIZE=80
C     WILL CHANGE THE LINESIZE TO 80.  LINESIZES BETWEEN 72 AND 133 ARE
C     PERMITTED.
C
C     PROGRAMMER
C     DR. A. RONALD GALLANT
C     DEPARTMENT OF STATISTICS
C     NORTH CAROLINA STATE UNIVERSITY
C     RALEIGH, NORTH CAROLINA  27695-8203
C
C
      SUBROUTINE DGMPNT(A,M,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 START,STOP,OUT
      REAL*8 A(M,N),F
      CHARACTER*1 DIGIT(10),DUMMY(8),COL(9)
      CHARACTER*8 TYPE(14),FMT(13)
      CHARACTER*104 CFMT
      COMMON /ZLNSIZ/ LNSIZE
      EQUIVALENCE (FMT(1),DUMMY(1)),(FMT(1),CFMT)
      SAVE
      DATA FMT   /'(   X,  ',11*'        ','       )'/
      DATA DIGIT /'0','1','2','3','4','5','6','7','8','9'/
      DATA COL   /6*' ','C','O','L'/
      DATA TYPE  /',0PF12.1',',1PD12.4',',0PF12.8',',0PF12.7',
     &            ',0PF12.6',',0PF12.5',',0PF12.4',',0PF12.3',
     &            ',0PF12.2',',0PF12.1',',0PF12.0',',9A1,I3 ',
     &            ' 6X     ','''ROW'',I3'/
      OUT=3
      LNSIZ=LNSIZE
      IF((LNSIZE.LT.72).OR.(LNSIZE.GT.133)) LNSIZ=133
      MAXCOL=(LNSIZ-8)/12
      IF(N.LT.MAXCOL) MAXCOL=N
      IPAD=(LNSIZ-8-12*MAXCOL)/2+1
      IPAD10=IPAD/10
      IPAD1=IPAD-10*IPAD10
      DUMMY(3)=DIGIT(IPAD10+1)
      DUMMY(4)=DIGIT(IPAD1+1)
      START=1
11    STOP=START-1+MAXCOL
      IF(STOP.GT.N) STOP=N
      K=2
      DO 13 J=START,STOP
      K=K+1
13    FMT(K)=TYPE(12)
      FMT(2)=TYPE(13)
      WRITE(OUT,3001)
      WRITE(OUT,3001)
      WRITE(OUT,CFMT) (COL,J,J=START,STOP)
      WRITE(OUT,3001)
      FMT(2)=TYPE(14)
      DO 19 I=1,M
      K=2
      DO 18 J=START,STOP
      K=K+1
      FMT(K)=TYPE(2)
      F=DABS(A(I,J))
      IF(F.LT.1.D+8 ) FMT(K)=TYPE(11)
      IF(F.LT.1.D+5 ) FMT(K)=TYPE(10)
      IF(F.LT.1.D+4 ) FMT(K)=TYPE( 9)
      IF(F.LT.1.D+3 ) FMT(K)=TYPE( 8)
      IF(F.LT.1.D+2 ) FMT(K)=TYPE( 7)
      IF(F.LT.1.D+1 ) FMT(K)=TYPE( 6)
C     IF(F.LT.1.D+0 ) FMT(K)=TYPE( 6)
      IF(F.LT.1.D-1 ) FMT(K)=TYPE( 5)
      IF(F.LT.1.D-2 ) FMT(K)=TYPE( 3)
      IF(F.LT.1.D-4 ) FMT(K)=TYPE( 2)
      IF(F.LT.1.D-55) FMT(K)=TYPE( 1)
18    CONTINUE
19    WRITE(OUT,CFMT) I,(A(I,J),J=START,STOP)
      IF(STOP.EQ.N) RETURN
      START=STOP+1
      GO TO 11
3001  FORMAT(' ')
      END
C.hr HEADER
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     HEADER   7/16/87
C
C     PURPOSE
C     TITLE OUTPUT.
C
C     USAGE
C     CALL HEADER(TITLE)
C
C     ARGUMENTS
C     TITLE - A CHARACTER STRING ENCLOSED IN QUOTES.  THE LAST CHARACTER
C             MUST BE AN UNDERSCORE.  THE STRING IS MADE UP OF AN UN-
C             LIMITED NUMBER OF CONCATONATED LINES.  A LINE IS MADE UP
C             OF 0-68 CHARACTERS FOLLOWED BY A SLASH.
C
C     REMARK
C     AN EXAMPLE IS:
C     CALL HEADER('//LINE 1/LINE 2/LINE 3///_')
C
C     COMMENT
C     THE DEFAULT LINESIZE IS 133, 132 PLUS CARRIAGE CONTROL CHARACTER.
C     THE USAGE:
C     COMMON /ZLNSIZ/ LNSIZE
C     LNSIZE=80
C     WILL CHANGE THE LINESIZE TO 80.  LINESIZES BETWEEN 72 AND 133 ARE
C     PERMITTED.
C
C
      SUBROUTINE HEADER(TITLE)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /ZLNSIZ/ LNSIZE
      CHARACTER*1 DIGIT(10),LDUMMY(8)
      CHARACTER*1 TITLE(1),PRINT(68),SLASH,UNC,LTR,BLANK
      CHARACTER*4 DUMMY(18)
      CHARACTER*8 FMT1(11),FMT2(3),FMT3(3),RDUMMY
      CHARACTER*24 CFMT2,CFMT3
      CHARACTER*88 CFMT1
      EQUIVALENCE (FMT1(1),CFMT1),(FMT2(1),CFMT2),(FMT3,CFMT3)
      EQUIVALENCE (RDUMMY,LDUMMY(1)),(DUMMY(1),PRINT(1))
      SAVE
      DATA FMT1/'(1X,    ','X,70H***',8*'********','***  )  '/
      DATA FMT2/'(1X,    ','X,1H*,17','A4,1H*) '/
      DATA FMT3/'(1X,    ','X,1H*,68','X,1H*)  '/
      DATA DIGIT/'0','1','2','3','4','5','6','7','8','9'/
      DATA SLASH/'/'/,UNC/'_'/,BLANK/' '/
      IOUT=3
      LNSIZ=LNSIZE
      IF((LNSIZE.LT.72).OR.(LNSIZE.GT.133)) LNSIZ=133
      RDUMMY=FMT1(1)
      IPAD=(LNSIZ-72)/2+1
      IPAD10=IPAD/10
      IPAD1=IPAD-10*IPAD10
      LDUMMY(7)=DIGIT(IPAD10+1)
      LDUMMY(8)=DIGIT(IPAD1+1)
      FMT1(1)=RDUMMY
      FMT2(1)=RDUMMY
      FMT3(1)=RDUMMY
      WRITE(IOUT,3001)
      WRITE(IOUT,3001)
      WRITE(IOUT,CFMT1)
      M=0
      L=0
10    M=M+L
      L=0
      DO 20 I=1,69
      LTR=TITLE(M+I)
      IF((LTR.EQ.SLASH).OR.(LTR.EQ.UNC)) GO TO 25
20    L=I
25    IF((L.EQ.0).AND.(LTR.EQ.UNC)) GO TO 50
      IF(L.EQ.0) GO TO 60
      IPAD=(68-L)/2
      DO 30 I=1,68
30    PRINT(I)=BLANK
      DO 40 I=1,L
      JADDR=I+IPAD
40    PRINT(JADDR)=TITLE(M+I)
      WRITE(IOUT,CFMT2) (DUMMY(I),I=1,17)
      L=L+1
      LTR=TITLE(M+L+1)
      IF(LTR.EQ.UNC) GO TO 50
      GO TO 10
50    WRITE(IOUT,CFMT1)
      RETURN
60    WRITE(IOUT,CFMT3)
      L=L+1
      GO TO 10
3001  FORMAT(' ')
      END
C.hr SPACER
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     SPACER   5/10/73
C
C     PURPOSE
C     SPACE OUTPUT
C
C     USAGE
C     CALL SPACER(I)
C
C     ARGUMENT
C     I - NUMBER OF LINES TO BE SKIPPED.  IF I=0 THEN THE PRINTER IS
C         ADVACED TO A NEW PAGE.
C         INTEGER*4
C
C
      SUBROUTINE SPACER(I)
      IF(I.EQ.0) WRITE(3,301)
      IF(I/66.GT.0) WRITE(3,301)
      IF(I.LE.0) RETURN
      M=MOD(I,66)
      DO 10 J=1,M
10    WRITE(3,302)
      RETURN
301   FORMAT('1')
302   FORMAT(' ')
      END
C.hr DSWEEP
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DSWEEP   10/4/72
C
C     PURPOSE
C     INVERT A POSITIVE DEFINITE MATRIX IN PLACE BY A DIAGONAL SWEEP.
C
C     USAGE
C     CALL DSWEEP(A,N,EPS,IER)
C
C     ARGUMENTS
C     A   - SYMMETRIC POSITIVE DEFINITE N BY N MATRIX STORED COLUMNWISE
C           (STORAGE MODE OF 0).  ON RETURN CONTAINS THE INVERSE OF A
C           STORED COLUMNWISE.
C           REAL*8
C     N   - NUMBER OF ROWS AND COLUMNS OF A.
C           INTEGER
C     EPS - INPUT CONSTANT USED AS A RELATIVE TOLERANCE IN TESTING FOR
C           DEGENERATE RANK.  A REASONABLE VALUE FOR EPS IS 1.D-13.
C           REAL*8
C     IER - ERROR PARAMETER CODED AS FOLLOWS
C           IER=0  NO ERROR, RANK OF A IS N.
C           IER.GT.0  A IS SINGULAR, RANK OF A IS N-IER.
C           INTEGER
C
C     REMARK
C     IF IER.GT.0 THEN DSWEEP RETURNS A G-INVERSE.
C
C     REFERENCE
C     SCHATZOFF, M. ET AL.  EFFICIENT CALCULATION OF ALL POSSIBLE REG-
C     RESSIONS.  TECHNOMETRICS, 10. 769-79 (NOVEMBER 1968)
C
C     PROGRAMMER
C     DR. A. RONALD GALLANT
C     DEPARTMANT OF STATISTICS
C     NORTH CAROLINA STATE UNIVERSITY
C     RALEIGH, NORTH CAROLINA  27696-8203
C
C
C
      SUBROUTINE DSWEEP(A,N,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 A(1)
      TOL=0.D0
      DO 5 I=1,N
      II=N*(I-1)+I
      TEST=A(II)
5     IF(TEST.GT.TOL) TOL=TEST
      TOL=TOL*EPS
      IER=0
      DO 50 K=1,N
      KK=N*(K-1)+K
      AKK=A(KK)
      IF(AKK.GT.TOL) GO TO 20
      DO 10 J=K,N
      KJ=N*(J-1)+K
10    A(KJ)=0.D0
      IF(K.EQ.1) GO TO 16
      KLESS1=K-1
      DO 15 I=1,KLESS1
      IK=KK-I
15    A(IK)=0.D0
16    IER=IER+1
      GO TO 50
20    D=1.D0/AKK
      DO 25 I=1,N
      DO 25 J=I,N
      IF((I.EQ.K).OR.(J.EQ.K)) GO TO 25
      IJ=N*(J-1)+I
      IF(I.LT.K) AIK=A(N*(K-1)+I)
      IF(I.GT.K) AIK=A(N*(I-1)+K)
      IF(K.LT.J) AKJ=A(N*(J-1)+K)
      IF(K.GT.J) AKJ=-A(N*(K-1)+J)
      A(IJ)=A(IJ)-AIK*AKJ*D
25    CONTINUE
      DO 30 J=K,N
      KJ=N*(J-1)+K
30    A(KJ)=A(KJ)*D
      IF(K.EQ.1) GO TO 36
      KLESS1=K-1
      DO 35 I=1,KLESS1
      IK=KK-I
35    A(IK)=-A(IK)*D
36    A(KK)=D
50    CONTINUE
      DO 55 I=1,N
      DO 55 J=I,N
      IF(I.EQ.J) GO TO 55
      IJ=N*(J-1)+I
      JI=N*(I-1)+J
      A(JI)=A(IJ)
55    CONTINUE
      RETURN
      END
C.hr DGMPRD
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMPRD  6/15/83
C
C     PURPOSE
C     MULTIPLY TWO MATRICES: R = A*B
C
C     USAGE
C     CALL DGMPRD(A,B,R,N,M,L)
C
C     ARGUMENTS
C     A - INPUT N BY M MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         REAL*8
C     B - INPUT M BY L MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         REAL*8
C     R - OUTPUT N BY L MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         REAL*8
C     N - NUMBER OF ROWS IN A
C         INTEGER*4
C     M - NUMBER OF COLUMNS IN A AND NUMBER OF ROWS IN A
C         INTEGER*4
C     L - NUMBER OF COLUMNS IN B
C         INTEGER*4
C
C     COMMENT
C     THE ROUTINE IS WRITTEN SO AS TO MINIMIZE THE NUMBER OF PAGE
C     FAULTS IN A VIRTUAL ENVIRONMENT.
C
      SUBROUTINE DGMPRD(A,B,R,N,M,L)
      REAL*8 A(1),B(1),R(1)
C     A(N,M),B(M,L),R(N,L)
      NL=N*L
      DO 10 I=1,NL
10    R(I)=0.D0
      DO 20 J=1,L
      DO 20 K=1,M
      DO 20 I=1,N
20    R(N*(J-1)+I)=R(N*(J-1)+I)+A(N*(K-1)+I)*B(M*(J-1)+K)
      RETURN
      END
C.hr DMCVT
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DMCVT    2/28/74
C
C     PURPOSE
C     CONVERT SYMMETRIC MATRICES FROM STORAGE MODE 0 TO MODE 1 OR FROM
C     STORAGE MODE 1 TO MODE 0.
C
C     USAGE
C     CALL DMCVT(A,R,N,MA,MR)
C
C     ARGUMENTS
C     A  - INPUT SYMMETRIC N BY N MATRIX STORED IN MODE MA.
C          REAL*8
C     R  - OUTPUT COPY OF A STORED IN MODE MR.
C          REAL*8
C     N  - ORDER OF A AND R.
C          INTEGER*4
C     MA - INPUT STORAGE MODE OF A.  MA MUST BE 0 OR 1 ONLY.
C          INTEGER*4
C     MR - INPUT STORAGE MODE OF R.  MR MUST BE 0 OR 1 ONLY.
C          INTEGER*4
C
C     REMARK
C     FOR IN-PLACE CONVERSTION THE USAGE IS:
C     CALL DMCVT(A,A,N,MA,MR)
C
C
C
      SUBROUTINE DMCVT(A,R,N,MA,MR)
      REAL*8 A(1),R(1)
      IF((MA.EQ.0).AND.(MR.EQ.0)) GO TO 40
      IF((MA.EQ.0).AND.(MR.EQ.1)) GO TO 10
      IF((MA.EQ.1).AND.(MR.EQ.0)) GO TO 20
      IF((MA.EQ.1).AND.(MR.EQ.1)) GO TO 50
      RETURN
C
C     CONVERT MODE 0  TO MODE 1
C
10    DO 15 J=1,N
      DO 15 I=1,J
15    R((J*(J-1))/2+I)=A(N*(J-1)+I)
      RETURN
C
C     CONVERT MODE 1 TO MODE 0
C
20    DO 25 JJ=1,N
      J=N-(JJ-1)
      DO 25 II=1,J
      I=J-(II-1)
25    R(N*(J-1)+I)=A((J*(J-1))/2+I)
      DO 30 J=1,N
      DO 30 I=1,J
30    R(N*(I-1)+J)=R(N*(J-1)+I)
      RETURN
C
C     COPY MODE 0 TO MODE 0
C
40    NN=N*N
      DO 45 I=1,NN
45    R(I)=A(I)
      RETURN
C
C     COPY MODE 1 TO MODE 1
C
50    NN=(N*(N-1))/2+N
      DO 55 I=1,NN
55    R(I)=A(I)
      RETURN
      END
C.hr DMFSD
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     SUBROUTINE DMFSD
C
C     PURPOSE
C     FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
C
C     USAGE
C     CALL DMFSD(A,N,EPS,IER)
C
C     ARGUMENTS
C     A   - UPPER TRIANGULAR PART OF GIVEN SYMMETRIC POSITIVE DEFINITE N
C           BY N MATRIX; ON RETURN, UPPER TRIANGULAR R, A=R'R.
C           REAL*8
C     N   - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX
C           INTEGER*4
C     EPS - INPUT CONSTANT WHICH IS USED AS RELATIVE TOLERANCE FOR TEST
C           ON LOSS OF SIGNIFICANCE, A REASONABLE VALUE IS 1.D-13.
C           REAL*8
C     IER - ERROR PARAMETER, CODED AS FOLLOWS
C           IER=0  - NO ERROR
C           IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER N OR
C                    BECAUSE SOME RADICAND IS NON-POSITIVE (MATRIX A IS
C                    NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS OF
C                    SIGNIFICANCE)
C           IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE.  THE
C                    RADICAND FORMED AT FACTORIZATION STEP K+1 WAS STILL
C                    POSITIVE BUT NO LONGER GREATER THAN
C                    ABS(EPS*A(K+1,K+1)
C
C     REMARK
C     THIS IS A COPY OF SSP ROUTINE DMFSD
C
      SUBROUTINE DMFSD(A,N,EPS,IER)
      DOUBLE PRECISION EPS,DPIV,DSUM
      DOUBLE PRECISION A(1)
      IF(N-1) 12,1,1
    1 IER=0
      KPIV=0
      DO 11 K=1,N
      KPIV=KPIV+K
      IND=KPIV
      LEND=K-1
      TOL=ABS(SNGL(EPS*A(KPIV)))
      DO 11 I=K,N
      DSUM=0.D0
      IF(LEND) 2,4,2
    2 DO 3 L=1,LEND
      LANF=KPIV-L
      LIND=IND-L
    3 DSUM=DSUM+A(LANF)*A(LIND)
    4 DSUM=A(IND)-DSUM
      IF(I-K) 10,5,10
    5 IF(SNGL(DSUM)-TOL) 6,6,9
    6 IF(DSUM) 12,12,7
    7 IF(IER) 8,8,9
    8 IER=K-1
    9 DPIV=DSQRT(DSUM)
      A(KPIV)=DPIV
      DPIV=1.D0/DPIV
      GO TO 11
   10 A(IND)=DSUM*DPIV
   11 IND=IND+I
      RETURN
   12 IER=-1
      RETURN
      END
C.hr FACTOR
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     FACTOR      12/20/87
C
C     PURPOSE
C     FACTOR A SYMMETRIC POSITIVE DEFINITE MATRIX A = RR' WHERE R IS
C     UPPER TRIANGULAR
C
C     USAGE
C     CALL FACTOR(A,DEL,M,EPS,IER)
C
C     ARGUMENTS
C     A   - ON INPUT CONTAINS THE UPPER TRIANGULAR PART OF A STORED
C           COLUMNWISE (STORAGE MODE 1).
C           ON OUTPUT CONTAINS THE UPPER TRIANGULAR PART OF R STORED
C           COLUMNWISE (STORAGE MODE 1).
C           THE (I,J) ELEMENT OF A IS IN PHYSICAL LOCATION J*(J-1)/2+I
C           WITH I.LE.J; THE TOTAL LENGTH IS M*(M+1)/2. THE SAME FOR R.
C           REAL*8,
C     DEL - ON OUTPUT CONTAINS THE JACOBIAN OF R WITH RESPECT TO A, AN
C           M*(M+1)/2 BY M*(M+1)/2 MATRIX STORED COLUMNWISE (STORAGE
C           MODE 0).  THE DERIVATIVE OF THE (I,J) ELEMENT OF R WITH
C           RESPECT TO THE (K,L) ELEMENT OF A WITH I.LE.J AND K.LE.L
C           IS IN CONCEPTUAL LOCATION (J*(J-1)/2+I,L*(L-1)/2+K) AND IN
C           PHYSICAL LOCATION (M*(M+1)/2)*(J*(J-1)/2+I-1)+L*(L-1)/2+K.
C           REAL*8
C     M   - INPUT, THE NUMBER OF ROWS (COLUMNS) IN A.
C           INTEGER*4
C     EPS - INPUT CONSTANT BETWEEN ZERO AND ONE WHICH IS USED AS
C           RELATIVE TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE, A
C           REASONABLE VALUE IS 1.D-13.
C           REAL*8
C     IER - OUTPUT ERROR PARAMETER, CODED AS FOLLOWS:
C           IER=0  - NO ERROR.
C           IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER M,EPS OR
C                    BECAUSE SOME RADICAND IS NOT POSITIVE (MATRIX A IS
C                    NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS OF
C                    SIGNIFICANCE).
C           IER=J  - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE.  THE
C                    RADICAND FORMED AT FACTORIZATION STEP J+1 WAS STILL
C                    POSITIVE BUT NO LONGER GREATER THAN
C                    DABS(EPS*A(J*(J+1)/2)).
C
C     NOTE
C     THE IMPLICIT LOWER PART OF A IS THE TRANSPOSE OF THE UPPER PART.
C     THE IMPLICIT LOWER PART OF R IS ZEROES BELOW THE DIAGONAL.
C
      SUBROUTINE FACTOR(A,DEL,M,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 A(M*(M+1)/2),DEL(M*(M+1)/2,M*(M+1)/2)
C
      IF (M.LE.0) THEN
      IER=-1
      RETURN
      ENDIF
C
      IF ((EPS.LE.0).OR.(EPS.GE.1.D0)) THEN
      IER=-1
      RETURN
      ENDIF
C
      LA=M*(M+1)/2
C
C     PUT DEL TO THE IDENTITY
C
      DO 10 I=1,LA
      DO 10 J=1,LA
10    DEL(I,J)=0.D0
      DO 20 I=1,LA
20    DEL(I,I)=1.D0
C
C     FACTORIZE; THIS IS A CHOLESKY FACTORIZATION GOTTEN BY WORKING FROM
C     SOUTHEAST TO NORTHWEST TO GET AN UPPER TRIANGULAR MATRIX INSTEAD
C     OF THE TRADITIONAL LOWER TRIANGULAR MATRIX GOTTEN BY WORKING FROM
C     NORTHWEST TO SOUTHEAST.
C
      IER=0
      DO 100 J1=0,M-1
      J=M-J1
      JJ=J*(J+1)/2
C
C     TAKE THE SQUARE ROOT OF AJJ, PUT IT IN AJJ
C
      R=A(JJ)
      IF (R.GT.0.D0) THEN
      R=DSQRT(R)
      A(JJ)=R
      DO 30 L=1,M
      DO 30 K=1,L
      KL=L*(L-1)/2+K
30    DEL(JJ,KL)=(0.5D0/R)*DEL(JJ,KL)
      ELSE
      IER=-1
      RETURN
      ENDIF
C
C     WE ARE FINISHED IF JJ.EQ.1, EQUIVALENTLY IF J.EQ.1
C
      IF (JJ.EQ.1) RETURN
C
C     DIVIDE THE ELEMENTS OF THE COLUMN ABOVE AJJ BY AJJ
C
      DO 50 I=1,J-1
      A(J*(J-1)/2+I)=A(J*(J-1)/2+I)/R
      DO 40 L=1,M
      DO 40 K=1,L
      KL=L*(L-1)/2+K
40    DEL(J*(J-1)/2+I,KL)=DEL(J*(J-1)/2+I,KL)/R
     &                    -(A(J*(J-1)/2+I)/R)*DEL(JJ,KL)
50    CONTINUE
C
C     SWEEP THE NORTHEAST SUBMATRIX
C
      SAVE=A((J-1)*(J)/2)
      DO 70 J0=1,J-1
      DO 70 I0=1,J0
      IJ=J0*(J0-1)/2+I0
      A(IJ)=A(IJ)-A(J*(J-1)/2+I0)*A(J*(J-1)/2+J0)
      DO 60 L=1,M
      DO 60 K=1,L
      KL=L*(L-1)/2+K
60    DEL(IJ,KL)=DEL(IJ,KL)
     &           -A(J*(J-1)/2+I0)*DEL(J*(J-1)/2+J0,KL)
     &           -DEL(J*(J-1)/2+I0,KL)*A(J*(J-1)/2+J0)
70    CONTINUE
C
C     TOLERANCE CHECK
C
      IF(SAVE.LE.0) THEN
      IER=-1
      RETURN
      ENDIF
      IF(A((J-1)*(J)/2).LT.DABS(EPS*SAVE)) IER=J1
C
100   CONTINUE
C
      RETURN
      END
C.hr DGMABP
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMABP   8/20/73
C
C     PURPOSE
C     POSTMULTIPLY A BY THE TRANSPOSE OF B (R=AB').
C
C     USAGE
C     CALL DGMABP(A,B,R,N,M,L)
C
C     ARGUMENTS
C     A - INPUT N BY M MATRIX STORED COLUMNWISE (STORAGE MODE 0).
C         REAL*8
C     B - INPUT L BY M MATRIX STORED COLUMNWISE (STORAGE MODE 0).
C         REAL*8
C     R - OUTPUT N BY L MATRIX STORED COLUMNWISE (STORAGE MODE 0).
C         REAL*8
C     N - NUMBER OF ROWS OF A AND R.
C         INTEGER*4
C     M - NUMBER OF COLUMNS OF A AND B.
C         INTEGER*4
C     L - NUMBER OF ROWS OF B AND COLUMNS OF R.
C         INTEGER*4
C
C
      SUBROUTINE DGMABP(A,B,R,N,M,L)
      REAL*8 A(N,M),B(L,M),R(N,L)
      DO 10 I=1,N
      DO 10 J=1,L
      R(I,J)=0.D0
      DO 10 K=1,M
10    R(I,J)=R(I,J)+A(I,K)*B(J,K)
      RETURN
      END
C.hr DRINV
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     SUBROUTINE DRINV
C
C     PURPOSE
C     INVERT AN UPPPER TRIANGULAR MATRIX.
C
C     USAGE
C     CALL DRINV(A,N)
C
C     ARGUMENTS
C     A - AN N*(N+1)/2-VECTOR CONTAINING THE UPPER TRIANGLE OF AN UPPER
C         TRIANGULAR MATRIX STORED COLUMNWISE (STORAGE MODE 1).
C         CONTAINS THE MATRIX TO BE INVERTED ON ENTRY AND THE INVERSE,
C         WHICH IS UPPER TRIANGULAR, ON RETURN; REAL*8
C     N - THE NUMBER OF ROWS AND COLUMNS IN A
C         INPUT, INTEGER*4
C
C     REMARK
C     THE ROUTINE ASSUMES THAT THE MATRIX TO BE INVERTED WAS RETURNED BY
C     DMFSD AND DOES NO ERROR CHECKING.  THE CODE WAS LIFTED FROM THE
C     MIDSECTION OF SSP ROUTINE DSINV
C
      SUBROUTINE DRINV(A,N)
      DIMENSION A(1)
      DOUBLE PRECISION A,DIN,WORK
    1 IPIV=N*(N+1)/2
      IND=IPIV
      DO 6 I=1,N
      DIN=1.D0/A(IPIV)
      A(IPIV)=DIN
      MIN=N
      KEND=I-1
      LANF=N-KEND
      IF(KEND) 5,5,2
    2 J=IND
      DO 4 K=1,KEND
      WORK=0.D0
      MIN=MIN-1
      LHOR=IPIV
      LVER=J
      DO 3 L=LANF,MIN
      LVER=LVER+1
      LHOR=LHOR+L
    3 WORK=WORK+A(LVER)*A(LHOR)
      A(J)=-WORK*DIN
    4 J=J-MIN
    5 IPIV=IPIV-MIN
    6 IND=IND-1
    9 RETURN
      END
C.hr DGMFAC
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMFAC      1/16/89 (Drew's twelfth birthday)
C
C     PURPOSE
C     FACTOR A SYMMETRIC, POSITIVE DEFINITE MATRIX AS A = RR' AND FACTOR
C     ITS INVERSE AS INV(A)=P'P WHERE R AND P ARE UPPER TRIANGULAR
C
C     USAGE
C     CALL DGMFAC(A,M,R,P,EPS,IER)
C
C     SUBROUTINES CALLED
C     DMFSD, DRINV
C
C     ARGUMENTS
C     A   - AN M BY M SYMMETRIC, POSITIVE DEFINITE MATRIX STORED
C           COLUMNWISE (STORAGE MODE 0).
C           INPUT, REAL*8
C     R   - FACTOR OF A, A=RR', AN M BY M MATRIX STORED COLUMNWISE
C           (STORAGE MODE 0).
C           OUTPUT, REAL*8
C     P   - FACTOR OF A INVERSE, INV(A)=P'P, AN M BY M MATRIX STORED
C           COLUMNWISE (STORAGE MODE 0).
C           OUTPUT, REAL*8
C     M   - THE NUMBER OF ROWS (COLUMNS) IN A.
C           INPUT, INTEGER*4
C     EPS - CONSTANT BETWEEN ZERO AND ONE WHICH IS USED AS RELATIVE
C           TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE, A REASONABLE
C           VALUE IS 1.D-13.
C           INPUT, REAL*8
C     IER - OUTPUT ERROR PARAMETER, CODED AS FOLLOWS:
C           IER=0  - NO ERROR.
C           IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER M,EPS OR
C                    BECAUSE SOME RADICAND IS NOT POSITIVE (MATRIX A IS
C                    NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS OF
C                    SIGNIFICANCE).
C           IER=J  - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE.  THE
C                    RADICAND FORMED AT FACTORIZATION STEP J+1 WAS STILL
C                    POSITIVE BUT NO LONGER GREATER THAN
C                    DABS(EPS*A(J*(J+1)/2)).
C
      SUBROUTINE DGMFAC(A,M,R,P,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 A(1),R(1),P(1)
C
C     PERMUTE THE ROWS AND COLUMNS OF A AND STORE IN P IN STORAGE MODE 1
C     THE PERMUTAIONS ARE BECAUSE WE FACTOR AS UPPER TRAINGULAR INSTEAD
C     OF THE STANDARD LOWER TRIANGULAR CHOLESKY FACTORIZATION.
C
      DO 10 J=1,M
      JJ=M+1-J
      DO 10 I=J,M
      II=M+1-I
10    P(JJ*(JJ-1)/2+II)=A(M*(J-1)+I)
C
C     COMPUTE THE CHOLESKY FACTORIZATION
C
      CALL DMFSD(P,M,EPS,IER)
      IF (IER.LT.0) RETURN
C
C     UNDO THE PERMUTATIONS AND PUT THE FACTOR IN R IN STORAGE MODE 1
C
      DO 20 J=1,M
      JJ=M+1-J
      DO 20 I=1,J
      II=M+1-I
20    R(J*(J-1)/2+I)=P(II*(II-1)/2+JJ)
C
C     COPY R TO P AND INVERT P
C
      DO 30 I=1,M*(M+1)/2
30    P(I)=R(I)
      CALL DRINV(P,M)
C
C     R AND P ARE IN COMPRESSED FORM, STORAGE MODE 1.  EXPAND TO
C     STORAGE MODE 0
C
      DO 40 JJ=1,M
      J=M+1-JJ
      DO 40 II=JJ,M
      I=M+1-II
      R(M*(J-1)+I)=R(J*(J-1)/2+I)
40    P(M*(J-1)+I)=P(J*(J-1)/2+I)
C
      DO 50 J=1,M
      DO 50 I=J,M
      IF(I.NE.J) THEN
         R(M*(J-1)+I)=0.D0
         P(M*(J-1)+I)=0.D0
         ENDIF
50    CONTINUE
C
      RETURN
      END
C.hr DGMABA
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMABA  6/20/87 (ANNIVERSARY)
C
C     PURPOSE
C     COMPUTE THE QUADRATIC FORM: R = A'BA.
C
C     USAGE
C     CALL DGMABA(A,B,R,N,M)
C
C     ARGUMENTS
C     A - INPUT N BY M MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         REAL*8
C     B - INPUT N BY N MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         REAL*8
C     R - OUTPUT M BY M MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         REAL*8
C     N - NUMBER OF ROWS IN A
C         INTEGER*4
C     M - NUMBER OF COLUMNS IN A
C         INTEGER*4
C
      SUBROUTINE DGMABA(A,B,R,N,M)
      REAL*8 A(1),B(1),R(1),ACC
C     A(N,M),B(N,N),R(M,M)
      DO 20 JJ=1,M
      DO 20 II=1,JJ
      ACC=0.D0
      DO 10 J=1,N
      DO 10 I=1,N
10    ACC=ACC+A(N*(JJ-1)+I)*A(N*(II-1)+J)*B(N*(J-1)+I)
      R(M*(II-1)+JJ)=ACC
20    R(M*(JJ-1)+II)=ACC
      RETURN
      END
C.hr DGMTRA
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMTRA 8/23/73
C
C     PURPOSE
C     DOUBLE PRECISION VERSION OF SSP GMTRA
C
      SUBROUTINE DGMTRA(A,R,N,M)
      REAL*8 A(1),R(1)
      DO 10 I=1,N
      DO 10 J=1,M
 10   R((I-1)*M+J)=A((J-1)*N+I)
      RETURN
      END
C.hr DGMADD
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMADD  10/16/71
C
C     PURPOSE
C     DOUBLE PRECISION VERSION OF SSP GMADD
C
      SUBROUTINE DGMADD(A,B,R,N,M)
      REAL*8 A(1),B(1),R(1)
      NM=N*M
      DO 10 I=1,NM
10    R(I)=A(I)+B(I)
      RETURN
      END
C.hr DGMSUB
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMSUB  10/16/71
C
C     PURPOSE
C     DOUBLE PRECISION VERSION OF SSP GMSUB
C
      SUBROUTINE DGMSUB(A,B,R,M,N)
      REAL*8 A(1),B(1),R(1)
      NM=N*M
      DO 10 I=1,NM
10    R(I)=A(I)-B(I)
      RETURN
      END
C.hr DCOND2
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DCOND2    5/15/74
C
C     PURPOSE
C     IMPROVE ACCURACY OF DSWEEP WHEN DIAGONAL ELEMENTS OF A ARE NOT OF
C     THE SAME ORDER OF MAGNITUDE.
C
C     USAGE
C     CALL DCOND2(A,N,S,0)
C     CALL DSWEEP(A,N,EPS,IER)
C     CALL DCOND2(A,N,S,1)
C
C     ARGUMENTS
C     A   - SYMMETRIC POSITIVE DEFINITE N BY N MATRIX STORED COLUMNWISE.
C           (STORAGE MODE 0)
C           REAL*8
C     N   - NUMBER OF ROWS AND COLUMNS OF A.  LENGTH OF S.
C           INTEGER*4
C     S   - WORK VECTOR OF LENGTH N.
C           REAL*8
C     ISW - INTEGER SWITCH.
C           INTEGER*4
C
C     REMARK
C     THE USAGE
C     CALL DCOND2(A,N,S,0)
C     WILL RETURN A CORRELATION MATRIX IN A AND THE STANDARD ERRORS IN
C     S WHEN A IS A VARIANCE-COVARIANCE MATRIX.
C
C
C
      SUBROUTINE DCOND2(A,N,S,ISW)
      REAL*8 A(N,N),S(N)
      IF(ISW.EQ.1) GO TO 20
      DO 10 I=1,N
10    S(I)=DSQRT(A(I,I))
20    DO 30 I=1,N
      DO 30 J=1,N
30    IF(S(I)*S(J).NE.0.D0)
     &A(I,J)=A(I,J)/(S(I)*S(J))
      RETURN
      END
