      LOGICAL FUNCTION OPENN(A,NX,N,LENGTH) 
  
*** INITIALIZE NORMAL EQUATION MATRIX 
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/NR,N1,N2,ISTATE 
  
*** SET STATE TO UNREDUCED
  
      ISTATE=0
      NR=N
      N1=N+1
      N2=N1+N 
  
*** ZERO THE MATRIX 
  
      IF(NX(N1).LE.LENGTH) THEN 
        OPENN=.TRUE.
        DO 1 I=1,NX(N1) 
    1   A(I)=0.D0 
      ELSE
        OPENN=.FALSE. 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION GETA(II,JJ,VAL,A,NX) 
  
*** RETURN ELEMENT OF MATRIX
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
*** CONVERT EXTERNAL NUMBERS TO INTERNAL
  
      IF(II.LT.N1) THEN 
        I=NX(N2+II) 
      ELSE
        I=N1
      ENDIF 
      IF(JJ.LT.N1) THEN 
        J=NX(N2+JJ) 
      ELSE
        J=N1
      ENDIF 
  
      IF(I.LE.J) THEN 
        IF(J.EQ.N1) THEN
          VAL=A(NX(N)+I)
          GETA=.TRUE. 
        ELSEIF(I.EQ.N) THEN 
          VAL=A(NX(N))
          GETA=.TRUE. 
        ELSE
          INDEX=NX(I)-I+J 
          IF(INDEX.GE.NX(I+1)) THEN 
            VAL=0.D0
            GETA=.FALSE.
          ELSE
            VAL=A(INDEX)
            GETA=.TRUE. 
          ENDIF 
        ENDIF 
      ELSE

*** 20-dec-89 mod

        IF(I.EQ.N1) THEN
          VAL=A(NX(N)+J)
          GETA=.TRUE. 
        ELSEIF(J.EQ.N) THEN 
          VAL=A(NX(N))
          GETA=.TRUE. 
        ELSE
          INDEX=NX(J)-J+I 
          IF(INDEX.GE.NX(J+1)) THEN 
            VAL=0.D0
            GETA=.FALSE.
          ELSE
            VAL=A(INDEX)
            GETA=.TRUE. 
          ENDIF 
        ENDIF 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION PUTA(II,JJ,VAL,A,NX) 
  
*** STORE ELEMENT OF MATRIX 
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
*** CONVERT EXTERNAL NUMBERS TO INTERNAL
  
      IF(II.LT.N1) THEN 
        I=NX(N2+II) 
      ELSE
        I=N1
      ENDIF 
      IF(JJ.LT.N1) THEN 
        J=NX(N2+JJ) 
      ELSE
        J=N1
      ENDIF 
  
      IF(I.LE.J) THEN 
        IF(J.EQ.N1) THEN
          A(NX(N)+I)=VAL
          PUTA=.TRUE. 
        ELSEIF(I.EQ.N) THEN 
          A(NX(N))=VAL
          PUTA=.TRUE. 
        ELSE
          INDEX=NX(I)-I+J 
          IF(INDEX.GE.NX(I+1)) THEN 
            PUTA=.FALSE.
          ELSE
            A(INDEX)=VAL
            PUTA=.TRUE. 
          ENDIF 
        ENDIF 
      ELSE

*** 20-dec-89 mod

        IF(I.EQ.N1) THEN
          A(NX(N)+J)=VAL
          PUTA=.TRUE. 
        ELSEIF(J.EQ.N) THEN 
          A(NX(N))=VAL
          PUTA=.TRUE. 
        ELSE
          INDEX=NX(J)-J+I 
          IF(INDEX.GE.NX(J+1)) THEN 
            PUTA=.FALSE.
          ELSE
            A(INDEX)=VAL
            PUTA=.TRUE. 
          ENDIF 
        ENDIF 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION ADDOBS(C,IC,LEN,OBS,WT,A,NX) 
  
*** ACCUMULATE OBSERVATIONS INTO NORMALS
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION C(LEN),IC(LEN)
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
*** LOOP OVER COEFFICIENTS
  
      A(NX(N1))=A(NX(N1))+OBS*WT*OBS
      DO 2 I=1,LEN
      II=NX(N2+IC(I)) 
      CIWT=C(I)*WT
      A(NX(N)+II)=A(NX(N)+II)+CIWT*OBS
      DO 2 J=I,LEN
      JJ=NX(N2+IC(J)) 
  
*** CHECK IF COEFFICIENTS WITHIN PROFILE
  
      IF(.NOT.IN(II,JJ,NX,INDEX)) THEN
        ADDOBS=.FALSE.
        RETURN
      ENDIF 
  
*** ACCUMULATE
  
  2   A(INDEX)=A(INDEX)+CIWT*C(J) 
  
*** ACCUMULATION SUCCESSFUL 
  
      ADDOBS=.TRUE. 
  
      RETURN
      END 
      LOGICAL FUNCTION ADDCOR(C,IC,OBS,SIGOBS,LUNK,LOBS,A,NX,IFLAG) 
  
*** ACCUMULATE CORRELATED OBSERVATIONS INTO NORMALS 
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION C(*),IC(LUNK),SIGOBS(LOBS,LOBS),OBS(LOBS) 
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
      DATA TOL/1.D-10/
  
