c------------------------------------------------------------------
c------------------------------------------------------------------
C
C     GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified 13 November 1997.
C
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     The GURU package consists of the following routines:
C
C     GURC2  - generates matrix C;
C     GURINV - matrix inversion;
C     GURMAT - generates response matrix from Blobel's example;
C     GURGAU - generates a pure gaussian response matrix;
C     GURTST - generates test distribution;
C     GURDAT - generates data distribution;
C     GURROT - rescaling, rotation and other preliminary operations;
C     GURTAN - estimates the effective rank of the system;
C     GURUNF - actual unfolding;
C     GURCHI - calculates chi2's.
C
C     The GURU package also uses:
C
C     SVD    - Singular Value Decomposition (from CERNLIB);
C     ROTATE - elementary rotations - used by SVD;
C     RANNOR - normally distributed random number generator.
C
c------------------------------------------------------------------
c------------------------------------------------------------------
C
C      subroutine rannor(x,y)
C
C     returns two (roughly) normally distributed random numbers x,y
C     with mean=0.0 and sigma=1.0. RANDOM is a uniformly distributed
C     random number in the interval [0,1).      
C
C      x=-6.0
C      y=-6.0
C      do 100 i=1,12
C      x=x+random(irs)
C 100  y=y+random(irs)
C
C      return
C      end
C
C---------------------------------------------------------------------

      real function tauint(nx,s,k)
C
C     returns the value of tau corresponding to an integer number
C     of effective equations k.
C
C     ff=SUM_i=1^nx{ s(i)**2/(s(i)**2+tau) }
C
C     solving equation ff-k=0 usig Newton's method.
C


      dimension s(nx)
      integer nx,k
      real t,ff,fp,s2,s3,s4

      parameter(ggt=1.e9)
      parameter(eps=1.e-4)

      tauint=-1.0
      if(k.lt.1.or.k.gt.nx) return

      tauint=s(nx)/ggt
      if(k.eq.nx) return

c      tauint=s(1)*ggt
c      if(k.eq.1) return


C     start of iterations

C     starting value
 
      if(k.eq.nx) then
        t=0.0
      else
        t=s(k+1)**2
      endif


 100  ff=-float(k)
      fp=0.0
      do 200 i=1,nx
        s2=s(i)**2
        s3=s2/(s2+t)
        s4=s3/(s2+t)
        ff=ff+s3                 ! summing up F
        fp=fp+s4                 ! summing up -F'
 200  continue

      delta=ff/fp

      t=t+delta

C      write(*,*) k, ff, fp, delta, t

      if(abs(delta).gt.eps*t) goto 100

      tauint=t

      return
      end

C---------------------------------------------------------------------

      subroutine gurc2(nx,ici,c)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Makes the "second derivative" matrix C(nx,nx).
C
C     Input:
C
C     nx      - dimension of the matrix;
C     ici     - 0 for length, 2 for second derivative.
C
C     Output:
C
C     c(nx,nx)- the matrix itself.
C
C     For ici=0, C=unit matrix. 
C
C     For ici=2 :
C     C(i,i)=-2.0 + xi; C(i,i+1)=C(i+1,i)=1.0; C(i,j)=0.0 otherwise.
C     The small parameter xi=0.001 makes the matrix non-singular.
C     Another parameter, endp=1.0, determines the behaviour near the
C     endpoints.
C
      dimension c(nx,nx)

      do 90 j=1,nx
        do 90 i=1,nx
         c(i,j)=0.
   90 continue

      if(ici.eq.0)then
        do 92, i=1,nx
          c(i,i)=1.0
   92   continue
      elseif(ici.eq.2)then
        xi=1.e-3
        endp=1.0
        do 100 j=1,nx
          c(j,j)=-2.
          if(j.eq.1.or.j.eq.nx)c(j,j)=c(j,j)+endp
          if(j.ne.nx)c(j,j+1)=1.
          if(j.ne.nx)c(j+1,j)=1.
  100   continue
        do 180 i=1,nx
  180   c(i,i)=c(i,i)+xi
      endif

      return
      end

c---------------------------------------------------------------------

      subroutine gurinv(n,u,s,v,nep,p)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Calculates the (pseudo)inverse of a square n*n matrix u,
C     using its Singular Value Decomposition (subroutine SVD).
C     The input matrix is decomposed into the product u*s*v^T,
C     where u and v are orthogonal matrices, and s is diagonal,
C     with non-negative diagonal elements called singular values
C     (^T denotes the transpose matrix). The inverse matrix p
C     is then calculated as p=v*s^-1*u^T, where s-1 is a diagonal 
C     matrix with diagonal elements (s^-1)_ii=(s_ii)^-1. If and when
C     any of the singular values is smaller than the parameter eps, the
C     matrix is considered to be too singular to be inverted properly.
C     In these cases the small singular values are discarded and the
C     resulting matrix p is called effective Moore-Penrose pseudoinverse.
C     Parameter nep counts the singular values which are larger than eps.
C     If on output nep=n, than p is the exact inverse (up to machine
C     accuracy), which satisfies the following conditions:
C
C       u*p*u=u;    p*u*p=p;   p*u and u*p are symmetric.  
C
C     If nep<n, p is the effective pseudoinverse, for which the first
C     condition is replaced by
C
C       |u*p*u-u|<eps.  
C
C     Input:
C
C     n      - dimension of the matrices/vectors involved;
C     u(n,n) - the matrix to be inverted.
C
C     Output:
C
C     u(n,n) - the left orthogonal matrix of the SVD;
C     s(n)   - the vector of singular values, of u;
C     v(n,n) - the right orthogonal matrix of the SVD;
C     nep    - number of significant singular values;
C     p(n,n) - the (pseudo)inverse.
C
      dimension u(n,n),s(n),v(n,n),p(n,n)

