*     Matrix subprograms
*
*     The subprograms are in alphabetical order, with the exception:
*      DOUBLE PRECISION subprograms are immediately behind the single
*      precision version
*
*     GMAB       product of general matrices
*        DGMAB   ...
*     GMPR       print general matrix
*     GMTR       transpose general matrix
*     SMAVA      product of general and symmetric matrix AVAT
*        DSMAVA   ... (error propagation)
*     SMHHL      Householder transformation
*        DSMHHL  ...
*     SMINV      symmetric matrix inversion and solution of equations
*        DSMINV  ...
*     SMTOG      convert symmetric to general matrix
*     SMPR       print symmetric matrix
*     SMPRV      print (symmetric) covariance matrix and vector
*        DSMPRV  ...
*     SMTOS      copy (part of) symmetric matrix to symmetric matrix
*     SMQU       square of symmetric matrix 
*     SMVA       product of symmetric and general matrix
*        DSMVA   ...
*     SOLVA      solution of symmetric matrix equation 
*     SOLVB      second part of solution ...
*        DSOLVA  ...
*        DSOLVB  ...
*     TQHHL      utility routine for SMHHL
*        DTQHHL  ...
*     VXY        scalar product of two vectors 
*

      SUBROUTINE GMAB(A,B,C,N,M,L)
*     Multiply general n-by-m matrix and general m-by-l matrix  to  form
*     general n-by-l matrix
*
*                  - -   - - -
*        CALL GMAB(A,B,C,N,M,L)         C  :=    A   *    B
*                      -               N*L      N*M      M*L
*
*
      REAL A(1),B(1),C(1)
      DOUBLE PRECISION SUM
*     ...
      IJ=0
      IK=0
      DO K=1,N
       DO I=1,L
        JK=I
        SUM=0.0
        DO J=1,M
         SUM=SUM+A(IJ+J)*B(JK)
         JK=JK+L
        END DO
        IK=IK+1
        C(IK)=SUM
       END DO
       IJ=IJ+M
      END DO
      END

      SUBROUTINE DGMAB(A,B,C,N,M,L)
*     Multiply general n-by-m matrix and general m-by-l matrix  to  form
*     general n-by-l matrix
*
*                   - -   - - -
*        CALL DGMAB(A,B,C,N,M,L)         C  :=    A   *    B
*                       -               N*L      N*M      M*L
*
*
      DOUBLE PRECISION A,B,C 
      DIMENSION A(*),B(*),C(*)
      DOUBLE PRECISION SUM
*     ...
      IJ=0
      IK=0
      DO K=1,N
       DO I=1,L
        JK=I
        SUM=0.0
        DO J=1,M
         SUM=SUM+A(IJ+J)*B(JK)
         JK=JK+L
        END DO
        IK=IK+1
        C(IK)=SUM
       END DO
       IJ=IJ+M
      END DO
      END

      SUBROUTINE GMPR(A,N,M,TEXT)
*     prints the N-by-M general matrix A
*                  - - - ----
*        CALL GMPR(A,N,M,TEXT)
*
      REAL A(*)
      CHARACTER*(*) TEXT
*
*     general matrix or vector
*
*     ...
      IF(N.EQ.1.OR.M.EQ.1) THEN
*        special case of vector
         WRITE(*,101) N,M,TEXT
         WRITE(*,102) (A(L),L=1,N*M)
      ELSE 
*        matrix
         WRITE(*,103) N,M,TEXT
         IJ=0
         DO I=1,N
          WRITE(*,104) I,(A(IJ+L),L=1,M)
          IJ=IJ+M
         END DO
      END IF
      WRITE(*,102)
  100 RETURN
  101 FORMAT(1X,I4,' *',I3,'  Vector ',A/)
  102 FORMAT(10X,10G12.5)
  103 FORMAT(1X,I4,' *',I3,'  Matrix ',A/)
  104 FORMAT(1X,I4,5X,10G12.5/(10X,10G12.5/))
      END

      SUBROUTINE GMTR(A,B,N,M)
*     Transpose N-by-M matrix A --> B
*
      REAL A(M,N),B(N,M)
      DO I=1,N
       DO J=1,M
        B(I,J)=A(J,I)
       END DO
      END DO
      END

      SUBROUTINE SMAVA(V,A,W,N,M)
*     Multiply symmetric N-by-N matrix from the left with general M-by-N
*     matrix and from the right with the transposed of the same  general
*     matrix  to  form  symmetric  M-by-M   matrix   (used   for   error
*     propagation).
*                   - -   - -
*        CALL SMAVA(V,A,W,N,M)
*                       -
*                                  T
*         W   =   A   *   V   *   A
*        M*M     M*N     N*N     N*M
*
*
*        where V = symmetric N-by-N matrix
*              A = general N-by-M matrix
*              W = symmetric M-by-M matrix
*
      REAL V(1),A(1),W(1)
      DOUBLE PRECISION CIK
*     ...
      DO I=1,(M*M+M)/2
       W(I)=0.0
      END DO