*** CHOLESKY FACTOR WEIGHT MATRIX 
  
      DO 1 I=1,LOBS 
      I1=I-1
      IF(I1.GT.0) THEN
        DO 2 K=1,I1 
    2   SIGOBS(I,I)=SIGOBS(I,I)-SIGOBS(K,I)*SIGOBS(K,I) 
      ENDIF 
  
*** SINGULARITY TEST
  
      IF(SIGOBS(I,I).LT.TOL) THEN 
        ADDCOR=.FALSE.
        IFLAG=1 
        RETURN
      ENDIF 
  
*** RESUME CHOLESKY FACTOR
  
      SIGOBS(I,I)=DSQRT(SIGOBS(I,I))
      IF(I+1.LE.LOBS) THEN
        DO 6 J=I+1,LOBS 
        IF(I1.GT.0) THEN
          DO 4 K=1,I1 
    4     SIGOBS(I,J)=SIGOBS(I,J)-SIGOBS(K,I)*SIGOBS(K,J) 
        ENDIF 
        SIGOBS(I,J)=SIGOBS(I,J)/SIGOBS(I,I) 
    6   SIGOBS(J,I)=SIGOBS(I,J) 
      ENDIF 
    1 CONTINUE
  
*** TRANSFORM OBSERVATIONS
  
      CALL CHSOLV(SIGOBS,OBS,LOBS)
  
*** TRANSFORM COEFFICIENTS
  
      DO 3 I=1,LUNK 
      J=(I-1)*LOBS+1
    3 CALL CHSOLV(SIGOBS,C(J),LOBS) 
  
*** ACCUMULATE TRANSFORMED OBSERVATIONS -- UNIT WEIGHT
  
      DO 100 K=1,LOBS 
  
*** LOOP OVER COEFFICIENTS
  
        A(NX(N1))=A(NX(N1))+OBS(K)*OBS(K) 
        DO 5 I=1,LUNK 
        II=NX(N2+IC(I)) 
        CI=C((I-1)*LOBS+K)
        A(NX(N)+II)=A(NX(N)+II)+CI*OBS(K) 
        DO 5 J=I,LUNK 
        JJ=NX(N2+IC(J)) 
  
*** CHECK IF COEFFICIENTS WITHIN PROFILE
  
        IF(.NOT.IN(II,JJ,NX,INDEX)) THEN
          ADDCOR=.FALSE.
          IFLAG=2 
          RETURN
        ENDIF 
  
*** ACCUMULATE
  
    5   A(INDEX)=A(INDEX)+CI*C((J-1)*LOBS+K)
  100 CONTINUE
  
*** ACCUMULATION SUCCESSFUL 
  
      ADDCOR=.TRUE. 
      IFLAG=0 
  
      RETURN
      END 
      SUBROUTINE CHSOLV(A,X,LENG) 
  
*** SOLVE AY=X FOR Y AND STORE RESULT IN X (GIVEN CHOLESKY FACTOR, A) 
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LENG,LENG),X(LENG)
  
*** SUBSTITUTION
  
      DO 5 I=1,LENG 
      I1=I-1
      TEMP=0.D0 
      IF(I1.GT.0) THEN
        DO 6 K=1,I1 
    6   TEMP=TEMP+A(K,I)*X(K) 
      ENDIF 
    5 X(I)=(X(I)-TEMP)/A(I,I) 
  
      RETURN
      END 
      LOGICAL FUNCTION SOLVE(A,NX,TOL,ISING,GSING,LSING)
  
*** FACTOR AND SOLVE A LINEAR SYSTEM
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION ISING(*),GSING(*) 
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
      DATA BIG,TINY/1.D30,1.D-30/ 
  
      NM1=N-1 
      LSING=0 
      IF(ISTATE.EQ.0) THEN
        ISTATE=1
        SOLVE=.TRUE.
  
*** CHOLESKY FACTOR 
  
        DO 1 I=N,1,-1 
          IDIAG=NX(I) 
          DIAG=A(IDIAG) 
          IF(I.NE.N) THEN 
            ITOP=NX(I+1)-1
            ILEN=ITOP-IDIAG 
            IF(ILEN.GT.0) THEN
              ICOL=ILEN+I 
              DO 2 J=I-1,1,-1 
                JDIAG=NX(J) 
                JTOP=NX(J+1)-1
                JLEN=JTOP-JDIAG 
                JCOL=JLEN+J 
                IF(JCOL.GT.I) THEN
                  JI=JDIAG-J+I
                  IF(ICOL.LT.JCOL) THEN 
                    NP=ILEN 
                    IP=ITOP-NP+1
                    JP=JTOP-(JCOL-ICOL)-NP+1
                  ELSE
                    NP=JLEN+J-I 
                    JP=JTOP-NP+1
                    IP=ITOP-(ICOL-JCOL)-NP+1
                  ENDIF 
                  CALL INNERP(A(IP),A(JP),NP,PROD)
                  A(JI)=A(JI)-PROD
                ENDIF 
    2         CONTINUE
            ENDIF 
            IA=IDIAG+1
            IB=NX(I+1)-1
            DO 4 JI=IA,IB 
    4       A(IDIAG)=A(IDIAG)-A(JI)*A(JI) 
          ENDIF 
  