c      eps=1.e-6
c      eps=1.0
c      eps=1.0E6

      call svd(u,s,v,n,n,n,n,0,.true.,.true.)

      eps=1.e-6*s(1)

      nep=0
      do 20 i=1,n
        if(s(i).gt.eps)nep=nep+1
        do 20 j=1,n
          t=0.0
          do 10 k=1,n
            t=t+v(i,k)*u(j,k)*s(k)/(s(k)**2+eps**2)
   10     continue
          p(i,j)=t
   20 continue
      return
      end

c---------------------------------------------------------------------

      subroutine gurinc(n,u,s,v,nep,p)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Calculates the (pseudo)inverse of a square n*n matrix u,
C     using its Singular Value Decomposition (subroutine SVD).
C     The input matrix is decomposed into the product u*s*v^T,
C     where u and v are orthogonal matrices, and s is diagonal,
C     with non-negative diagonal elements called singular values
C     (^T denotes the transpose matrix). The inverse matrix p
C     is then calculated as p=v*s^-1*u^T, where s-1 is a diagonal 
C     matrix with diagonal elements (s^-1)_ii=(s_ii)^-1. If and when
C     any of the singular values is smaller than the parameter eps, the
C     matrix is considered to be too singular to be inverted properly.
C     In these cases the small singular values are discarded and the
C     resulting matrix p is called effective Moore-Penrose pseudoinverse.
C     Parameter nep counts the singular values which are larger than eps.
C     If on output nep=n, than p is the exact inverse (up to machine
C     accuracy), which satisfies the following conditions:
C
C       u*p*u=u;    p*u*p=p;   p*u and u*p are symmetric.  
C
C     If nep<n, p is the effective pseudoinverse, for which the first
C     condition is replaced by
C
C       |u*p*u-u|<eps.  
C
C     Input:
C
C     n      - dimension of the matrices/vectors involved;
C     u(n,n) - the matrix to be inverted.
C
C     Output:
C
C     u(n,n) - the left orthogonal matrix of the SVD;
C     s(n)   - the vector of singular values, of u;
C     v(n,n) - the right orthogonal matrix of the SVD;
C     nep    - number of significant singular values;
C     p(n,n) - the (pseudo)inverse.
C
      dimension u(n,n),s(n),v(n,n),p(n,n)

      eps=1.e-6

      call svd(u,s,v,n,n,n,n,0,.true.,.true.)

      nep=0
      do 20 i=1,n
        if(s(i).gt.eps)nep=nep+1
        do 20 j=1,n
          t=0.0
          do 10 k=1,n
            if(s(k).gt.eps)t=t+v(i,k)*u(j,k)/s(k)
   10     continue
          p(i,j)=t
   20 continue
      return
      end

c---------------------------------------------------------------------

      subroutine gurmat(nb,rbl,rbu,nx,rxl,rxu,
     &                   apro,atru,xini)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Builds sample response matrix from Blobel example, and also
C     corresponding true matrix and initial x-vector.
C
C     Input:
C
C     nb        - number of bins in measured histogram;
C     rbl,rbu   - lower and upper values of the measured variable;
C     nx        - number of bins in generated/unfolded histograms;
C     rxl,rxu   - lower and upper values of the unfolded variable.
C
C     Output:
C
C     apro(nb,nx) - probability response matrix, calculated acording
C                   Blobel formulas;
C     atru(nb,nx) - number-of-events response matrix;
C     xini(nx)    - distribution of unfolded variable used for generating
C                   the response matrix. 
C
C-----------------------------------------------------------------------  
C in apro(i,k),  i must be MEASURED,     k must be  TRUE (GENERATED)
C    atru(i,k): {b, right axis in PAW}     {x, left axis in PAW}
C-----------------------------------------------------------------------  

      dimension apro(nb,nx),atru(nb,nx),xini(nx)                

C-----bin widths:

      xw=(rxu-rxl)/real(nx) 
      bw=(rbu-rbl)/real(nb) 

C-----Probability response matrix from Blobel's example:

      do 120 j=1,nx
        xi   = rxl+xw*(real(j)-0.5)                
        anorm= (1.-0.5*(1.-xi)**2)
        amu  = xi-0.05*xi**2
        do 110 i=1,nb
          bi   = rbl+bw*(real(i)-0.5) 
          apro(i,j)=anorm*3.6*exp(-(bi-amu)**2/0.02)*bw
          if(apro(i,j).lt.1.0e-20)apro(i,j)=0.0
  110   continue
  115   format(10e8.2)
  120 continue

C-----Uniform initial distribution:

      do 180 j=1,nx
        xini(j)= 250.0
 180  continue

C-----True number-of-events response matrix:

      do 520 i=1,nb                        
        do 520 j=1,nx
          atru(i,j)=apro(i,j)*xini(j)
 520  continue

      return
      end
                        

c---------------------------------------------------------------------

      subroutine gurgau(nb,rbl,rbu,nx,rxl,rxu,
     &                   apro,atru,xini)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Builds sample response matrix from Blobel example, and also
C     corresponding true matrix and initial x-vector.
C
C     Input:
C
C     nb        - number of bins in measured histogram;
C     rbl,rbu   - lower and upper values of the measured variable;
C     nx        - number of bins in generated/unfolded histograms;
C     rxl,rxu   - lower and upper values of the unfolded variable.
C
C     Output:
C
C     apro(nb,nx) - probability response matrix for a pure gaussian
C                   reflected against lower bound;
C     atru(nb,nx) - equals apro(nb,nx);
C     xini(nx)    - filled with 1.00 for all i=1,nx.
C
C-----------------------------------------------------------------------  
C in apro(i,k),  i must be MEASURED,     k must be  TRUE (GENERATED)
C    atru(i,k): {b, right axis in PAW}     {x, left axis in PAW}
C-----------------------------------------------------------------------  

      dimension apro(nb,nx),atru(nb,nx),xini(nx)                
      parameter(ss=0.00085)

      anorm=1./sqrt(2.*3.14159)/ss

C-----bin widths:

      xw=(rxu-rxl)/real(nx) 
      bw=(rbu-rbl)/real(nb) 