*
      IL=-N
      IJS=0
      DO I=1,M
       IJS=IJS+I-1
       IL=IL+N
       LKL=0
       DO K=1,N
        CIK=0.0
        LKL=LKL+K-1
        LK=LKL
        DO L=1,K
         LK=LK+1
         CIK=CIK+A(IL+L)*V(LK)
        END DO
        DO L=K+1,N
         LK=LK+L-1
         CIK=CIK+A(IL+L)*V(LK)
        END DO
        JK=K
        IJ=IJS
        DO J=1,I
         IJ=IJ+1
         W(IJ)=W(IJ)+CIK*A(JK)
         JK=JK+N
        END DO 
       END DO
      END DO  
      END

      SUBROUTINE DSMAVA(V,A,W,N,M)
*     Multiply symmetric N-by-N matrix from the left with general M-by-N
*     matrix and from the right with the transposed of the same  general
*     matrix  to  form  symmetric  M-by-M   matrix   (used   for   error
*     propagation).
*
*                    - -   - -
*        CALL DSMAVA(V,A,W,N,M)
*                        -
*                                  T
*         W   =   A   *   V   *   A
*        M*M     M*N     N*N     N*M
*
*
*        where V = symmetric N-by-N matrix
*              A = general N-by-M matrix
*              W = symmetric M-by-M matrix
*
      DOUBLE PRECISION V,A,W
      DIMENSION V(*),A(*),W(*)
      DOUBLE PRECISION CIK
*     ...
      DO I=1,(M*M+M)/2
       W(I)=0.0
      END DO
*
      IL=-N
      IJS=0
      DO I=1,M
       IJS=IJS+I-1
       IL=IL+N
       LKL=0
       DO K=1,N
        CIK=0.0
        LKL=LKL+K-1
        LK=LKL
        DO L=1,K
         LK=LK+1
         CIK=CIK+A(IL+L)*V(LK)
        END DO
        DO L=K+1,N
         LK=LK+L-1
         CIK=CIK+A(IL+L)*V(LK)
        END DO
        JK=K
        IJ=IJS
        DO J=1,I
         IJ=IJ+1
         W(IJ)=W(IJ)+CIK*A(JK)
         JK=JK+N
        END DO 
       END DO
      END DO  
      END

      SUBROUTINE SMHHL(V,N,D,U,IDIR)
*     Determination  of  eigenvalues  and  eigenvectors  by  Householder
*     method
*                 - -     ----
*     CALL SMHHL(V,N,D,U,IDIR)
*                     - -
*     Method:
*     Reduce symmetric matrix V to tridiagonal  form  using  Householder
*     transformations, and, by  call  to  TQHHL,  diagonalization   of
*     tridiagonal matrix
*        D(I)   = diagonal elements
*        U(I,J) = transformations matrix
*        IDIR   <  0  means eigenvalues in decreasing order
*        IDIR   >= 0  means eigenvalues in increasing order
*
*
*     SMHHL uses a work array of N words in COMMON/MATCOM/. For N > 400
*     the user has to define COMMON/MATCOM/ with N words.
*
*     Subroutine called: TQHHL
*
      REAL V(1),D(1),U(N,N),TOL/1.0E-15/
      COMMON/MATCOM/E(400)
*     ...
      IJ=0
      DO I=1,N
       DO J=1,I
        IJ=IJ+1
        U(I,J)=V(IJ)
        V(IJ)=0.0
       END DO
      END DO
*
      DO II=2,N
       I=N+2-II
       L=I-2
       F=U(I,I-1)
       G=0.0
       IF(L.NE.0) THEN
          DO K=1,L
           IF(ABS(U(I,K)).GT.TOL) G=G+U(I,K)*U(I,K)
          END DO
          H=G+F*F
       END IF
*
*     If  G  too  small  for  orthogonality  to   be   guaranteed,   the
*     transformation is skipped
*
       IF(G.LT.TOL) THEN
          E(I)=F
          H   =0.0
       ELSE
          L=L+1
          SH=SQRT(H)
          IF(F.GE.0.0) SH=-SH
          G=SH
          E(I)=SH
          H=H-F*G
          U(I,I-1)=F-G
          F=0.0
*
          DO J=1,L
           U(J,I)=U(I,J)/H
           G=0.0
*          form element of a u
           DO K=1,J
            IF(ABS(U(J,K)).GT.TOL.AND.ABS(U(I,K)).GT.TOL) THEN
               G=G+U(J,K)*U(I,K)
            END IF 
           END DO
           J1=J+1
           DO K=J1,L
            IF(ABS(U(K,J)).GT.TOL.AND.ABS(U(I,K)).GT.TOL) THEN
               G=G+U(K,J)*U(I,K)
            END IF
           END DO
           E(J)=G/H
           F=F+G*U(J,I)
          END DO
*         form k
          HH=F/(H+H)
*         form reduced a
          DO J=1,L
           F=U(I,J)
           E(J)=E(J)-HH*F
           G=E(J)
           DO K=1,J
            U(J,K)=U(J,K)-F*E(K)-G*U(I,K)
           END DO
          END DO
       END IF
       D(I)=H
      END DO
*
      D(1)=0.0
      E(1)=0.0