*** SINGULARITY TEST
  
          IF(DIAG.LE.0.D0) THEN 
            GOOGE=0.D0
          ELSE
            GOOGE=A(IDIAG)/DIAG 
          ENDIF 
          IF(GOOGE.LT.TOL) THEN 
            LSING=LSING+1 
            ISING(LSING)=NX(N1+I) 
            IF(GOOGE.LE.0.D0) THEN
              GSING(LSING)=0.D0 
            ELSE
              GSING(LSING)=GOOGE
            ENDIF 
            SOLVE=.FALSE. 
            A(IDIAG)=BIG
          ENDIF 
  
*** RESUME CHOLESKY FACTOR
  
          DIAG=DSQRT(A(IDIAG))
          A(IDIAG)=DIAG 
          RDIAG=1.D0/DIAG 
          DO 5 J=1,I-1
            JI=NX(J)-J+I
            IF(JI.LT.NX(J+1)) THEN
              TEMP=A(JI)*RDIAG
              IF(DABS(TEMP).LT.TINY) THEN 
                A(JI)=0.D0
              ELSE
                A(JI)=TEMP
              ENDIF 
            ENDIF 
    5     CONTINUE
    1   CONTINUE
  
*** COMPLETE CHOLESKY FACTOR
  
*** FORWARD SUBSTITUTION
  
        NM1=N-1 
        A(NX(N1)+N)=A(NX(N)+N)/A(NX(N)) 
        DO 7 I=NM1,1,-1 
          TEMP=0.D0 
          IA=NX(I)+1
          IB=NX(I+1)-1
          K=I 
          DO 8 IK=IA,IB 
            K=K+1 
            TEMP=TEMP+A(IK)*A(NX(N1)+K) 
    8     CONTINUE
    7   A(NX(N1)+I)=(A(NX(N)+I)-TEMP)/A(NX(I))
  
*** BASEMENT WINDOW 
  
        TEMP=A(NX(N1))
        DO 10 I=1,N 
          Y=A(NX(N1)+I) 
   10   TEMP=TEMP-Y*Y 
        A(NX(N1))=TEMP
  
*** BACKWARD SOLUTION 
  
        DO 9 I=1,N
        I1=I-1
        TEMP=0.D0 
        DO 6 K=1,I1 
          KI=NX(K)-K+I
          IF(KI.LT.NX(K+1)) TEMP=TEMP+A(KI)*A(NX(N)+K)
    6   CONTINUE
    9   A(NX(N)+I)=(A(NX(N1)+I)-TEMP)/A(NX(I))
      ELSE
        SOLVE=.FALSE. 
      ENDIF 
  
      RETURN
      END 
      SUBROUTINE INNERP(X,Y,N,PROD) 
  
*** COMPUTE INNER PRODUCT 
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(*),Y(*) 
  
      PROD=0.D0 
      DO 1 I=1,N
    1 PROD=PROD+X(I)*Y(I) 
  
      RETURN
      END 
      LOGICAL FUNCTION RESOLV(A,NX) 
  
*** SOLVE THE FACTORED NORMAL EQUATIONS FOR A NEW R.H. SIDE 
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
      IF(ISTATE.NE.1) THEN
        RESOLV=.FALSE.
      ELSE
  
*** FORWARD SUBSTITUTION
  
        RESOLV=.TRUE. 
  
        NM1=N-1 
        A(NX(N1)+N)=A(NX(N)+N)/A(NX(N)) 
        DO 7 I=NM1,1,-1 
          TEMP=0.D0 
          IA=NX(I)+1
          IB=NX(I+1)-1
          K=I 
          DO 8 IK=IA,IB 
            K=K+1 
            TEMP=TEMP+A(IK)*A(NX(N1)+K) 
    8     CONTINUE
    7   A(NX(N1)+I)=(A(NX(N)+I)-TEMP)/A(NX(I))
  
*** BASEMENT WINDOW 
  
        TEMP=A(NX(N1))
        DO 10 I=1,N 
          Y=A(NX(N1)+I) 
   10   TEMP=TEMP-Y*Y 
        A(NX(N1))=TEMP
  
*** BACKWARD SOLUTION 
  
        DO 9 I=1,N
        I1=I-1
        TEMP=0.D0 
        DO 6 K=1,I1 
          KI=NX(K)-K+I
          IF(KI.LT.NX(K+1)) TEMP=TEMP+A(KI)*A(NX(N)+K)
    6   CONTINUE
    9   A(NX(N)+I)=(A(NX(N1)+I)-TEMP)/A(NX(I))
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION INVERT(A,NX) 
  
*** INVERT THE CHOLESKY FACTOR
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
      NM1=N-1 
      IF(ISTATE.NE.1) THEN
        INVERT=.FALSE.
      ELSE
        ISTATE=2
        INVERT=.TRUE. 
  
*** PUT CHOLESKY RECIPROCALS ON DIAGONAL
  
        DO 1 I=1,N
    1   A(NX(I))=1.D0/A(NX(I))
  
*** COMPUTE GAUSSIAN FACTOR FROM CHOLESKY FACTOR
  
        DO 2 I=N,2,-1 
          DIAG=A(NX(I)) 
          DO 7 J=1,I
            JI=NX(J)-J+I
            IF(JI.LT.NX(J+1)) A(JI)=DIAG*A(JI)
    7     CONTINUE
    2   CONTINUE
        A(1)=A(1)*A(1)
  