C-----Probability response matrix from Blobel's example:

      do 120 j=1,nx
        xi   = rxl+xw*(real(j)-0.5)                
        do 110 i=1,nb
          bi   = rbl+bw*(real(i)-0.5) 
          zz   = (bi-xi)/ss
          cc=0.0
          if(ABS(zz).lt.10.0) cc   = anorm*exp(-zz*zz/2.)
          apro(i,j) = cc
          bi   = rbl+bw*(real(i+2*j-1)-0.5) 
          zz   = (bi-xi)/ss
          cc=0.0
          if(ABS(zz).lt.10.0) cc   = anorm*exp(-zz*zz/2.)
          apro(i,j) = apro(i,j) + cc
  110   continue
  120 continue

C-----Uniform initial distribution:

      do 180 j=1,nx
        xi   = rxl+xw*(real(j)-0.5)                
ccv        xini(j)= 1.0/xi/(xi*xi+(2.5/91.2/2)**2)
ccv        xini(j)= 1.0/(xi*xi+(0.002)**2)
        xini(j)= 1.0
 180  continue

C-----True number-of-events response matrix:

      do 520 i=1,nb                        
        do 520 j=1,nx
          atru(i,j)=apro(i,j)*xini(j)
 520  continue

      return
      end
                        

C-------------------------------------------------------------------------

      subroutine gurtst(nb,rbl,rbu,nx,rxl,rxu,apro,
     &                             xtst,btst,btster,xi)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified 13 November 1997.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Builds sample test and measured distributions from Blobel example.
C
C     Input:
C
C     nb          - number of bins in measured histogram;
C     rbl,rbu     - lower and upper values of the measured variable;
C     nx          - number of bins in generated/unfolded histograms;
C     rxl,rxu     - lower and upper values of the unfolded variable.
C     apro(nb,nx) - probability response matrix, as calculated by
C                   subroutine gurmat.
C
C     Output:
C
C     xtst(nx)    - calculated test distribution from Blobel example;
C     btst(nb)    - "measured" test distribution. The measurement
C                   process is simulated by the response matrix apro,
C                   so that btst=apro*xtst 
C     btster(nb)  - errors in btst. The errors are assumed to be statis-
C                   tical only and independent, so that btster(i)=
C                   sqrt(btst(i));
C     xi(nx)      - trial constant distrib. used in reg. term.
C

      dimension apro(nb,nx),xtst(nx),xi(nx)
      dimension btst(nb),btster(nb)

      evnum=51000                                   ! number of events 

      funt=0.                                       ! choice of distribution

C-----bin widths:

      xw=(rxu-rxl)/real(nx) 
      bw=(rbu-rbl)/real(nb) 

      do 100 i=1,nx
         xx=rxl+xw*(float(i)-0.5)
         xblob = (4.0/(4.  +(xx-0.4)**2)+                 
     +            0.4/(0.04+(xx-0.8)**2)+
     +            0.2/(0.04+(xx-1.5)**2))*0.1       ! Blobel example

         if(funt.eq.1) xblob=exp(-xx)               ! various distributions
         if(funt.eq.2) xblob=1.5+sin(4*xx)          ! to be unfolded
         if(funt.eq.3) xblob=exp((xx-2.)**2) 
         if(funt.eq.4) xblob=(xx-1.)**2 
         if(funt.eq.5) xblob=xx**0.5 

         xtst(i)=xblob*xw*evnum

 100  continue

      do 82 i=1,nx
        xi(i)=0.0
c        xi(i)=400*.25*(float(i)-4.-float(nx)/2.)**2
c        xi(i)=xtst(i)
  82  continue

      do 540 i=1,nb                      
         btst(i)=0.0
         do 530 k=1,nx
            btst(i)=btst(i)+apro(i,k)*xtst(k)
 530     continue
         if(btst(i).gt.0.)then
            btster(i)=sqrt((btst(i)))
         else
            btster(i)=0.0
         endif
 540  continue

      return
      end

C-------------------------------------------------------------------------

      subroutine gurdat(nb,btst,bdat,bdater)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Builds sample measured distribution bdat by smearing btst.
C
C     Input:
C
C     nb          - number of bins in measured histogram;
C     btst(nb)    - "measured" test distribution. The errors in btst
C                   are assumed to be statistical only and independent.
C
C     Output:
C
C     bdat(nb)    - a sample measured data distribution which is obtained
C                   from btst by adding to it random gaussian fluctuations 
C                   with standard deviations given by btster;
C     bdater(nb)  - errors in bdat. The errors are again assumed to be 
C                   statistical only and independent, so that bdater(i)=
C                   sqrt(bdat(i)).
C

      dimension btst(nb),bdat(nb),bdater(nb)



c-------------------real measured data should be input as bdat------------

      do 542 i=1,nb                          ! bdat = smeared btst
        call rannor(x,y)
        bdat(i)=btst(i)+x*sqrt(btst(i))
  542 continue

c------- errors in bdat are assumed statistical and independent ----------

      do 560 i=1,nb
         if(bdat(i).gt.0.)then
            bdater(i)=sqrt((bdat(i)))
         else
            bdat(i)=0.0
            bdater(i)=0.0
         endif
  560 continue

      return
      end

C----------------------------------------------------------

      subroutine gurrot(nb,nx,a,xini,bdat,bdater,c,cp,xi,
     &       atil,btil,uort,sing,vort,vreg,vrm1,xcmdm1,ddat)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified 13 November 1997.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Rescaling of equations, Singular Value Decomposition of the