*     accumulation of transformation matrices
      DO I=1,N
       L=I-1
       IF(D(I).NE.0.0) THEN
          DO J=1,L
           G=0.0
           DO K=1,L
            G=G+U(I,K)*U(K,J)
           END DO
           DO K=1,L
            U(K,J)=U(K,J)-G*U(K,I)
           END DO
          END DO
       END IF
       D(I)=U(I,I)
       U(I,I)=1.0
       DO J=1,L
        U(I,J)=0.0
        U(J,I)=0.0
       END DO
      END DO
      CALL TQHHL(N,D,E,U,IDIR)
      END

      SUBROUTINE DSMHHL(V,N,D,U,IDIR)
*     Determination  of  eigenvalues  and  eigenvectors  by  Householder
*     method
*                 - -     ----
*     CALL DSMHHL(V,N,D,U,IDIR)
*                     - -
*     Method:
*     Reduce symmetric matrix V to tridiagonal  form  using  Householder
*     transformations, and, by  call  to  TQHHL,  diagonalization   of
*     tridiagonal matrix
*        D(I)   = diagonal elements
*        U(I,J) = transformations matrix
*        IDIR   <  0  means eigenvalues in decreasing order
*        IDIR   >= 0  means eigenvalues in increasing order
*
*
*     SMHHL uses a work array of N words in COMMON/MATCOM/. For N > 400
*     the user has to define COMMON/MATCOM/ with N words.
*
*     Subroutine called: DTQHHL

*
      DOUBLE PRECISION V(1),D(1),U(N,N),TOL,E
      COMMON/MATCOM/E(400)
      DATA TOL/1.0D-30/
*     ...
      IJ=0
      DO I=1,N
       DO J=1,I
        IJ=IJ+1
        U(I,J)=V(IJ)
        V(IJ)=0.0
       END DO
      END DO
*
      DO II=2,N
       I=N+2-II
       L=I-2
       F=U(I,I-1)
       G=0.0
       IF(L.NE.0) THEN
          DO K=1,L
           G=G+U(I,K)*U(I,K)
          END DO
          H=G+F*F
       END IF
*
*     If  G  too  small  for  orthogonality  to   be   guaranteed,   the
*     transformation is skipped
*
       IF(G.LT.TOL) THEN
          E(I)=F
          H   =0.0
       ELSE
          L=L+1
          SH=SQRT(H)
          IF(F.GE.0.0) SH=-SH
          G=SH
          E(I)=SH
          H=H-F*G
          U(I,I-1)=F-G
          F=0.0
*
          DO J=1,L
           U(J,I)=U(I,J)/H
           G=0.0
*          form element of a u
           DO K=1,J
            G=G+U(J,K)*U(I,K)
           END DO
           J1=J+1
           DO K=J1,L
            G=G+U(K,J)*U(I,K)
           END DO
           E(J)=G/H
           F=F+G*U(J,I)
          END DO
*         form k
          HH=F/(H+H)
*         form reduced a
          DO J=1,L
           F=U(I,J)
           E(J)=E(J)-HH*F
           G=E(J)
           DO K=1,J
            U(J,K)=U(J,K)-F*E(K)-G*U(I,K)
           END DO
          END DO
       END IF
       D(I)=H
      END DO
*
      D(1)=0.0
      E(1)=0.0

*     accumulation of transformation matrices

      DO I=1,N
       L=I-1
       IF(D(I).NE.0.0) THEN
          DO J=1,L
           G=0.0
           DO K=1,L
            G=G+U(I,K)*U(K,J)
           END DO
           DO K=1,L
            U(K,J)=U(K,J)-G*U(K,I)
           END DO
          END DO
       END IF
       D(I)=U(I,I)
       U(I,I)=1.0
       DO J=1,L
        U(I,J)=0.0
        U(J,I)=0.0
       END DO
      END DO
      CALL DTQHHL(N,D,E,U,IDIR)
      END

      SUBROUTINE SMINV(V,B,NARG,M,NRANK)
*     Obtain solution of a system of linear equations V *  X  =  B  with
*     symmetric matrix V and inverse (for M =  1)  or  matrix  inversion
*     only (for M = 0).
*
*                   - - - -
*        CALL SMINV(V,B,N,M,NRANK)
*                   - -     -----
*
*           V = symmetric N-by-N matrix in symmetric storage mode
*               V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, . . .
*               replaced by inverse matrix
*           B = N-vector   (for M = 0 use a dummy argument)
*               replaced by solution vector
*           M = see above
*
*
*     Method of solution is by elimination selecting the  pivot  on  the
*     diagonal each stage. The rank of the matrix is returned in  NRANK.
*     For NRANK ne N, all remaining  rows  and  cols  of  the  resulting
*     matrix V and the corresponding elements of  B  are  set  to  zero.
*     SMINV can be used for a dimension up to 100 (see INVCDR).
*
      REAL V(*),B(*),VKK,D,E
      COMMON/INVCDR/DR(2,100)
      DATA EPS/1.E-6/
*
*     construct table
*
      N=IABS(NARG)
      IF(NARG.GT.0) THEN
*        reset flags
         DO I=1,N
          DR(1,I)=1.0
         END DO 
         NI=N
      ELSE
*        call with negative N - matrix partially inverted
         NI=0
         DO I=1,N
          IF(DR(1,I).NE.0.0) NI=NI+1
         END DO
      END IF
      II=0
      DO I=1,N
       II=II+I
       DR(2,I)=ABS(V(II))
      END DO 
