subroutine istationary(n,phi,sig, a)
!
!  ISTATIONARY 
!     is the inverse map of stationary.m which maps a general square
!     matrix to one with eigenvalues inside the unit circle.
!
!     istationary(n,Phi,Sig, A) takes the matrix Phi, which has
!     eigenvalues inside the unit circle, and Sig, which is a
!     symmetric positive definite matrix and computes:
!
!           T = chol(X)',  X = Phi*X*Phi' + Sig
!           P = inv(T)* Phi *T
!           B = inv( chol( I-P*P' )' )
!           A = B * P
!
!     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)          :: iwork,ipiv
  real, dimension(n,n), intent(in)   :: phi,sig
  real, dimension(n,n), intent(out)  :: a
  real, dimension(n,n)               :: p,pt,x,tem
  real, dimension(lwork)             :: work

  !
  ! Solve for X in X=Phi*X*Phi'+Sig
  !
  p    = phi
  pt   = -transpose(phi)
  x    = sig
  call sb04qd(n,n,p,n,pt,n,x,n,tem,n,iwork,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY.f90: problem with sb04qd'
    write(*,*) 'info = ',info
    return
  endif
  !
  ! Compute the Cholesky decomposition of the solution 
  !
  call dpotrf('L',n,x,n,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY: problem with dpotrf'
    write(*,*) 'info = ',info
    return
  endif
  do j=2,n
    do i=1,j-1
      x(i,j) = 0.
    enddo
  enddo
  p = x
  !
  ! Invert it
  !
  call dgetrf(n,n,x,n,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY: problem with dgetrf'
    write(*,*) 'info = ',info
    return
  endif 
  call dgetri(n,x,n,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY: problem with dgetri'
    write(*,*) 'info = ',info
    return
  endif
  !
  ! And use it to compute P = inv(T)* Phi *T which has singular
  ! values less than 1
  !
  
  p =  matmul(phi,p)
  p =  matmul(x,p)
  !
  ! Next, map P to the general matrix A = inv(chol(I-P*P')')*P
  !

  a = 0.
  do i = 1,n
    a(i,i) = 1.
  enddo
  a  = a - matmul(p,transpose(p))
  call dpotrf('L',n,a,n,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY: problem with second dpotrf'
    write(*,*) 'info = ',info
    return
  endif
  do j=2,n
    do i=1,j-1
      a(i,j) = 0.
    enddo
  enddo
  call dgetrf(n,n,a,n,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY: problem with second dgetrf'
    write(*,*) 'info = ',info
    return
  endif 
  call dgetri(n,a,n,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in ISTATIONARY: problem with second dgetri'
    write(*,*) 'info = ',info
    return
  endif
  a  = matmul(a,p)

end subroutine istationary