subroutine stationary(n,a,sig, phi)
!
!  STATIONARY 
!     maps a square matrix to one with eigenvalues inside the unit 
!     circle; the inverse map is in istationary.f90.
!
!     stationary(n,A,Sig, Phi) takes the matrix A, which is a general
!     square matrix, and Sig, which is a symmetric positive definite 
!     matrix and computes:
!
!           P   = inv( chol( I+A*A' )' ) * A
!           B   = inv( chol( I-P*P' )' )
!           T   = chol(Sig)'* B
!           Phi = T* P *inv(T)
!
!     Reference: C.F. Ansley and R. Kohn, "A Note on Reparameterizing
!        a Vector Autoregressive Moving Average Model to Enforce
!        Staionarity," J. of Statistical Computational Simulations,
!        1986, Vol 24, 99-106.
!                                                                               
!     Ellen McGrattan, 5-13-06
!     Revised, 5-30-06

  implicit none
  integer, parameter                 :: lwork=5000
  integer                            :: i,j,info
  integer, intent(in)                :: n
  integer, dimension(lwork)          :: ipiv
  real, dimension(n,n), intent(in)   :: a,sig
  real, dimension(n,n), intent(out)  :: phi
  real, dimension(n,n)               :: p,b,t
  real, dimension(lwork)             :: work

  !
  ! Compute P, which has singular values less <= 1
  !
  p  = 0.
  do i = 1,n
    p(i,i) = 1.
  enddo
  b   = p
  p   = p + matmul(a,transpose(a))
  call dpotrf('L',n,p,n,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dpotrf'
    write(*,*) 'info = ',info
    return
  endif
  do j=2,n
    do i=1,j-1
      p(i,j) = 0.
    enddo
  enddo
  call dgetrf(n,n,p,n,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dgetrf'
    write(*,*) 'info = ',info
    return
  endif 
  call dgetri(n,p,n,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dgetri'
    write(*,*) 'info = ',info
    return
  endif
  p   = matmul(p,a)

  !
  ! Compute B=inv(chol(I-P*P')')
  !
  b   = b - matmul(p,transpose(p))
  call dpotrf('L',n,b,n,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with second dpotrf'
    write(*,*) 'info = ',info
    return
  endif
  do j=2,n
    do i=1,j-1
      b(i,j) = 0.
    enddo
  enddo
  call dgetrf(n,n,b,n,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with second dgetrf'
    write(*,*) 'info = ',info
    return
  endif 
  call dgetri(n,b,n,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with second dgetri'
    write(*,*) 'info = ',info
    return
  endif

  !
  ! Compute T=chol(Sig)'* B
  ! 
  t  = sig
  call dpotrf('L',n,t,n,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with third dpotrf'
    write(*,*) 'info = ',info
    return
  endif
  do j=2,n
    do i=1,j-1
      t(i,j) = 0.
    enddo
  enddo
  t   = matmul(t,b)
  
  !
  ! Finally, compute Phi = T*P*inv(T)
  !
  phi = matmul(t,p)
  call dgetrf(n,n,t,n,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with third dgetrf'
    write(*,*) 'info = ',info
    return
  endif 
  call dgetri(n,t,n,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with third dgetri'
    write(*,*) 'info = ',info
    return
  endif
  phi = matmul(phi,t)

end subroutine stationary