*     loop begin, with loop on all remaining rows/cols
      NRANK=N-NI
      DO I=1,NI
*      search for pivot and test for linearity and zero matrix
       K=0
       JJ=0
       KK=0
       VKK=0.0
       DO J=1,N
        JJ=JJ+J
        IF(DR(1,J).NE.0.0) THEN
           IF(ABS(V(JJ)).GT.VKK.AND.ABS(V(JJ)).GT.EPS*DR(2,J)) THEN
              VKK=ABS(V(JJ))
              K=J
              KK=JJ
           END IF   
        END IF
       END DO
c       if(k.ne.0) then
c       write(*,*) 'sminv',i,k,vkk,dr(2,k)
c       end if 
       IF(K.EQ.0.OR.VKK.LT.1.0E-16) THEN
*         no pivot found, clear of matrix
          IJ=0
          DO III=1,N
           IF(M.EQ.1.AND.DR(1,III).NE.0.0) B(III)=0.0
           DO J=1,III
            IJ=IJ+1
            IF(DR(1,III)+DR(1,J).NE.0.0) V(IJ)=0.0
            V(IJ)=-V(IJ)
           END DO 
          END DO
          GOTO 100
       END IF
*      preparation for elimination
       NRANK=NRANK+1
       DR(1,K)=0.0
       D=1.0/V(KK)
       IF(ABS(D).LT.1.0E-16) D=0.0
       V(KK)=-D
       IF(M.EQ.1) B(K)=B(K)*D
c       write(*,*) '      d and Bk ',D,B(K) 
       JK=KK-K
       JL=0
*      elimination 
       DO J=1,N
        IF(J.EQ.K) THEN
           JK=KK
           JL=JL+J
        ELSE
           IF(J.LT.K) THEN
              JK=JK+1
           ELSE
              JK=JK+J-1
           END IF
           E    =V(JK)
c           write(*,*) 'factor e',e,d 
           IF(ABS(E).LT.1.0E-16) E=0.0
           
           V(JK)=D*E
           IF(M.EQ.1) B(J) =B(J)-B(K)*E
           LK   =KK-K
           DO L=1,J
            JL=JL+1
            IF(L.EQ.K) THEN
               LK=KK
            ELSE
               IF(L.LT.K) THEN
                  LK=LK+1
               ELSE
                  LK=LK+L-1
               END IF
               V(JL)=V(JL)-V(LK)*E
            END IF
           END DO 
        END IF
       END DO 
      END DO 
*     change sign
      IJ=0
      DO I=1,N
       DO J=1,I
        IJ=IJ+1
        V(IJ)=-V(IJ)
       END DO
      END DO
  100 RETURN
      END

      SUBROUTINE DSMINV(V,B,NARG,M,NRANK)
*     Obtain solution of a system of linear equations V *  X  =  B  with
*     symmetric matrix V and inverse (for M =  1)  or  matrix  inversion
*     only (for M = 0).
*
*                   - - - -
*        CALL SMINV(V,B,N,M,NRANK)
*                   - -     -----
*
*           V = symmetric N-by-N matrix in symmetric storage mode
*               V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, . . .
*               replaced by inverse matrix
*           B = N-vector   (for M = 0 use a dummy argument)
*               replaced by solution vector
*           M = see above
*
*
*     Method of solution is by elimination selecting the  pivot  on  the
*     diagonal each stage. The rank of the matrix is returned in  NRANK.
*     For NRANK ne N, all remaining  rows  and  cols  of  the  resulting
*     matrix V and the corresponding elements of  B  are  set  to  zero.
*     SMINV can be used for a dimension up to 100 (see INVCDR).
*
      DOUBLE PRECISION V,B,VKK,D,E,EPS
      DIMENSION V(*),B(*)
      COMMON/INVCDR/DR(2,512)
      DATA EPS/1.D-6/
*
*     construct table
*
      N=IABS(NARG)
      IF(NARG.GT.0) THEN
*        reset flags
         DO I=1,N
          DR(1,I)=1.0
         END DO 
         NI=N
      ELSE
*        call with negative N - matrix partially inverted
         NI=0
         DO I=1,N
          IF(DR(1,I).NE.0.0) NI=NI+1
         END DO
      END IF
      II=0
      DO I=1,N
       II=II+I
       DR(2,I)=ABS(V(II))
      END DO 
*     loop begin, with loop on all remaining rows/cols
      NRANK=N-NI
      DO I=1,NI
*      search for pivot and test for linearity and zero matrix
       K=0
       JJ=0
       KK=0
       VKK=0.0
       DO J=1,N
        JJ=JJ+J
        IF(DR(1,J).NE.0.0) THEN
           IF(ABS(V(JJ)).GT.VKK.AND.ABS(V(JJ)).GT.EPS*DR(2,J)) THEN
              VKK=ABS(V(JJ))
              K=J
              KK=JJ
           END IF   
        END IF
       END DO
       IF(K.EQ.0) THEN