C     response matrix, and system diagonalization. 
C
C     Input:
C
C     nb         - number of bins in the measured distribution;
C     nx         - number of bins in the unfolded distribution;
C     a(nb,nx)   - the response matrix;
C     xini(nx)   - initial x-distribution, used to generate the response
C                  matrix;
C     bdat(nb)   - measured distribution;
C     bdater(nb) - errors in the measured distribution;
C     c(nx,nx)   - the second derivative matrix;
C     cp(nx,nx)  - the inverse of matrix C;
C     xi(nx)      - constant part of reg. term.
C
C     Output:
C
C     atil(nb,nx)   - rescaled a, each row divided by bdater;
C     btil(nx)      - bdat divided by bdater;
C     uort(nb,nx)   - left orthogonal matrix of SVD;
C     sing(nx)      - vector of singular values of the matrix product A*C^-1;
C     vort(nb,nx)   - right orthogonal matrix of SVD;
C     vreg(nx,nx)   - matrix performing the rotation from the diagonalized
C                     system back to x. Prepared for subroutine gurunf;
C     vrm1(nx,nx)   - matrix performing the rotation to diagonalize the
C                     system; not used outside this subroutine;
C     xcmdm1(nx,nx) - inverse of the covariance matrix of the unfolded 
C                     distribution. 
C     ddat(nx)      - vector containing the r.h.s. of the diagonalized
C                     system, prepared for subroutines gurtan and gurunf;
C
      dimension a(nb,nx),xini(nx),bdat(nb),bdater(nb),xi(nx)
      dimension c(nx,nx),cp(nx,nx)
      dimension atil(nb,nx),btil(nb)
      dimension uort(nb,nx),sing(nx),vort(nx,nx),ddat(nx)
      dimension vreg(nx,nx),vrm1(nx,nx),xcmdm1(nx,nx)

      do 80 i=1, nb
        if (bdater(i).gt.0.)then
          btil(i)=bdat(i)/bdater(i)
          do 70 j=1,nx
  70      atil(i,j)=a(i,j)/bdater(i)
        else
          btil(i)=0.
          do 75 j=1,nx
  75      atil(i,j)=0.
        endif   
  80  continue

      do 82 i=1,nb
         do 82 j=1,nx
            if(xini(j).gt.0)then
               btil(i)=btil(i)-atil(i,j)*xi(j)/xini(j)
            endif
 82   continue

C      do 82 i=1,nb
C        do 82 j=1,nx
C          btil(i)=btil(i)-atil(i,j)*xi(j)/xini(j)
C  82  continue
      
      do 85 i=1,nb                   ! multiply Atil*Cp      
        do 85 j=1,nx
          uort(i,j)=0.0
          do 85 k=1,nx
            uort(i,j)=uort(i,j)+atil(i,k)*cp(k,j)
  85  continue

      call svd(uort,sing,vort,nb,nx,nb,nx,0,.true.,.true.)

      do 830 i=1,nx                          ! VR=xini*Cp*V 
        do 830 j=1,nx
          vreg(i,j)=0.0
          do 825 k=1,nx
            vreg(i,j)=vreg(i,j)+cp(i,k)*vort(k,j)
  825     continue
          vreg(i,j)=xini(i)*vreg(i,j)
  830 continue

      do 832 i=1,nx                          ! VRm1=VT*C/xini 
        do 832 j=1,nx
          vrm1(i,j)=0.0
          do 832 k=1,nx
            vrm1(i,j)=vrm1(i,j)+vort(k,i)*c(k,j)/xini(j)
  832 continue

      do 858 i=1,nx            ! X-1 = (1/xini*VT*C*S)T*(1/xini*VT*C*S)
        do 858 j=1,nx                        
          xcmdm1(i,j)=0.0
          do 858 k=1,nx
            xcmdm1(i,j)=xcmdm1(i,j)+vrm1(k,i)*vrm1(k,j)*sing(k)**2
  858 continue

      do 840 i=1,nx                          ! d = UT*b    
          ddat(i)=0.0
          do 835 k=1,nb
            ddat(i)=ddat(i)+uort(k,i)*btil(k) 
  835     continue
  840 continue

      return
      end

C----------------------------------------------------------

      subroutine guruni(nb,nx,atru,bdater,cp,uort,sing,vort,tau,
     &                        ainv,unitb,unitx)
C
C     Part of GURU package, version 0.93. Written by V.Kartvelishvili
C     in October 1995. Last modified on 27 January 1997.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Calculates the pseudo-inverse  AINV of the true response matrix
C     ATRU, with the regularization parameter TAU.
c
c     AINV = CP*V*(s/(s*s+tau))*UT*(bdater)-1
c
C
C     Input:
C
C     nb            - number of bins in the measured distribution;
C     nx            - number of bins in the unfolded distribution;
C     atru(nb,nx)   - the response matrix;
C     bdater(nb)    - errors in the measured distribution;
C     cp(nx,nx)     - the inverse of matrix c.
C     btil(nx)      - bdat divided by bdater;
C     uort(nb,nx)   - left orthogonal matrix of SVD;
C     sing(nx)      - vector of singular values of the matrix product A*C^-1;
C     vort(nb,nx)   - right orthogonal matrix of SVD;
C     tau           - regulatization parameter. For tau=0., one should
C                     get the exact iverse if it exists, i.e. if all
C                     singular values are nonzero.
C
C     Output:
C
C     ainv(nx,nb)   - matrix performing the rotation from the diagonalized
C                     system back to x. Prepared for subroutine gurunf;
C     unitb(nb,nb)  - product AINV*ATRU;
C     unitx(nx,nx)  - product ATRU*AINV.
C
C
      dimension atru(nb,nx),bdater(nb)
      dimension cp(nx,nx)
      dimension ainv(nx,nb)
      dimension uort(nb,nx),sing(nx),vort(nx,nx)
      dimension unitb(nb,nb), unitx(nx,nx)


      do 4110 i=1,nx
        do 4110 l=1,nb
          ainv(i,l)=0.0
          if(bdater(l).gt.0.)then
            do 4108 j=1,nx
              do 4108 k=1,nx
                ainv(i,l)=ainv(i,l)+cp(i,j)*vort(j,k)*
     *          (sing(k)/(sing(k)**2+tau))*uort(l,k)/bdater(l)
 4108       continue
          endif
 4110 continue 


C Calculate the product UNITB = ATRU * AINV. This is nb x nb matrix.

      do 4120 i=1,nb
        do 4120 k=1,nb
          unitb(i,k)=0.0
          do 4120 j=1,nx
            unitb(i,k)=unitb(i,k)+atru(i,j)*ainv(j,k)
 4120 continue 

