subroutine axlbx2(n,a,ia,b,ib)
c
c this subroutine reduces the symmetric positive definite
c tridiagonal matrix b to the identity matrix while keeping the
c symmetric tridiagonal matrix a tridiagonal. 
c The diagonal of the matrix a
c is stored in the last(second) column of the a array.
c similar with b. 
c Elements 2 through n of their first columns contains the
c the off diagonal elements
c This code is based on Crawford's algorithm, but it incorpoprates
c a burn at both ends technique which halves
c operations.
c ia and ib should be at least n.
c
c on output the second column of a is the diagonal d for tql1
c and the first column is the off diagonal e
        integer n,ia,ib,iz,nz
        real a(ia,2),b(ib,2)
        real s,c,e,g,ce,sg,sh,ch,dia,x,y,p
        integer i,j,k,k1,np1,ibb,nm1
c
c do i=n,n-1,...1
c
       np1=n+1
       k=2
       k1=1
       nm1=n-1
       nov2=n/2
       do 100 ibb=1,nov2
          i=np1-ibb
c
c apply cholesky to b
c
         x=sqrt(b(i,k))
         y=b(i,k1)/x
         p=a(i,k1)/x
         a(i,k)=a(i,k)/b(i,k)
         a(i,k1)=p-y*a(i,k)
         a(i-1,k)=a(i-1,k)-y*(p+a(i,k1))
         b(i-1,k)=b(i-1,k)-y*y
         if (i.eq.ibb+1) go to 10
         x2=sqrt(b(ibb,k))
         ibbp1=ibb+1
         y2=b(ibbp1,k1)/x2
         p2=a(ibbp1,k1)/x2
         a(ibb,k)=a(ibb,k)/b(ibb,k)
         a(ibbp1,k1)=p2-y2*a(ibb,k)
         a(ibbp1,k)=a(ibbp1,k)-y2*(p2+a(ibbp1,k1))
         b(ibbp1,k)=b(ibbp1,k)-y2*y2
10        if (i.eq.n) go to 100
         a(i+1,k1)=a(i+1,k1)/x
c
c generate unwanted element in f and chase it down the matrix
c
         f=-y*a(i+1,k1)
           if (ibb+1.eq.i) go to 15
         a(ibb,k1)=a(ibb,k1)/x2
         f2=-y2*a(ibb,k1)
         jbb=ibb
15      continue
         do 50 j=i,nm1
            sc=f
            if(abs(f).le.abs(a(j,k1)))go to 20
                  r=a(j,k1)/f
            go to 30
  20         sc=a(j,k1)
             r=f/a(j,k1)
  30        continue
            dia=sc*sqrt(1.0+r*r)
            c=a(j,k1)/dia
            s=f/dia
            a(j,k1)=dia
            e=a(j,k)
            g=a(j+1,k)
            ce=c*e
            sg=s*g
            sh=s*a(j+1,k1)
            ch=c*a(j+1,k1)
            a(j,k)=c*ce+2*s*ch+s*sg
            a(j+1,k)=e+g-a(j,k)
            a(j+1,k1)=s*(c*(g-e)-sh)+c*ch
            if (i.eq.ibb+1) go to 45
            sc2=f2
            if(abs(f2).le.abs(a(jbb+1,k1)))go to 22
                  r2=a(jbb+1,k1)/f2
            go to 32
  22         sc2=a(jbb+1,k1)
             r2=f2/a(jbb+1,k1)
  32        continue
            dia2=sc2*sqrt(1.0+r2*r2)
            c2=a(jbb+1,k1)/dia2
            s2=-f2/dia2
            a(jbb+1,k1)=dia2
            e=a(jbb-1,k)
            g=a(jbb,k)
            c2e=c2*e
            s2g=s2*g
            s2h=s2*a(jbb,k1)
            c2h=c2*a(jbb,k1)
            a(jbb-1,k)=c2*c2e+2*s2*c2h+s2*s2g
            a(jbb,k)=e+g-a(jbb-1,k)
            a(jbb,k1)=s2*(c2*(g-e)-s2h)+c2*c2h
            jbb=jbb-1
45          if (j.eq.nm1) go to 50
               f=s*a(j+2,k1)
               a(j+2,k1)=c*a(j+2,k1)
               if (i.eq.ibb+1) go to 50
               f2=-s2*a(jbb,k1)
               a(jbb,k1)=c2*a(jbb,k1)
50          continue
100         continue
        n2=(n+1)/2
        x=sqrt(b(n2,k))
        a(n2,k1)=a(n2,k1)/x
        a(n2+1,k1)=a(n2+1,k1)/x
        a(n2,k)=a(n2,k)/b(n2,k)
        return
        end

.