*         no pivot found, clear of matrix
          IJ=0
          DO III=1,N
           IF(M.EQ.1.AND.DR(1,III).NE.0.0) B(III)=0.0
           DO J=1,III
            IJ=IJ+1
            IF(DR(1,III)+DR(1,J).NE.0.0) V(IJ)=0.0
            V(IJ)=-V(IJ)
           END DO 
          END DO
          GOTO 100
       END IF
*      preparation for elimination
       NRANK=NRANK+1
       DR(1,K)=0.0
       D=1.0/V(KK)
       V(KK)=-D
       IF(M.EQ.1) B(K)=B(K)*D
       JK=KK-K
       JL=0
*      elimination 
       DO J=1,N
        IF(J.EQ.K) THEN
           JK=KK
           JL=JL+J
        ELSE
           IF(J.LT.K) THEN
              JK=JK+1
           ELSE
              JK=JK+J-1
           END IF
           E    =V(JK)
           V(JK)=D*E
           IF(M.EQ.1) B(J) =B(J)-B(K)*E
           LK   =KK-K
           DO L=1,J
            JL=JL+1
            IF(L.EQ.K) THEN
               LK=KK
            ELSE
               IF(L.LT.K) THEN
                  LK=LK+1
               ELSE
                  LK=LK+L-1
               END IF
               V(JL)=V(JL)-V(LK)*E
            END IF
           END DO 
        END IF
       END DO 
      END DO 
*     change sign
      IJ=0
      DO I=1,N
       DO J=1,I
        IJ=IJ+1
        V(IJ)=-V(IJ)
       END DO
      END DO
  100 RETURN
      END

      SUBROUTINE SMTOG(V,A,N)
*
*     copy symmetric n-by-n matrix V to general n-by-n matrix A 
*
*                   -   -
*        CALL SMTOG(V,A,N)        A := V
*                     -
      REAL V(1),A(1)
*     ...         
      IJV=0
      DO I=1,N
       IJA=I
       DO J=1,I
        IJV=IJV+1
        A(IJA)=V(IJV)
        A(IJA+(N-1)*(I-J))=V(IJV)
        IJA=IJA+N
       END DO
      END DO
      END 

      SUBROUTINE SMPR(V,N,TEXT)
*     prints the symmetric n-by-n matrix v
*                       ----
*        CALL SMPRT(V,N,TEXT)
*
      REAL V(*)
      CHARACTER*(*) TEXT
*     ...
      WRITE(*,101) N,N,TEXT
      IF(N.LE.0) GOTO 100
      II=0
      DO I=1,N
       WRITE(*,102) I,(V(II+J),J=1,I)
       II=II+I
      END DO
      WRITE(*,102)
  100 RETURN
  101 FORMAT(1X,I4,' by',I3,' symmetric matrix ',A/)
  102 FORMAT(1X,I4,5X,10G12.5/(10X,10G12.5/))
      END

      SUBROUTINE SMPRV(X,V,N)
*     Prints the n-vector X and the symmetric N-by-N  covariance  matrix
*     V, the latter as a correlation matrix.
*
*                    - - -
*        CALL SMPRTC(X,V,N)     X = vector of parameters
*                               V = covariance matrix
*
      REAL X(*),V(*)
      INTEGER MC(15)
*
      WRITE(*,103)
      WRITE(*,101)
      II=0
      DO I=1,N
       IJ=II
       II=II+I
       ERR=0.0
       IF(V(II).GT.0.0) ERR=SQRT(V(II))
       L=0
       JJ=0
       DO J=1,I
        JJ=JJ+J
        IJ=IJ+1
        RHO=0.0
        PD=V(II)*V(JJ)
        IF(PD.GT.0.0) RHO=V(IJ)/SQRT(PD)
        L=L+1
        MC(L)=100.0*ABS(RHO)+0.5
        IF(RHO.LT.0.0) MC(L)=-MC(L)
        IF(J.EQ.I.OR.L.EQ.15) THEN
           IF(J.LE.15) THEN
              IF(J.EQ.I) THEN
                 WRITE(*,102) I,X(I),ERR,(MC(M),M=1,L-1)
              ELSE
                 WRITE(*,102) I,X(I),ERR,(MC(M),M=1,L)
              END IF
           ELSE
              IF(J.EQ.I) THEN
                 WRITE(*,103) (MC(M),M=1,L-1)
              ELSE
                 WRITE(*,103) (MC(M),M=1,L)
              END IF
              L=0
           END IF
        END IF
       END DO
      END DO
      WRITE(*,104)
  100 RETURN
  101 FORMAT(9X,'Param',7X,'error',7X,'correlation coefficients'/)
  102 FORMAT(1X,I5,2G12.4,1X,15I5)
  103 FORMAT(31X,15I5)
  104 FORMAT(33X,'(correlation coefficients in percent)')
      END

      SUBROUTINE DSMPRV(X,V,N)
*     Prints the n-vector X and the symmetric N-by-N  covariance  matrix
*     V, the latter as a correlation matrix.
*
*                    - - -
*        CALL DPRPAC(X,V,N)     X = vector of parameters
*                               V = covariance matrix
*
      DOUBLE PRECISION X,V
      DIMENSION X(*),V(*)
      INTEGER MC(15)