*** PROCEED BY COLUMNS--TOP TO BOTTOM 
  
        DO 6 I=2,N
  
*** SAVE I-TH COL IN SCRATCH SPACE
  
          DO 3 J=1,I-1
            JI=NX(J)-J+I
            IF(JI.LT.NX(J+1)) A(NX(N1)+J)=-A(JI)
    3     CONTINUE
  
*** GET INVERSE OF OFF DIAGONAL ELEMENTS OF I-TH COL
  
          DO 5 J=1,I-1
            JI=NX(J)-J+I
            IF(JI.LT.NX(J+1)) THEN
              A(JI)=0.D0
              DO 4 K=1,I-1
                KI=NX(K)-K+I
                IF(KI.LT.NX(K+1)) THEN
                  YK=A(NX(N1)+K)
                  IF(K.LT.J) THEN 
                    KJ=NX(K)-K+J
                    IF(KJ.LT.NX(K+1)) A(JI)=A(JI)+YK*A(KJ)
                  ELSE
                    JK=NX(J)-J+K
                    IF(JK.LT.NX(J+1)) A(JI)=A(JI)+YK*A(JK)
                  ENDIF 
                ENDIF 
    4         CONTINUE
            ENDIF 
    5     CONTINUE
  
*** GET INVERSE OF DIAGONAL AT I-TH COL 
  
          DO 8 J=1,I-1
            JI=NX(J)-J+I
            IF(JI.LT.NX(J+1)) THEN
              A(NX(I))=A(NX(I))+A(JI)*A(NX(N1)+J) 
            ENDIF 
    8     CONTINUE
    6   CONTINUE
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION PROP(C,IC,LEN,VAR,A,NX,IFLAG)
  
*** COMPUTE VARIANCE BY LINEAR ERROR PROPAGATION
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION C(LEN),IC(LEN)
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
*** VERIFY PARAMETER COVARIANCE AVAILABLE 
  
      IF(ISTATE.NE.2) THEN
        IFLAG=1 
        PROP=.FALSE.
      ELSE
        VAR=0.D0
  
        DO 20 I=1,LEN 
        II=NX(N2+IC(I)) 
        VAR=VAR+C(I)*C(I)*A(NX(II)) 
        I1=I+1
        IF(I1.LE.LEN) THEN
          C2=C(I)+C(I)
          DO 10 J=I1,LEN
          JJ=NX(N2+IC(J)) 
          IF(.NOT.IN(II,JJ,NX,INDEX)) THEN
            IFLAG=2 
            PROP=.FALSE.
            RETURN
          ENDIF 
   10     VAR=VAR+C2*C(J)*A(INDEX)
        ENDIF 
   20   CONTINUE
  
        IFLAG=0 
        PROP=.TRUE. 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION PROPCV(C1,IC1,LEN1,C2,IC2,LEN2,COV,A,NX,IFLAG) 
  
*** COMPUTE COVARIANCE BY LINEAR ERROR PROPAGATION
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GETA
      DIMENSION C1(LEN1),IC1(LEN1)
      DIMENSION C2(LEN2),IC2(LEN2)
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
*** VERIFY PARAMETER COVARIANCES AVAILABLE
  
      IF(ISTATE.NE.2) THEN
        IFLAG=1 
        PROPCV=.FALSE.
      ELSE
        COV=0.D0
        DO 20 I=1,LEN1
        CI=C1(I)
        II=IC1(I) 
        DO 10 J=1,LEN2
          JJ=IC2(J) 
          IF(.NOT.GETA(II,JJ,SIG,A,NX)) THEN
            IFLAG=2 
            PROPCV=.FALSE.
            RETURN
          ENDIF 
   10     COV=COV+CI*C2(J)*SIG
   20   CONTINUE
  
        IFLAG=0 
        PROPCV=.TRUE. 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION CORR(A,NX) 
  
*** COMPUTE CORRELATION MATRIX FROM COVARIANCE
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IN
      DIMENSION A(*),NX(*)
      COMMON/DENNIS/N,N1,N2,ISTATE
  
      NM1=N-1 
      IF(ISTATE.NE.2) THEN
        CORR=.FALSE.
      ELSE
        ISTATE=3
        CORR=.TRUE. 
  
*** CONVERT VARIANCES TO STANDARD DEVIATIONS
  
        DO 1 I=1,N
    1   A(NX(I))=DSQRT(A(NX(I)))
  
*** DIVIDE EACH OFF-DIAGONAL ROW ELEMENT BY DIAGONAL
  
        DO 5 I=1,NM1
        RDIAG=1.D0/A(NX(I)) 
        IA=NX(I)+1
        IB=NX(I+1)-1
        IF(IA.LE.IB) THEN 
          DO 2 IJ=IA,IB 
    2     A(IJ)=A(IJ)*RDIAG 
        ENDIF 
    5   CONTINUE
  
*** DIVIDE EACH OFF-DIAGONAL COLUMN ELEMENT BY DIAGONAL 
  
        DO 3 J=N,2,-1 
        RDIAG=1.D0/A(NX(J)) 
        DO 3 I=1,J-1
        IF(IN(I,J,NX,IJ)) A(IJ)=A(IJ)*RDIAG 
    3   CONTINUE
  