C Calculate the product UNITX = AINV * ATRU. This is nx x nx matrix.

      do 4130 i=1,nx
        do 4130 k=1,nx
          unitx(i,k)=0.0
          do 4130 j=1,nb
            unitx(i,k)=unitx(i,k)+ainv(i,j)*atru(j,k)
 4130 continue 

      return
      end

c---------------------------------------------------------------------

      subroutine gurtan(nx,ddat,aneta,neta,dav,dis)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Determines the effective rank of the system, neta.
C
C     nx          - number of bins to be unfolded;
C     ddat(nx)    - vector containing the r.h.s. of the diagonalized
C                   system, prepared in subroutine gurrot.
C
C     Output:
C
C     aneta(nx)   - the functional to be minimized;
C     neta        - estimated effective rank of the system;
C     dav         - average of ddat(i) for neta>=i>=nx;
C     dis         - variance of ddat(i) for neta>=i>=nx.
C
C     Tentative optimal value of tau should be then calculated as:
C                      
C                  tau=(sing(neta))**2
C
C     The last (nx-neta) components of ddat should be compatible with
C     gaussian distribution with zero mean and unit variance, i.e.
C     dav should be close to zero, dis should be close to unity.
C 

      dimension ddat(nx),aneta(nx)
      eps=1.0e-6
      err=1.                         ! use ddat with equal weights
C      err=0.                         ! use ddat weighted as ddat**2
 
      anemin=999999.     
      do 300 mm=2,nx     ! perform linear fit of ln(ddat) 
      s11=0.
      s12=0.
      s22=0.
      g1 =0.
      g2 =0.
      do 200 i=1,mm
        d2=ddat(i)**2
        dl=alog(d2+eps)
        if(err.eq.1.)d2=1.
        s11=s11+d2
        s12=s12+d2*i
        s22=s22+d2*i*i
        g1 =g1 +d2*dl
        g2 =g2 +d2*dl*i
 200  continue                           ! 2*ln(ddat(i))=a1+i*a2
      dd=s11*s22-s12*s12                 
      a1=(g1*s22-g2*s12)/dd
      a2=(g2*s11-g1*s12)/dd
c      da1=sqrt(s22/dd)                   ! error in a1
c      da2=sqrt(s11/dd)                   ! error in a2
c      c12=-s12/dd/da1/da2                ! correlation coefficient

      aneta(mm)=(real(mm)+a1/a2)**2
 300  continue

      anemin=999999.     
      do 400 mm=3,nx-1
        temp=(aneta(mm-1)+2.*aneta(mm)+aneta(mm+1))/4.0
        if(temp.lt.anemin)then
          anemin=temp
          neta=mm
        endif
 400  continue

      ss1=0.0
      ss2=0.0
      do 500 i=neta,nx
        ss1=ss1+ddat(i)
        ss2=ss2+ddat(i)**2
  500 continue
      dav=ss1/(nx-neta+1)
      dis=sqrt(ss2/(nx-neta+1)-dav*dav)    
 
      return
      end

C---------------------------------------------------------

      subroutine gurunf(nx,sing,ddat,vreg,tau,isin,xi,
     &                                   xdat,xdatcm)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified 13 November 1997.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Performs the actual unfolding.
C
C     Input:
C
C     nx          - number of bins to be unfolded;
C     sing(nx)    - vector of singular values of the matrix product
C                   A*C^-1, prepared in subroutine gurrot;
C     ddat(nx)    - vector containing the r.h.s. of the diagonalized
C                   system, prepared in subroutine gurrot;
C     vreg(nx,nx) - matrix performing the rotation from the diagonalized
C                   system back to x. Prepared in subroutine gurrot;
C     tau         - the regularization parameter used as a smooth cutoff
C                   for small singular values. MUST BE POSITIVE;
C     isin        - number of sing.values exempt from regularization;
C     xi(nx)      - constant part of reg.term.
C
C     Output:
C
C     xdat(nx)      - regularized unfolded distribution;
C     xdatcm(nx,nx) - regularized covariance matrix of the unfolded 
C                     distribution. 
C
C     Note that the inverse of this matrix, xcmdm1(nx,nx), does not 
C     require regularization and is therefore calculated in gurrot.
C
      dimension sing(nx),ddat(nx),vreg(nx,nx),xi(nx)
      dimension xdat(nx),xdatcm(nx,nx)

      do 100 i=1,nx
        xdat(i)=0.0
        do 100 k=1,nx
          xdatcm(i,k)=0.0
  100 continue

      if(tau.le.0.0)return

C-------------------Calculate the unfolded distribution xdat---------


      do 854 i=1,nx                      ! x = xini*Cp*V*z    
        wdat=0.0
        do 852 k=1,nx
          if(k.le.isin)then
            reg=0.
          else
            reg=1.
          endif
          zdat = ddat(k)*sing(k)/(sing(k)**2+tau*reg)
          wdat = wdat + vreg(i,k)*zdat 
  852   continue
        xdat(i)=wdat+xi(i)
  854 continue

C-------------------Calculate the covariance matrix xdatcm-----------

      do 858 i=1,nx                      ! X = (xini*Cp*V)*Z*(xini*Cp*V)T
        do 858 j=1,nx                        
          wdatcm=0.0
          do 856 k=1,nx
            if(k.le.isin)then
              reg=0.
            else
              reg=1.
            endif
            zdater=sing(k)/(sing(k)**2+tau*reg)
            wdatcm=wdatcm + vreg(i,k)*vreg(j,k)*zdater**2
  856     continue
          xdatcm(i,j)=wdatcm 
  858 continue

      return
      end