*
      WRITE(*,103)
      WRITE(*,101)
      II=0
      DO I=1,N
       IJ=II
       II=II+I
       ERR=0.0
       IF(V(II).GT.0.0) ERR=SQRT(V(II))
       L=0
       JJ=0
       DO J=1,I
        JJ=JJ+J
        IJ=IJ+1
        RHO=0.0
        PD=V(II)*V(JJ)
        IF(PD.GT.0.0) RHO=V(IJ)/SQRT(PD)
        L=L+1
        MC(L)=100.0*ABS(RHO)+0.5
        IF(RHO.LT.0.0) MC(L)=-MC(L)
        IF(J.EQ.I.OR.L.EQ.15) THEN
           IF(J.LE.15) THEN
              IF(J.EQ.I) THEN
                 WRITE(*,102) I,X(I),ERR,(MC(M),M=1,L-1)
              ELSE
                 WRITE(*,102) I,X(I),ERR,(MC(M),M=1,L)
              END IF
           ELSE
              IF(J.EQ.I) THEN
                 WRITE(*,103) (MC(M),M=1,L-1)
              ELSE
                 WRITE(*,103) (MC(M),M=1,L)
              END IF
              L=0
           END IF
        END IF
       END DO
      END DO
      WRITE(*,104)
  100 RETURN
  101 FORMAT(9X,'Param',7X,'error',7X,'correlation coefficients'/)
  102 FORMAT(1X,I5,2G12.4,1X,15I5)
  103 FORMAT(31X,15I5)
  104 FORMAT(33X,'(correlation coefficients in percent)')
      END

      SUBROUTINE SMTOS(V,I,W,J,N)
*     Copy symmetric N-by-N matrix or a N-by-N submatrix of a  symmetric
*     matrix to another symmetric matrix
*
*                     - -     -
*          CALL SMTOS(V,I,W,J,N)
*                         - -
*
*     N rows and columns of the matrix V, starting from diagonal element
*     (i,i), are copied to the matrix W, starting  in  diagonal  element
*     (j,j). thus if a complete symmetric matrix has to be copied, i=j=1
*     has to be used.
*    
      REAL V(*),W(*)
*     ...
      IM=(I*I+I)/2-1
      JM=(J*J+J)/2-1
      DO K=1,N
       DO L=1,K
        IM=IM+1
        JM=JM+1
        W(JM)=V(IM)
       END DO
       IM=IM+I-1
       JM=JM+J-1
      END DO 
      END

      SUBROUTINE SMQU(V,VQ,N)
*     calculate square of symmetric matrix V:   VQ:= V * V
*
      REAL V(1),VQ(1)
      DOUBLE PRECISION SUM
*     ...
      IJ=0
      DO I=1,N
       DO J=1,I
        IJ=IJ+1
        SUM=0.0
        IK=(I*I-I)/2
        JK=(J*J-J)/2
        DO K=1,N
         IF(K.LE.I) IK=IK+1
         IF(K.LE.J) JK=JK+1
         SUM=SUM+V(IK)*V(JK)
         IF(K.GE.I) IK=IK+K
         IF(K.GE.J) JK=JK+K
        END DO
        VQ(IJ)=SUM
       END DO
      END DO
      END

      SUBROUTINE SMVA(V,A,B,N,M)
*     multiply symmetric N-by-N matrix and general N-by-M matrix
*
*                  - -   - -
*        CALL SMVA(V,A,B,N,M)
*                      -
*
*         B   :=    V   *    A
*        N*M       N*N      N*M
*
      REAL V(1),A(1),B(1)
      DOUBLE PRECISION SUM
      IJS=1
      IKB=0
      DO I=1,N
       DO K=1,M
        SUM=0.0
        IJ=IJS
        JK=K
        DO J=1,N
         SUM=SUM+V(IJ)*A(JK)
         IF(J.LT.I) THEN
            IJ=IJ+1
         ELSE 
            IJ=IJ+J
         END IF 
         JK=JK+M
        END DO 
        IKB=IKB+1
        B(IKB)=SUM
       END DO 
       IJS=IJS+I
      END DO 
      END

      SUBROUTINE DSMVA(V,A,B,N,M)
*     multiply symmetric N-by-N matrix and general N-by-M matrix
*
*                   - -   - -
*        CALL DSMVA(V,A,B,N,M)
*                       -
*
*         B   :=    V   *    A
*        N*M       N*N      N*M
*
      DOUBLE PRECISION V,A,B
      DIMENSION V(*),A(*),B(*)
      DOUBLE PRECISION SUM
      IJS=1
      IKB=0
      DO I=1,N
       DO K=1,M
        SUM=0.0
        IJ=IJS
        JK=K
        DO J=1,N
         SUM=SUM+V(IJ)*A(JK)
         IF(J.LT.I) THEN
            IJ=IJ+1
         ELSE 
            IJ=IJ+J
         END IF 
         JK=JK+M
        END DO 
        IKB=IKB+1
        B(IKB)=SUM
       END DO 
       IJS=IJS+I
      END DO 
      END

      SUBROUTINE SOLVA(V,R,N,U,D)