*** PLACE ONE'S ALONG DIAGONAL
  
        DO 4 I=1,N
    4   A(NX(I))=1.D0 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION IN(I,J,NX,INDEX) 
  
*** RETURN ELEMENT INDEX OF MATRIX PROFILE
*** CAUTION: I AND J ARE INTERNAL NUMBERS !!
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION NX(*) 
      COMMON/DENNIS/N,N1,N2,ISTATE
  
      IF(I.LE.J) THEN 
        IF(J.EQ.N1) THEN
          INDEX=NX(N)+I 
          IN=.TRUE. 
        ELSEIF(I.EQ.N) THEN 
          INDEX=NX(N) 
          IN=.TRUE. 
        ELSE
          INDEX=NX(I)-I+J 
          IF(INDEX.GE.NX(I+1)) THEN 
            IN=.FALSE.
          ELSE
            IN=.TRUE. 
          ENDIF 
        ENDIF 
      ELSE
        IF(J.EQ.N) THEN 
          INDEX=NX(N) 
          IN=.TRUE. 
        ELSEIF(I.EQ.N1) THEN
          INDEX=NX(N)+J 
          IN=.TRUE. 
        ELSE
          INDEX=NX(J)-J+I 
          IF(INDEX.GE.NX(J+1)) THEN 
            IN=.FALSE.
          ELSE
            IN=.TRUE. 
          ENDIF 
        ENDIF 
      ENDIF 
  
      RETURN
      END 
      SUBROUTINE OPENG(N,NX,NXLENG) 
  
*** INITIALIZE LINKED NEIGHBOR LIST 
  
      DIMENSION NX(*) 
      COMMON/KATHY/NR,N1,N2,IEMX,IFREE
  
*** COMPUTE INDICIES
  
      NR=N
      N1=N+1
      N2=N1+N 
  
*** COMPUTE MAXIMUM NUMBER OF EDGES 
*** INSURE IEMX (2*EDGES) IS EVEN 
  
      IEMX=((NXLENG-N-N-1)/6)*2 
  
      DO 1 I=1,N
    1 NX(I)=-I
      IFREE=1 
  
      RETURN
      END 
      LOGICAL FUNCTION ADDCON(IC,LEN,NX)
  
*** ADD OBSERVATION EQUATION CONNECTIONS TO NEIGHBOR LIST 
  
      LOGICAL ADDG
      DIMENSION IC(*),NX(*) 
      COMMON/KATHY/N,N1,N2,IEMX,IFREE 
  
      IF(LEN.LE.1) THEN 
        WRITE(*,2) LEN
    2   FORMAT('LEN ERROR IN ADDCON',I5)
        STOP
      ENDIF 
      L1=LEN-1
  
*** LOOP OVER ALL CONNECTIONS 
  
      DO 1 I=1,L1 
      DO 1 J=I+1,LEN
      IF(.NOT.ADDG(IC(I),IC(J),NX(1),NX(N1),NX(N1+IEMX))) THEN
        ADDCON=.FALSE.
        RETURN
      ENDIF 
    1 CONTINUE
      ADDCON=.TRUE. 
  
      RETURN
      END 
      LOGICAL FUNCTION ADDG(I,J,IHEAD,NBR,LINK) 
  
*** ADD A CONNECTION
  
      DIMENSION IHEAD(*),NBR(*),LINK(*) 
      LOGICAL CONECT
      COMMON/KATHY/N,N1,N2,IEMX,IFREE 
  
      ADDG=.TRUE. 
  
*** ADD CONNECTION IF IT DOES NOT EXIST 
  
      IF(.NOT.CONECT(I,J,IPOINT,IHEAD,NBR,LINK)) THEN 
        IF(IPOINT.EQ.0) THEN
          IHEAD(I)=IFREE
        ELSE
          IF(LINK(IPOINT).NE.-I) THEN 
            WRITE(*,1) I,J,IPOINT,LINK(IPOINT)
    1       FORMAT(4I5,' -- GRAPH SUBROUTINE LINK ERROR') 
            STOP
          ELSE
            LINK(IPOINT)=IFREE
          ENDIF 
        ENDIF 
        NBR(IFREE)=J
        LINK(IFREE)=-I
        IFREE=IFREE+1 
        IF(IFREE.GT.IEMX) THEN
          ADDG=.FALSE.
          RETURN
        ENDIF 
      ENDIF 
  
*** ADD CONNECTION IN OTHER DIRECTION 
  
      IF(.NOT.CONECT(J,I,IPOINT,IHEAD,NBR,LINK)) THEN 
        IF(IPOINT.EQ.0) THEN
          IHEAD(J)=IFREE
        ELSE
          IF(LINK(IPOINT).NE.-J) THEN 
            WRITE(*,1) I,J,IPOINT,LINK(IPOINT)
            STOP
          ELSE
            LINK(IPOINT)=IFREE
          ENDIF 
        ENDIF 
        NBR(IFREE)=I
        LINK(IFREE)=-J
        IFREE=IFREE+1 
        IF(IFREE.GT.IEMX) THEN
          ADDG=.FALSE.
          RETURN
        ENDIF 
      ENDIF 
  
      RETURN
      END 
      LOGICAL FUNCTION CONECT(I,J,IPNTX,IHEAD,NBR,LINK) 
  