C---------------------Calculating Chi2's--------------------

      subroutine gurchi(nx,xtst,xdat,xcmdm1,chi2,chi3)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1995. Last modified on 15 May 1996.
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Calculates  the deviation of the resulting unfolded distribution
C     xdat from the input test distribution xtst.
C
C     Input:
C
C     nx           - number of bins in test and unfolded histograms;
C     xtst(nx)     - test distribution used initially for generating
C                    btst. The "measured" distribution bdat was then
C                    obtained by addidng random gaussian errors to
C                    the latter. Errors in xtst are assumed to be 
C                    purely statistical and independent;
C     xdat(nx)     - the result obtained after unfolding bdat,
C                    output of subroutine gurunf;
C     xcmdm1(nx,nx)- the inverse of the covariance matrix of xdat,
C                    output of gurrot.
C
C     Output:
C
C     chi2         - chi-squared of xtst versus xdat. The diagonal
C                    error matrix was used with diagonal elements
C                    equal to xtst(i);
C     chi3         - chi-squared of xdat versus xtst. The full inverse
C                    covariance matrix xcmdm1 was used. 
C
C     The unfolding procedure generally results in an increase of errors
C     compared to pure statistical errors in xtst (i.e. compared to the 
C     case if the response matrix were equal to identity). The quantity 
C                      
C                        q=sqrt(chi3/chi2)
C  
C     is the average increase of errors in xdat due to unfolding.
C     

      dimension xtst(nx),xdat(nx),xcmdm1(nx,nx)

      chi2=0.0           ! chi2 of xtst vs xdat with errors of xtst.
      chi3=0.0           ! chi2 of xdat vs xtst with inv.cov.mat. of xdat
      do 9120 i=1,nx                        
        if(xtst(i).ne.0.)chi2=chi2+((xdat(i)-xtst(i))**2)/abs(xtst(i))
        do 9100 j=1,nx                        
          chi3=chi3+(xdat(i)-xtst(i))*xcmdm1(i,j)*(xdat(j)-xtst(j))
 9100   continue
 9120 continue

c      chi2=chi2/nx
c      chi3=chi3/nx

      return
      end


      subroutine guresi(nb,nx,c,apro,xdat,xini,
     +               bdat,bdater,ddat,sing,tau,isin,xi,
     +               resib,resic,resix,resid1,resid2,resid3,xeff)
C
C     Part of GURU package, version 0.97. Written by V.Kartvelishvili
C     in October 1997. Last modified 13 November 1997. 
C     GURU (General and Universal Routine for Unfolding) is the
C     Fortran realisation of the algorithm desicribed in:
C     A.Hoecker and V.Kartvelishvili, "SVD Approach to Data Unfolding",
C     hep-ph/9509307, August 95; published in NIM A 372 (1996) 469.
C
C     Calculates two versions of the residue length.
C
C     Input:
C
C     nx          - number of bins to be unfolded;
C     nb          - number of bins measured;
C     c(nx,nx)    - regularization matrix;
C     apro(nb,nx) - probability response matrix;
C     xdat(nx)    - regularized unfolded distribution;
C     xini(nx)    - initial MC distribution;
C     bdat(nb)    - measured distribution;
C     bdater(nb)  - errors in measured distribution;
C     ddat(nx)    - vector containing the rotated r.h.s. of the 
C                   diagonalized system, prepared in subroutine gurrot;    
C     sing(nx)    - vector of singular values of the matrix product
C                   A*C^-1, prepared in subroutine gurrot;
C     tau         - the regularization parameter used as a smooth cutoff
C                   for small singular values. MUST BE POSITIVE;
C     isin        - number of sing.val. exempt from regularization;
C     xi(nx)      - constant part of the regularization term.
C
C     Output:
C
C     resib(nb)   - contribution of each initial equation to
C                   the overall residue length squared;
C     resic(nx)   - partial contribution of curvature term to
C                   the overall residue length squared;
C     resix(nx)   - contribution of each diagonalized equation to
C                   the overall residue length squared;
C     resid1      - part of residue squared from |Ax-b|**2;
C     resid2      - part of residue squared from tau*|Cx|**2;
C     resid3      - same as resid1 but after diagonalization;
C     xeff        - effective number of d.o.f.
C
      dimension c(nx,nx),apro(nb,nx),xdat(nx),xini(nx),bdat(nb),xi(nx)
      dimension bdater(nb),sing(nx),ddat(nx)
      dimension resib(nb),resic(nx),resix(nx)

      real tau,resid1,resid2,resid3,delta,sum
      integer i,j,nx,nb,isin

      resid1=0.0
      do 500 i=1,nb
        sum=0.0
        do 400 j=1,nx
          sum=sum+apro(i,j)*xdat(j)
 400    continue
        if(bdater(i).gt.0.0) then                
          delta=(sum-bdat(i))/bdater(i)
        else
          delta=0.0
        endif
        resib(i)=delta**2
        resid1=resid1+delta**2
 500  continue

      resid2=0.0
      do 580 i=1,nx
        sum=0.0
        do 540 j=1,nx
          sum=sum+c(i,j)*(xdat(j)-xi(j))/xini(j)
 540    continue
        resic(i)=tau*sum**2
        resid2=resid2+tau*sum**2
 580  continue

      resid3=0.0
      xeff=0.0
      do 600 j=1,nx
        if(sing(j).ne.0.0.or.tau.ne.0.0) then
          if(j.le.isin)then
            reg=0.
          else
            reg=1
          endif
          delta3=tau*ddat(j)**2/(sing(j)**2+tau*reg)
          delta4=sing(j)**2/(sing(j)**2+tau*reg)
        else
          delta3=0.0
          delta4=0.0
        endif
        resix(j)=delta3
        resid3=resid3+delta3
        xeff=xeff+delta4
 600  continue

      return
      end


C----------------------------------------------------------------

      SUBROUTINE SVD (A, S, V, MMAX, NMAX, M, N, P, WITHU, WITHV)
C     ---------------------  START OF SVD  -----------------------------
C
      REAL       A(MMAX,1), S(NMAX), V(NMAX,NMAX)
      INTEGER    M, N, P
      LOGICAL    WITHU, WITHV
