SUBROUTINE PROVAR(N,P,LDP,SIGMA,LDSIG)
C  This routine constructs the variance matrix, from a multivariate projection
C  matrix (for description of latter see VARPRO).  This routine is the inverse
C  of VARPRO.
C  Inputs:
C  N         Order of multivariate distribution (= number of variables)
C  P         Projection matrix
C  LDP       Leading dimension of P
C  LDSIG     Leading dimension of SIGMA
C  Outputs:
C  SIGMA     Variance matrix
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION P(LDP,N),SIGMA(LDSIG,N)
      CALL UM0SET(N,N,A1,LD)
      DO 20 I=1,N
      AA=1.0D0/P(I,I)
      A1(I,I)=AA
      IF(I.EQ.1)GO TO 20
      DO 10 J=1,I-1
   10 A1(I,J)=-AA*P(I,J)
   20 CONTINUE
      CALL DLINRT(N,A1,LD,1,A1,LD)
      DO 40 I=1,N
      DO 40 J=1,I
      AA=0.0D0
      DO 30 K=1,J
   30 AA=AA+A1(I,K)*A1(J,K)
      SIGMA(I,J)=AA
   40 SIGMA(J,I)=AA
      RETURN
      END
      SUBROUTINE VARPRO(N,SIGMA,LDSIG,P,LDP)
C  This routine constructs a recursive projection matrix P, given a variance
C  matrix.  The linear projection of x(i) on x(1), ... x(i-1) is
C    x(i)=p(i,1)x(1)+...+p(i,i-1)x(i-1)+u(i); s.d.[u(i)]=p(i,i).
C  Inputs:
C  N         Order of multivariate distribution (= number of variables)
C  SIGMA     Variance matrix
C  LDSIG     Leading dimension of SIGMA
C  LDP       Leading dimension of P
C  Outputs:
C  P         Projection matrix
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION SIGMA(LDSIG,N),P(LDP,N)
      CALL DLFTDS(N,SIGMA,LDSIG,A1,LD)
      CALL DLINRT(N,A1,LD,2,A1,LD)
      CALL UM0SET(N,N,P,LDP)
      DO 20 I=1,N
      AA=1.0D0/A1(I,I)
      P(I,I)=AA
      IF(I.EQ.1)GO TO 20
      DO 10 J=1,I-1
   10 P(I,J)=-AA*A1(J,I)
   20 CONTINUE
      RETURN
      END