*** IS CONNECTION IN THE GRAPH? 
*** IF NOT, FIND POINTER TO END OF LINKED LIST
  
      DIMENSION IHEAD(*),NBR(*),LINK(*) 
  
      IPNTX=0 
      IPOINT=IHEAD(I) 
  
    1 IF(IPOINT.LT.0) THEN
        CONECT=.FALSE.
        RETURN
      ELSE
        IF(NBR(IPOINT).EQ.J) THEN 
          CONECT=.TRUE. 
          RETURN
        ENDIF 
        IPNTX=IPOINT
        IPOINT=LINK(IPOINT) 
        GO TO 1 
      ENDIF 
  
      RETURN
      END 
      SUBROUTINE REORDR(NX) 
  
*** REORDER UNKNOWNS INTO CUTHILL-MCKEE ORDER 
*** FOR MINIMUM PROFILE, ROW STORAGE, UPPER TRIANGLE
  
      DIMENSION NX(*) 
      COMMON/KATHY/N,N1,N2,IEMX,IFREE 
  
*** PRODUCE ADJACENCY STRUCTURE FROM LINKED LIST
  
      CALL TRNSLT(NX(1),NX(N1),NX(N1+IEMX),NX(N1+3*IEMX),NX(N1+2*IEMX)) 
  
*** PRODUCE CUTHILL-MCKEE ORDER 
  
      CALL GENCM(N,NX(N1+3*IEMX),NX(N1+2*IEMX),NX(N1+1),NX(2*N1), 
     *            NX(3*N+2))
  
*** INVERT ORDER
  
      DO 777 I=1,N
      J=NX(N1+I)
  777 NX(N2+J)=I
  
*** PRODUCE PROFILE 
  
      CALL PROFIL(NX(1),NX(N1+3*IEMX),NX(N1+2*IEMX))
  
      RETURN
      END 
      SUBROUTINE TRNSLT(IHEAD,NBR,LINK,NBRPNT,NBRLST) 
  
*** TRANSLATE LINKED NEIGHBOR LIST INTO ADJACENCY STRUCTURE 
  
      DIMENSION IHEAD(*),NBR(*),LINK(*),NBRPNT(*),NBRLST(*) 
      COMMON/KATHY/N,N1,N2,IEMX,IFREE 
  
      IFREE=1 
      DO 100 I=1,N
      IPOINT=IHEAD(I) 
      NBRPNT(I)=IFREE 
    1 IF(IPOINT.GT.0) THEN
        NBRLST(IFREE)=NBR(IPOINT) 
        IPOINT=LINK(IPOINT) 
        IFREE=IFREE+1 
        GO TO 1 
      ENDIF 
  100 CONTINUE
      NBRPNT(N1)=IFREE
  
      RETURN
      END 
      SUBROUTINE PROFIL(NX,NBRPNT,NBRLST) 
  
*** FIND VARIABLE COLUMN PROFILE FROM ADJACENCY STRUCTURE 
  
      DIMENSION NX(*),NBRPNT(*),NBRLST(*) 
      COMMON/KATHY/N,N1,N2,IEMX,IFREE 
  
*** STORE ROW WIDTHS IN NX(2 THRU N1) 
  
      DO 20 IR=1,N
      I=NX(N2+IR) 
      LBEGIN=NBRPNT(IR) 
      LEND=NBRPNT(IR+1)-1 
      LROW=0
      IF(LBEGIN.LE.LEND) THEN 
        DO 10 IPT=LBEGIN,LEND 
        JC=NBRLST(IPT)
        J=NX(N2+JC) 
        IF(J.GT.I) THEN 
          LEN=J-I 
          IF(LEN.GT.LROW) LROW=LEN
        ENDIF 
   10   CONTINUE
      ENDIF 
   20 NX(I+1)=LROW+1
  
*** CONVERT ROW WIDTHS TO PROFILE 
  
      NX(1)=1 
      DO 30 I=2,N 
   30 NX(I)=NX(I)+NX(I-1) 
      NX(N1)=NX(N)+N1 
  
      RETURN
      END 
      SUBROUTINE GENCM(NEQNS,NBRPNT,NBRLST,NEWORD,MASK,LEVPTR)
  
*** FIND CUTHILL-MCKEE ORDERING FOR A GENERAL GRAPH 
*** FOR EACH COMPONENT, ORDERING IS FOUND BY SUBROUTINE CM
  
      DIMENSION NBRLST(*),MASK(*),NEWORD(*),LEVPTR(*),NBRPNT(*) 
  
*** INITIALIZE COMPONENT MASK 
  
      DO 100 I=1,NEQNS
  100 MASK(I)=1 
      NUM=1 
  
*** LOOP OVER ALL EQUATIONS 
  
      DO 200 I=1,NEQNS
      IF(MASK(I).EQ.0) GO TO 200
      IROOT=I 
  
*** FIND PSEUDO-PERIPHERAL NODE 
  
      CALL FNROOT(IROOT,NBRPNT,NBRLST,MASK,NLVL,LEVPTR,NEWORD(NUM)) 
  