*     solve ill-conditioned linear system of equations V X = R
*     
*     CALL SOLVA(V,R,N,U,D)
*
*     define positive limit DLIM for eigenvalues
*     DO I=1,N
*      IF(D(I).GT.DLIM) THEN
*         D(I)=1.0/D(I)
*      ELSE
*         D(I)=0.0
*      END IF
*     END DO
*     CALL SOLVB(V,R,N,U,D) 
*
      REAL V(*),R(*),U(N,N),D(*)
*     ...
*     transformation to diagonal matrix ... 
      CALL SMHHL(V,N,D,U,-1)   
*     with eigenvalues D in descending order
      END 

      SUBROUTINE SOLVB(V,R,N,U,D,DLIM)
*     ... see SOLVA 
      REAL V(*),R(*),U(N,N),D(*)
      DOUBLE PRECISION SUM 
*     ...   
      DO I=1,N
       IF(D(I).GT.ABS(DLIM)) THEN
          D(I)=1.0/D(I)
       ELSE
          D(I)=0.0
       END IF
      END DO 
*     inverse matrix V by product UT D U      
      IK=0 
      DO I=1,N
       DO K=1,I
        SUM=0.0
        DO J=1,N
         SUM=SUM+U(I,J)*D(J)*U(K,J)
        END DO
        IK=IK+1
        V(IK)=SUM
       END DO
      END DO 
*     copy R to D        
      DO I=1,N
       D(I)=R(I)
      END DO 
*     solution by product V R => R
      II=0 
      DO I=1,N
       SUM=0.0
       IJ=II
       DO J=1,N
        IF(J.LE.I) IJ=IJ+1 
        SUM=SUM+D(J)*V(IJ)
        IF(J.GE.I) IJ=IJ+J
       END DO
       R(I)=SUM
       II=II+I
      END DO
      END 

      SUBROUTINE DSOLVA(V,R,N,U,D)
*     solve ill-conditioned linear system of equations V X = R
*     
*     CALL DSOLVA(V,R,N,U,D)
*
*     define positive limit DLIM for eigenvalues
*     CALL DSOLVB(V,R,N,U,D,DLIM) 
*
      DOUBLE PRECISION V(*),R(*),U(N,N),D(*)
*     ...
*     transformation to diagonal matrix ... 
      CALL DSMHHL(V,N,D,U,-1)   
*     with eigenvalues D in descending order
      END 

      SUBROUTINE DSOLVB(V,R,N,U,D,DLIM)
*     .. see DSOLVA 
      DOUBLE PRECISION V(*),R(*),U(N,N),D(*),DLIM
      DOUBLE PRECISION SUM 
*     ...   
      DO I=1,N
       IF(D(I).GT.ABS(DLIM)) THEN
          D(I)=1.0/D(I)
       ELSE
          D(I)=0.0
       END IF
      END DO 
*     inverse matrix V by product UT D U      
      IK=0 
      DO I=1,N
       DO K=1,I
        SUM=0.0
        DO J=1,N
         SUM=SUM+U(I,J)*D(J)*U(K,J)
        END DO
        IK=IK+1
        V(IK)=SUM
       END DO
      END DO 
*     copy R to D        
      DO I=1,N
       D(I)=R(I)
      END DO 
*     solution by product V R => R
      II=0 
      DO I=1,N
       SUM=0.0
       IJ=II
       DO J=1,N
        IF(J.LE.I) IJ=IJ+1 
        SUM=SUM+D(J)*V(IJ)
        IF(J.GE.I) IJ=IJ+J
       END DO
       R(I)=SUM
       II=II+I
      END DO
      END 

      SUBROUTINE TQHHL(N,D,E,U,IDIR)
*     Auxiliary routine for SMHHL
*
      REAL D(1),E(1),U(N,N)
      DATA EPS/1.0E-6/
*     ...
      IF(IDIR.GE.0) THEN
*        increasing order of eigenvalues
         IREV=1
      ELSE
*        decreasing order of eigenvalues
         IREV=0
      END IF
      DO I=2,N
       E(I-1)=E(I)
      END DO
      E(N)=0.0
      B=0.0
      F=0.0
*
      DO 60 L=1,N
      J=0
      H=EPS*(ABS(D(L))+ABS(E(L)))
      IF(B.LT.H) B=H
*
*     look for small sub-diagonal element
*
      DO M=L,N
       IF(ABS(E(M)).LE.B) GOTO 25
      END DO
      M=L
   25 IF(M.EQ.L) GOTO 59
*
*     next iteration
*
   30 IF(J.EQ.30) GOTO 100
      J=J+1
      G=D(L)
      P=(D(L+1)-G)/(2.0*E(L))
      R=SQRT(1.0+P*P)
      D(L)=E(L)
      IF(P.LT.0.0) D(L)=D(L)/(P-R)
      IF(P.GE.0.0) D(L)=D(L)/(P+R)
      H=G-D(L)
      DO I=L+1,N
       D(I)=D(I)-H
      END DO
      F=F+H
*     QL transformation
      P=D(M)
      C=1.0
      S=0.0
      M1=M-1
      DO 50 II=L,M1
      I=M1+L-II
      G=C*E(I)
      H=C*P
      IF(ABS(P).GE.ABS(E(I))) THEN
         C=E(I)/P
         R=SQRT(1.0+C*C)
         E(I+1)=S*P*R
         S=C/R
         C=1.0/R
      ELSE
         C=P/E(I)
         R=SQRT(1.0+C*C)
         E(I+1)=S*E(I)*R
         S=1.0/R
         C=C/R
      END IF
      P=C*D(I)-S*G
      D(I+1)=H+S*(C*G+S*D(I))