C
C     ADDITIONAL SUBROUTINE NEEDED: ROTATE
C
C     THIS SUBROUTINE COMPUTES THE SINGULAR VALUE DECOMPOSITION
C     OF A REAL M*N MATRIX A, I.E. IT COMPUTES MATRICES U, S, AND V
C     SUCH THAT
C                  A = U * S * VT ,
C     WHERE
C              U IS AN M*N MATRIX AND UT*U = I, (UT=TRANSPOSE
C                                                    OF U),
C              V IS AN N*N MATRIX AND VT*V = I, (VT=TRANSPOSE
C                                                    OF V),
C        AND   S IS AN N*N DIAGONAL MATRIX.
C
C     DESCRIPTION OF PARAMETERS:
C     A = REAL ARRAY. A CONTAINS THE MATRIX TO BE DECOMPOSED.
C         THE ORIGINAL DATA ARE LOST. IF WITHU=.TRUE., THEN
C         THE MATRIX U IS COMPUTED AND STORED IN THE ARRAY A.
C     MMAX = INTEGER VARIABLE. THE NUMBER OF ROWS IN THE
C            ARRAY A.
C     NMAX = INTEGER VARIABLE. THE NUMBER OF ROWS IN THE
C            ARRAY V.
C     M,N = INTEGER VARIABLES. THE NUMBER OF ROWS AND COLUMNS
C           IN THE MATRIX STORED IN A. (N<=M<=100. IF IT IS
C           NECESSARY TO SOLVE A LARGER PROBLEM, THEN THE
C           AMOUNT OF STORAGE ALLOCATED TO THE ARRAY T MUST
C           BE INCREASED ACCORDINGLY.) IF M<N, THEN EITHER
C           TRANSPOSE THE MATRIX A OR ADD ROWS OF ZEROS TO
C           INCREASE M TO N.
C     P = INTEGER VARIABLE. IF P>0, THEN COLUMNS N+1, . . . ,
C         N+P OF A ARE ASSUMED TO CONTAIN THE COLUMNS OF AN M*P
C         MATRIX B. THIS MATRIX IS MULTIPLIED BY UT, AND UPON
C         EXIT, A CONTAINS IN THESE SAME COLUMNS THE N*P MATRIX
C         UT*B. (P-=0)
C     WITHU, WITHV = LOGICAL VARIABLES. IF WITHU=.TRUE., THEN
C         THE MATRIX U IS COMPUTED AND STORED IN THE ARRAY A.
C         IF WITHV=.TRUE., THEN THE MATRIX V IS COMPUTED AND
C         STORED IN THE ARRAY V.
C     S = REAL ARRAY.  S(1), . . . , S(N) CONTAIN THE DIAGONAL
C         ELEMENTS OF THE MATRIX S ORDERED SO THAN S(I)>=S(I+1),
C         I=1, . . . , N-1.
C     V = REAL ARRAY.  V CONTAINS THE MATRIX V.  IF WITHU
C         AND WITHV ARE NOT BOTH =.TRUE., THEN THE ACTUAL
C         PARAMETER CORRESPONDING TO A AND V MAY BE THE SAME.
C
C     THIS SUBROUTINE IS A REAL VERSION OF A FORTRAN SUBROUTINE
C     BY BUSINGER AND GOLUB, ALGORITHM 358:  SINGULAR VALUE
C     DECOMPOSITION OF A COMPLEX MATRIX, COMM. ACM, V. 12,
C     NO. 10, PP. 564-565 (OCT. 1969).
C     WITH REVISIONS BY RC SINGLETON, MAY 1972.
C     ------------------------------------------------------------------
CCC      DIMENSION  T(100)
      DIMENSION  T(10000)
C
C+SELF,IF=CDC.
C      DATA ETA,TOL/1.E-14,1.E-279/
C+SELF,IF=IBM.
C      DATA ETA,TOL/1.E-6,1.E-72/
C+SELF,IF=-CDC,IF=-IBM.
      DATA ETA,TOL / 1.E-6,1.E-37/
C+SELF.
C     ETA AND TOL ARE MACHINE DEPENDENT CONSTANTS
C     ETA IS THE MACHINE EPSILON (RELATIVE ACCURACY);
C     TOL IS THE SMALLEST REPRESENTABLE REAL DIVIDED BY ETA.
C
      NP = N + P
      N1 = N + 1
C
C     HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM
      G = 0.0
      EPS = 0.0
      L = 1
   10 T(L) = G
      K = L
      L = L + 1
C
C     ELIMINATION OF A(I,K), I=K+1, . . . , M
      S(K) = 0.0
      Z = 0.0
      DO 20 I = K,M
   20    Z = Z + A(I,K)**2
      IF (Z.LT.TOL) GOTO 50
      G = SQRT(Z)
      F = A(K,K)
      IF (F.GE.0.0) G = - G
      S(K) = G
      H = G * (F - G)
      A(K,K) = F - G
      IF (K.EQ.NP) GOTO 50
      DO 40 J = L,NP
      F = 0
      DO 30 I = K,M
   30 F= F + A(I,K)*A(I,J)
      F = F/H
      DO 40 I = K,M
   40 A(I,J) = A(I,J) + F*A(I,K)
C
C     ELIMINATION OF A(K,J), J=K+2, . . . , N
   50 EPS = MAX(EPS,ABS(S(K)) + ABS(T(K)))
      IF (K.EQ.N) GOTO 100
      G = 0.0
      Z = 0.0
      DO 60 J = L,N
   60 Z = Z + A(K,J)**2
      IF (Z.LT.TOL) GOTO 10
      G = SQRT(Z)
      F = A (K,L)
      IF (F.GE.0.0) G = - G
      H = G * (F - G)
      A(K,L) = F - G
      DO 70 J = L,N
   70 T(J) = A(K,J)/H
      DO 90 I = L,M
      F = 0
      DO 80 J = L,N
   80 F = F + A(K,J)*A(I,J)
      DO 90 J = L,N
   90 A(I,J) = A(I,J) + F*T(J)
C
      GOTO 10