*** NUMBER COMPONENT IN CUTHILL-MCKEE ORDER 
  
      CALL CM(IROOT,NBRPNT,NBRLST,MASK,NEWORD(NUM),LCOMP,LEVPTR)
      NUM=NUM+LCOMP 
      IF(NUM.GT.NEQNS) RETURN 
  200 CONTINUE
  
      RETURN
      END 
      SUBROUTINE FNROOT(IROOT,NBRPNT,NBRLST,MASK,NLVL,LEVPTR,LS)
  
*** FIND PSEUDO-PERIPHERAL NODE 
  
      DIMENSION NBRLST(*),LS(*),MASK(*),LEVPTR(*),NBRPNT(*) 
  
*** DETERMINE LEVEL STRUCTURE ROOTED AT IROOT 
  
      CALL ROOTLS(IROOT,NBRPNT,NBRLST,MASK,NLVL,LEVPTR,LS)
      LCOMP=LEVPTR(NLVL+1)-1
      IF(NLVL.EQ.1.OR.NLVL.EQ.LCOMP) RETURN 
  
*** PICK A NODE WITH MINIMUM DEGREE FROM THE LAST LEVEL 
  
  100 JSTRT=LEVPTR(NLVL)
      MINDEG=LCOMP
      IROOT=LS(JSTRT) 
      IF(LCOMP.EQ.JSTRT) GO TO 400
        DO 300 J=JSTRT,LCOMP
        NODE=LS(J)
        NDEG=0
        KSTRT=NBRPNT(NODE)
        KSTOP=NBRPNT(NODE+1)-1
        DO 200 K=KSTRT,KSTOP
        NABOR=NBRLST(K) 
        IF(MASK(NABOR).GT.0) NDEG=NDEG+1
  200   CONTINUE
        IF(NDEG.GE.MINDEG) GO TO 300
        IROOT=NODE
        MINDEG=NDEG 
  300   CONTINUE
  
*** GENERATE ITS ROOTED LEVEL STRUCTURE 
  
  400 CALL ROOTLS(IROOT,NBRPNT,NBRLST,MASK,NUNLVL,LEVPTR,LS)
      IF(NUNLVL.LE.NLVL) RETURN 
      NLVL=NUNLVL 
      IF(NLVL.LT.LCOMP) GO TO 100 
  
      RETURN
      END 
      SUBROUTINE ROOTLS(IROOT,NBRPNT,NBRLST,MASK,NLVL,LEVPTR,LS)
  
*** GENERATE LEVEL STRUCTURE ROOTED AT IROOT
  
      DIMENSION NBRLST(*),LS(*),MASK(*),LEVPTR(*),NBRPNT(*) 
  
*** INITIALIZATION
  
      MASK(IROOT)=0 
      LS(1)=IROOT 
      NLVL=0
      LVLEND=0
      LCOMP=1 
  
*** LBEGIN IS POINTER TO BEGINNING OF CURRENT LEVEL 
*** LVLEND IS POINTER TO END OF CURRENT LEVEL 
  
  200 LBEGIN=LVLEND+1 
      LVLEND=LCOMP
      NLVL=NLVL+1 
      LEVPTR(NLVL)=LBEGIN 
  
*** GENERATE NEXT LEVEL BY FINDING ALL MASKED NEIGHBORS OF NODES IN 
*** THE CURRENT LEVEL 
  
      DO 400 I=LBEGIN,LVLEND
      NODE=LS(I)
      JSTRT=NBRPNT(NODE)
      JSTOP=NBRPNT(NODE+1)-1
      IF(JSTOP.LT.JSTRT) GO TO 400
      DO 300 J=JSTRT,JSTOP
      NBR=NBRLST(J) 
      IF(MASK(NBR).EQ.0) GO TO 300
      LCOMP=LCOMP+1 
      LS(LCOMP)=NBR 
      MASK(NBR)=0 
  300 CONTINUE
  400 CONTINUE
  
*** COMPUTE CURRENT LEVEL WIDTH--IF NONZERO, GENERATE NEXT LEVEL
  
      LVSIZE=LCOMP-LVLEND 
      IF(LVSIZE.GT.0) GO TO 200 
  
*** RESET MASK TO 1 FOR NODES IN THE LEVEL STRUCTURE
  
      LEVPTR(NLVL+1)=LVLEND+1 
      DO 500 I=1,LCOMP
      NODE=LS(I)
  500 MASK(NODE)=1
  
      RETURN
      END 
      SUBROUTINE CM(IROOT,NBRPNT,NBRLST,MASK,NEWORD,LCOMP,LDEG) 
  
*** NUMBER CONNECTED COMPONENT USING CUTHILL-MCKEE ALGORITHM
  
      DIMENSION NBRLST(*),LDEG(*),MASK(*),NEWORD(*),NBRPNT(*) 
  
*** FIND DEGREES OF NODES IN COMPONENT GIVEN BY MASK & IROOT
  
      CALL DEGREE(IROOT,NBRPNT,NBRLST,MASK,LDEG,LCOMP,NEWORD) 
      MASK(IROOT)=0 
      IF(LCOMP.LE.1) RETURN 
      LVLEND=0
      LNBR=1
  