*     form vector
      DO K=1,N
       H=U(K,I+1)
       U(K,I+1)=S*U(K,I)+C*H
       U(K,I)=C*U(K,I)-S*H
      END DO
   50 CONTINUE
      E(L)=S*P
      D(L)=C*P
      IF(ABS(E(L)).GT.B) GOTO 30
   59 D(L)=D(L)+F
   60 CONTINUE
*
      DO 90 I=1,N
      K=I
      P=D(I)
*     search for largest/smallest eigenvalue
      DO 70 J=I,N
      IF(IREV.EQ.0.AND.D(J).LE.P) GOTO 70
      IF(IREV.EQ.1.AND.D(J).GE.P) GOTO 70
      K=J
      P=D(J)
   70 CONTINUE
      IF(K.NE.I) THEN
*        exchange
         D(K)=D(I)
         D(I)=P
         DO J=1,N
          P=U(J,I)
          U(J,I)=U(J,K)
          U(J,K)=P
         END DO
      END IF
   90 CONTINUE
  100 RETURN
      END

      SUBROUTINE DTQHHL(N,D,E,U,IDIR)
*     Auxiliary routine for SMHHL
*
      DOUBLE PRECISION D(*),E(*),U(N,N),EPS
      DATA EPS/1.0D-6/
*     ...
      IF(IDIR.GE.0) THEN
*        increasing order of eigenvalues
         IREV=1
      ELSE
*        decreasing order of eigenvalues
         IREV=0
      END IF
      DO I=2,N
       E(I-1)=E(I)
      END DO
      E(N)=0.0
      B=0.0
      F=0.0
*
      DO 60 L=1,N
      J=0
      H=EPS*(ABS(D(L))+ABS(E(L)))
      IF(B.LT.H) B=H
*
*     look for small sub-diagonal element
*
      DO M=L,N
       IF(ABS(E(M)).LE.B) GOTO 25
      END DO
      M=L
   25 IF(M.EQ.L) GOTO 59
*
*     next iteration
*
   30 IF(J.EQ.30) GOTO 100
      J=J+1
      G=D(L)
      P=(D(L+1)-G)/(2.0*E(L))
      R=SQRT(1.0+P*P)
      D(L)=E(L)
      IF(P.LT.0.0) D(L)=D(L)/(P-R)
      IF(P.GE.0.0) D(L)=D(L)/(P+R)
      H=G-D(L)
      DO I=L+1,N
       D(I)=D(I)-H
      END DO
      F=F+H
*     QL transformation
      P=D(M)
      C=1.0
      S=0.0
      M1=M-1
      DO 50 II=L,M1
      I=M1+L-II
      G=C*E(I)
      H=C*P
      IF(ABS(P).GE.ABS(E(I))) THEN
         C=E(I)/P
         R=SQRT(1.0+C*C)
         E(I+1)=S*P*R
         S=C/R
         C=1.0/R
      ELSE
         C=P/E(I)
         R=SQRT(1.0+C*C)
         E(I+1)=S*E(I)*R
         S=1.0/R
         C=C/R
      END IF
      P=C*D(I)-S*G
      D(I+1)=H+S*(C*G+S*D(I))
*     form vector
      DO K=1,N
       H=U(K,I+1)
       U(K,I+1)=S*U(K,I)+C*H
       U(K,I)=C*U(K,I)-S*H
      END DO
   50 CONTINUE
      E(L)=S*P
      D(L)=C*P
      IF(ABS(E(L)).GT.B) GOTO 30
   59 D(L)=D(L)+F
   60 CONTINUE
*
      DO 90 I=1,N
      K=I
      P=D(I)
*     search for largest/smallest eigenvalue
      DO 70 J=I,N
      IF(IREV.EQ.0.AND.D(J).LE.P) GOTO 70
      IF(IREV.EQ.1.AND.D(J).GE.P) GOTO 70
      K=J
      P=D(J)
   70 CONTINUE
      IF(K.NE.I) THEN
*        exchange
         D(K)=D(I)
         D(I)=P
         DO J=1,N
          P=U(J,I)
          U(J,I)=U(J,K)
          U(J,K)=P
         END DO
      END IF
   90 CONTINUE
  100 RETURN
      END

      FUNCTION VXY(X,Y,N)
*     Scalar product of two vectors
*                - - -                 T
*        S = VXY(X,Y,N)           S = X  * Y (scalar product)
*
      REAL X(*),Y(*)
      DOUBLE PRECISION SUM
*     ...
      SUM=0.0D0
      M=N/4
      I=0
      DO J=1,M
       SUM=SUM+X(I+1)*Y(I+1)+X(I+2)*Y(I+2)+X(I+3)*Y(I+3)+X(I+4)*Y(I+4)
       I=I+4
      END DO
      DO J=I+1,N
       SUM=SUM+X(J)*Y(J)
      END DO
      VXY=SUM
      END