subroutine ludcmp_dp(a,n,np,indx,d)

c       jul 2011 malv+fgg

      implicit none

      integer n, np
      real*8 a(np,np), d
      integer indx(n)
      
      integer nmax, i, j, k, imax
      real*8 tiny
      parameter (nmax=100,tiny=1.0d-20)                                       
      real*8 vv(nmax), aamax, sum, dum


      d=1.0d0
      do 12 i=1,n                                                             
        aamax=0.0d0
        do 11 j=1,n                                                           
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))                         
11      continue                                                              
        if (aamax.eq.0.0) pause 'singular matrix.'                          
        vv(i)=1.0d0/aamax                          
12    continue                                                                
      do 19 j=1,n                                                             
        if (j.gt.1) then                                                      
          do 14 i=1,j-1                                                       
            sum=a(i,j)                                                        
            if (i.gt.1)then                                                   
              do 13 k=1,i-1                                                   
                sum=sum-a(i,k)*a(k,j)                                         
13            continue                                                        
              a(i,j)=sum                                                      
            endif                                                             
14        continue                                                            
        endif                                                                 
        aamax=0.0d0                                                           
        do 16 i=j,n                                                           
          sum=a(i,j)                                                          
          if (j.gt.1)then                                                     
            do 15 k=1,j-1                                                     
              sum=sum-a(i,k)*a(k,j)                                           
15          continue                                                          
            a(i,j)=sum                                                        
          endif                                                               
          dum=vv(i)*abs(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(j.ne.n)then                                                        
          if(a(j,j).eq.0.0)a(j,j)=tiny                                     
          dum=1.0d0/a(j,j)                                                  
          do 18 i=j+1,n                                                       
            a(i,j)=a(i,j)*dum                                                 
18        continue                                                            
        endif                                                                 
19    continue                                                                
      if(a(n,n).eq.0.0)a(n,n)=tiny                                     
      return                                                                  
      end