*** LBEGIN POINTS TO BEGINNING OF CURRENT LEVEL 
*** LVLEND POINTS TO END OF CURRENT LEVEL 
  
  100 LBEGIN=LVLEND+1 
      LVLEND=LNBR 
      DO 600 I=LBEGIN,LVLEND
  
*** FOR EACH NODE IN CURRENT LEVEL
  
      NODE=NEWORD(I)
      JSTRT=NBRPNT(NODE)
      JSTOP=NBRPNT(NODE+1)-1
  
*** FIND UNNUMBERED NEIGHBORS OF NODE 
*** IFNBR POINTS TO FIRST UNNUMBERED NEIGHBOR OF CURRENT NODE IN NEWORD 
*** LNBR POINTS TO LAST UNNUMBERED NEIGHBOR OF CURRENT NODE IN NEWORD 
  
      IFNBR=LNBR+1
      DO 200 J=JSTRT,JSTOP
      NBR=NBRLST(J) 
      IF(MASK(NBR).EQ.0) GO TO 200
      LNBR=LNBR+1 
      MASK(NBR)=0 
      NEWORD(LNBR)=NBR
  200 CONTINUE
      IF(IFNBR.GE.LNBR) GO TO 600 
  
*** SORT NEIGHBORS OF NODE INCREASING BY DEGREE (LIN. INSERT.)
  
      K=IFNBR 
  300 L=K 
      K=K+1 
      NBR=NEWORD(K) 
  400 IF(L.LT.IFNBR) GO TO 500
      LPERM=NEWORD(L) 
      IF(LDEG(LPERM).LE.LDEG(NBR)) GO TO 500
      NEWORD(L+1)=LPERM 
      L=L-1 
      GO TO 400 
  500 NEWORD(L+1)=NBR 
      IF(K.LT.LNBR) GO TO 300 
  
*** END OF INSERTION SORT 
  
  600 CONTINUE
      IF(LNBR.GT.LVLEND) GO TO 100
  
      RETURN
      END 
      SUBROUTINE DEGREE(IROOT,NBRPNT,NBRLST,MASK,LDEG,LCOMP,LS) 
  
*** FIND DEGREES OF NODES OF CONNECTED COMPONENT
  
      DIMENSION NBRLST(*),LDEG(*),LS(*),MASK(*),NBRPNT(*) 
  
*** INITIALIZATION
*** NBRPNT USED AS TEMP MARKER TO SHOW WHICH NODES HAVE BEEN CONSIDERED 
  
      LS(1)=IROOT 
      NBRPNT(IROOT)=-NBRPNT(IROOT)
      LVLEND=0
      LCOMP=1 
  
*** LBEGIN IS POINTER TO BEGINNING OF CURRENT LEVEL 
*** LVLEND IS POINTER TO END OF CURRENT LEVEL 
  
  100 LBEGIN=LVLEND+1 
      LVLEND=LCOMP
  
*** FIND DEGREES OF NODES IN CURRENT LEVEL--ALSO GENERATE NEXT LEVEL
  
      DO 400 I=LBEGIN,LVLEND
      NODE=LS(I)
      JSTRT=-NBRPNT(NODE) 
      JSTOP=IABS(NBRPNT(NODE+1))-1
      IDEG=0
      IF(JSTOP.LT.JSTRT) GO TO 300
      DO 200 J=JSTRT,JSTOP
      NBR=NBRLST(J) 
      IF(MASK(NBR).EQ.0) GO TO 200
      IDEG=IDEG+1 
      IF(NBRPNT(NBR).LT.0) GO TO 200
      NBRPNT(NBR)=-NBRPNT(NBR)
      LCOMP=LCOMP+1 
      LS(LCOMP)=NBR 
  200 CONTINUE
  300 LDEG(NODE)=IDEG 
  400 CONTINUE
  
*** COMPUTE CURRENT LEVEL WIDTH--IF NONZERO, GENERATE ANOTHER LEVEL 
  
      LVSIZE=LCOMP-LVLEND 
      IF(LVSIZE.GT.0) GO TO 100 
  
*** RESET NBRPNT TO CORRECT SIGN
  
      DO 500 I=1,LCOMP
      NODE=LS(I)
  500 NBRPNT(NODE)=-NBRPNT(NODE)
  
      RETURN
      END 
      LOGICAL FUNCTION BAND(N,NX,NXLENG,LBAND)
  
*** GENERATE PROFILE AND INDEXING FOR BANDED STORAGE
  
      DIMENSION NX(*) 
  
      N1=N+1
      N2=N1+N 
      LENG=N2+N 
      IF(LBAND.LT.1) THEN 
        LBAND=1 
      ELSEIF(LBAND.GT.N) THEN 
        LBAND=N 
      ENDIF 
  
      IF(LENG.GT.NXLENG) THEN 
        BAND=.FALSE.
      ELSE
        DO 1 I=1,N
        NX(N1+I)=I
    1   NX(N2+I)=I
  
        NX(1)=1 
        LROW=N1 
        DO 2 I=2,N
          LROW=LROW-1 
          IF(LBAND.LT.LROW) THEN
            NX(I)=NX(I-1)+LBAND 
          ELSE
            NX(I)=NX(I-1)+LROW
          ENDIF 
    2     CONTINUE
        NX(N1)=NX(N)+N1 
        BAND=.TRUE. 
      ENDIF 
  
      RETURN
      END 