C
C     TOLERANCE FOR NEGLIGIBLE ELEMENTS
  100 EPS = EPS*ETA
C
C     ACCUMULATION OF TRANSFORMATIONS
      IF (.NOT.WITHV) GOTO 160
      K = N
      GOTO 140
  110 IF (T(L).EQ.0.0) GOTO 140
      H = A(K,L)*T(L)
      DO 130 J = L,N
      Q = 0
      DO 120 I = L,N
  120 Q = Q + A(K,I)*V(I,J)
      Q = Q/H
      DO 130 I = L,N
  130 V(I,J) = V(I,J) + Q*A(K,I)
  140 DO 150 J = 1,N
  150 V(K,J) = 0
      V(K,K) = 1.0
      L=K
      K = K - 1
      IF (K.NE.0) GOTO 110
C
  160 K = N
      IF (.NOT. WITHU) GO TO 230
      G = S (N)
      IF (G.NE.0.0) G = 1.0/G
      GO TO 210
  170 DO 180 J = L,N
  180 A(K,J) = 0
      G = S(K)
      IF (G.EQ.0.0) GOTO 210
      H = A(K,K)*G
      DO 200 J = L,N
      Q =0
      DO 190 I = L,M
  190 Q= Q + A(I,K)*A(I,J)
      Q = Q/H
      DO 200 I = K,M
  200 A(I,J) = A(I,J) + Q*A(I,K)
      G = 1.0/G
  210 DO 220 J = K,M
  220 A(J,K) = A(J,K)*G
      A(K,K) = A(K,K) + 1.0
      L = K
      K = K - 1
      IF (K.NE.0) GOTO 170
C
C     QR DIAGONALIZATION
      K = N
C
C     TEST FOR SPLIT
  230 L = K
  240 IF (ABS(T(L)).LE.EPS) GOTO 290
      L = L - 1
      IF (ABS(S(L)).GT.EPS) GOTO 240
C
C     CANCELLATION
      CS = 0.0
      SN = 1.0
      L1 = L
      L = L + 1
      DO 280 I = L,K
      F = SN*T(I)
      T(I) = CS*T(I)
      IF (ABS(F).LE.EPS) GOTO 290
      H = S(I)
      W = SQRT(F*F + H*H)
      S(I) = W
      CS = H/W
      SN = - F/W
      IF (WITHU) CALL ROTATE(A(1,L1),A(1,I), CS, SN, M)
      IF(NP.EQ.N) GO TO 280
      DO 270 J = N1,NP
      Q = A(L1,J)
      R = A(I,J)
      A(L1,J) = Q*CS + R*SN
  270 A(I,J) = R*CS - Q*SN
  280 CONTINUE
C
C     TEST FOR CONVERGENCE
  290 W = S(K)
      IF (L.EQ.K) GO TO 360
C
C     ORIGIN SHIFT
      X = S(L)
      Y = S(K-1)
      G = T(K-1)
      H = T(K)
      F = ((Y-W) * (Y+W) + (G-H) * (G+H)) / (2.0*H*Y)
      G=SQRT(F*F+1.0)
      IF (F.LT.0.0) G = - G
      F = ((X - W) * (X+W)+ (Y/(F+G) - H) *H)/X
C
C   QR STEP
      CS = 1.0
      SN = 1.0
      L1 = L + 1
      DO 350 I = L1,K
      G = T(I)
      Y = S(I)
      H = SN*G
      G = CS*G
      W = SQRT(H*H + F*F)
      T(I-1) = W
      CS = F/W
      SN = H/W
      F = X*CS + G*SN
      G = G*CS - X*SN
      H = Y*SN
      Y = Y*CS
      IF (WITHV) CALL ROTATE(V(1,I-1), V(1,I), CS, SN, N)
      W = SQRT(H*H + F*F)
      S(I-1) = W
      CS = F/W
      SN = H/W
      F = CS*G + SN*Y
      X = CS*Y - SN*G
      IF (WITHU) CALL ROTATE(A(1,I-1), A(1,I), CS, SN, M)
      IF (N.EQ.NP) GO TO 350
      DO 340 J = N1,NP
      Q = A(I-1,J)
      R = A(I,J)
      A(I-1,J) = Q*CS + R*SN
  340 A(I,J) = R*CS - Q*SN
  350 CONTINUE
C
      T(L) = 0.0
      T(K) = F
      S(K) = X
      GO TO 230
C
C   CONVERGENCE
  360 IF (W.GE.0.0) GO TO 380
      S(K) = - W
      IF (.NOT.WITHV) GO TO 380
      DO 370 J = 1,N
  370 V(J,K) = - V(J,K)
  380 K = K-1
      IF (K.NE.0) GO TO 230
C
C   SORT SINGULAR VALUES
      DO 450 K = 1,N
      G = -1.0
      DO 390 I = K,N
      IF (S(I).LT.G) GO TO 390
      G = S(I)
      J = I
  390 CONTINUE
      IF (J.EQ.K) GO TO 450
      S(J) = S(K)
      S(K) = G
      IF (.NOT.WITHV) GO TO 410
      DO 400 I = 1,N
      Q = V(I,J)
      V(I,J) = V(I,K)
  400 V(I,K) = Q
  410 IF (.NOT.WITHU) GO TO 430
      DO 420 I = 1,M
      Q = A(I,J)
      A(I,J) = A(I,K)
  420 A(I,K) = Q
  430 IF (N.EQ.NP) GO TO 450
      DO 440 I = N1,NP
      Q = A(J,I)
      A(J,I) = A(K,I)
  440 A(K,I) = Q
  450 CONTINUE
C
      RETURN
      END
C+DECK,ROTATE.
      SUBROUTINE ROTATE (X, Y, CS, SN, N)
      REAL X(N), Y (N), CS, SN
C
      DO 10 J = 1,N
      XX = X(J)
      X(J) = XX*CS + Y(J)*SN
   10 Y(J) = Y(J)*CS - XX*SN
      RETURN
C     --------------------  END OF SVD ---------------------------------
      END