	subroutine dgesv_compareLU(NMAX,a,b,x,indx,ipiv,tlu)
        implicit none
        integer NMAX
        real*8 b(NMAX),x(NMAX),a(NMAX,NMAX)
        real*8 err
        real*8 myrandom
        real*8 tzero,tnow
        real*8 tlapack
        real*8 csecond
        real*8 tmin,tmax,tave,titer
        integer i,j,info,lda,ldb,nrhs,n,k
        integer ipiv(NMAX)
        integer iseed
        integer niter
        integer indx(NMAX)
        real*8 d,tlu
       
        tlapack=0.0
        n=NMAX
         
!  initialize random numbers generator
        call init_randoms_generator(iseed)

! repeat few times : niter = 1
        niter = 1
        tmax=-1.0
        tmin=1.00000000e37
        do k = 1, niter
!  fill array A, vector x with random numbers
        do i = 1,NMAX
           do j = 1,NMAX
              a(j,i)=myrandom(iseed)
           enddo
           x(i)=myrandom(iseed)
        enddo
!   create vector b from known A and x : b=Ax
        do i = 1,NMAX
           b(i)=0.0000000
           do j = 1, NMAX
              b(i)=b(i) + a(j,i) * x(j)
           enddo
        enddo

        nrhs = 1
        lda = NMAX
        ldb = NMAX

        
        tzero=csecond()
            call ludcmp(a,NMAX,NMAX,indx,d)
            call lubksb(a,NMAX,NMAX,indx,b)
        tnow=csecond()
        titer=tnow-tzero
        if (titer.gt.tmax) then 
              tmax=titer
        endif
        if (titer.lt.tmin) then 
              tmin=titer
        endif
        tlapack=tlapack+titer
        enddo

        tave = tlapack/dble(niter)
        tlu = tave
        return
        end

      SUBROUTINE lubksb(a,n,np,indx,b)
	implicit real*8 (a-h,o-z)
      INTEGER n,np,indx(n)
      REAL*8 a(np,np),b(n)
      INTEGER i,ii,j,ll
      REAL*8 sum
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        do 13 j=i+1,n
          sum=sum-a(i,j)*b(j)
13      continue
        b(i)=sum/a(i,i)
14    continue
      return
      END

      SUBROUTINE ludcmp(a,n,np,indx,d)
	implicit real*8 (a-h,o-z)
      INTEGER n,np,indx(n),NMAX
      REAL*8 d,a(np,np),TINY
      PARAMETER (NMAX=16384,TINY=1.0e-20)
      INTEGER i,imax,j,k
      REAL*8 aamax,dum,sum,vv(NMAX)
      d=1.
      do 12 i=1,n
        aamax=0.
        do 11 j=1,n
          if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
11      continue
        if (aamax.eq.0.) pause 'singular matrix in ludcmp'
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.
        do 16 i=j,n
          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*dabs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0.)a(j,j)=TINY
        if(j.ne.n)then
          dum=1./a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      return
      END
