subroutine liapqr(nlags,ldder,der,n,lexpsv,lexpqr)

c Purpose: create the Jacobian and update the product
c       of Jacobians with the current vector of partial 
c       derivatives; using this product the
c       Lyapunov exponents are estimated
c
c On Entry:
c   nlags                number rows in der, number of lags in model
c   der(nlags,n)         matrix of partial derivatives
c   ldder                leading dimension (#rows) of der in calling routine
c   
c On Exit:
c   lexpsv(nlags)          arrays of ordered estimated Lyapunov exponents, from
c   lexpqr(nlags)             SVD and QR decompositions, respectively
c Subprograms Called Directly:
c       Linpack - dsvdc,dqrdc
c Should be linked with Double-precision BLAS

      parameter(LMAX=65,LDW=6*LMAX)
      INTEGER nlags,info

      double precision jac(lmax,lmax),s(lmax,lmax),t(lmax,lmax),
     * v(lmax,lmax),u(lmax,lmax),e(LMAX),epsln,
     * der(ldder,n),lexpsv(nlags),cnorm,w(LMAX),norm,work(LDW),
     * lexpqr(nlags),qraux(LMAX),one,zero
      integer jpvt(LMAX)
	parameter(one=1.d0, zero=0.d0)

      epsln= 1.0d-15
      if (LMAX.lt.ldder) then
          write(*,*) ' Not enough storage in liapest'
          stop
      endif

c---------------- Initialize S=identity matrix
      do 110 i=1,nlags
        do 100 j=1,nlags
            s(i,j)=zero
100     continue
        s(i,i)=one
110   continue


c----------------- Initialize fixed elements of Jacobian matrix
      DO 130 I=1,NLAGS
          DO 120 J=1,NLAGS                                          
              jac(I,J)=zero                                      
              IF (I.EQ.(J+1)) THEN                               
                 jac(I,J)=one                                     
              ENDIF                                            
120            CONTINUE                                        
130        CONTINUE                                           
      CNORM=zero                                            

c--------------------- Loop over kk=1,2,.. n to compute T_n
      do 200 kk=1,n

c--------------------- Put derivatives at time kk into Jac
      DO 10 I=1,NLAGS                                              
          jac(1,I)=DER(I,kk)                                           
10    CONTINUE

c--------------------- Update T(kk-1) by multiplying by Jac(kk)
      DO 13 I=1,NLAGS
          DO 12 J=1,NLAGS                                     
              T(I,J)=zero                                    
              DO 11 K=1,NLAGS                               
                  T(I,J)=jac(I,K)*S(K,J)+T(I,J)              
11            CONTINUE
12        CONTINUE
13    CONTINUE


c------------------ Normalize to Max(abs(T_i,j)) == 1
      NORM=zero
      DO 15 I=1,NLAGS                                             
          DO 14 J=1,NLAGS                                        
              NORM=MAX(NORM,ABS(T(I,J)))                        
14            CONTINUE                                         
15        CONTINUE                                            

      DO 17 I=1,NLAGS                                                 
          DO 16 J=1,NLAGS                                            
              T(I,J)=T(I,J)/NORM
              S(I,J)=T(I,J)
16            CONTINUE                                             
17        CONTINUE                                                

c---------------- Update Cnorm=cumulative log of renormalizations
      CNORM=LOG(NORM)+CNORM                                     

200   continue

c------------- Estimate of LE's based on SV decomposition
c     SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
      call dsvdc(t,LMAX,nlags,nlags,W,E,U,LMAX,V,LMAX,         
     *           work,00,info)

      if( info.gt.0) then
           write(*,*) 'info from svddc', info
      endif

      do 210 k=1, nlags
        if( abs(w(k)).gt.epsln) then
           lexpsv(k)= (LOG(ABS(w(k)))+CNORM)/dfloat(n)
        else
           lexpsv(k)=(log(epsln) +cnorm)/dfloat(n)
       endif
210   continue

c---------- Estimate of LE's based on QR decomposition
c     SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)
      call dqrdc(S,LMAX,NLAGS,NLAGS,QRAUX,JPVT,WORK,0)

      do 220 k=1,nlags
        if(abs(S(k,k)).gt.epsln) then
           lexpqr(k)= (LOG(ABS(S(k,k)))+CNORM)/dfloat(n)
        else
           lexpqr(k)= (log(epsln)+cnorm)/dfloat(n)
        endif
220   continue

      if(nlags.gt.1) then

c    -- sort lexpqr in ascending order
	call sort(nlags,lexpqr)

c    -- reverse into descending order
        do 230 k=1,nlags
          w(k)=lexpqr(nlags+1-k)
230     continue
        do 240 k=1,nlags
          lexpqr(k)=w(k)
240     continue
      endif
      RETURN                                                          
      END


      subroutine sort(n,x)
      implicit real*8(a-h,o-z)
      dimension x(n)
      l=n/2+1
      ix=n
5     continue
        if(l.gt.1)then
          l=l-1
          xx=x(l)
        else
          xx=x(ix)
          x(ix)=x(1)
          ix=ix-1
          if(ix.eq.1)then
            x(1)=xx
            return
          endif
        endif
        i=l
        j=l+l
10      if(j.le.ix)then
          if(j.lt.ix)then
            if(x(j).lt.x(j+1))j=j+1
          endif
          if(xx.lt.x(j))then
            x(i)=x(j)
            i=j
            j=j+j
          else
            j=ix+1
          endif
        go to 10
        endif
        x(i)=xx
      go to 5
      end