SUBROUTINE INV(MA,IA,N,TOL)                                     
      IMPLICIT REAL*8 (A-H,O-Z)                                       
      IMPLICIT INTEGER*4(I-N)
      REAL*8 MA,IA                                                    
      INTEGER*4 oo                                                      
      DIMENSION IA(N,N),MA(N,N)
      real*8, allocatable ::E(:,:)
      write(6,*)N*2*N/8,' bytes',N,2*N
      allocate(E(N,2*N))
      NMAT=N                                                          
      DO 1  ii=1,nmat                                                 
      DO 1  jj=1,nmat                                                 
      e(ii,jj)=ma(ii,jj)
      E(II,JJ+NMAT)=0.0D0                                             
      if (ii.EQ.jj)e(ii,jj+nmat)=1.0D0                                
1     CONTINUE                                                        
      DO 2 ii=1,nmat-1                                                
      IW=II                                                           
      if (DABS(e(iw,iw)).LT.TOL) then                                 
       DO 3 io=iw+1,nmat                                              
       oo=IO
       if (io.GT.nmat)  oo=io-nmat                                    
       if (DABS(e(oo,iw)).GE.TOL) goto 11                             
3      CONTINUE                                                       
       do i=1,N
       write(6,90)(MA(i,j),j=1,N)
       enddo
90     format(12g12.4)
       WRITE(*,*)' INVERSE MATRIX CAN NOT BE FOUND '                  
       STOP                                                           
11     CONTINUE                                                       
       DO 4 kk=1, 2*nmat                                              
       w=e(iw,kk)                                                     
       e(iw,kk)=e(oo,kk)                                              
4      e(oo,kk)=w                                                     
      ENDIF                                                           
      DO 5 jj=ii+1,nmat                                               
      DO 6 kk=ii+1, 2*nmat
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii)                    
5     e(jj,ii)=0.0D0                                                  
2     CONTINUE                                                        
      DO 7  ii=1, nmat-1
      i2=nmat-ii+1                                                    
      DO 8 jj=1,i2-1                                                  
      j2=i2-jj                                                        
      DO 9 kk=1, nmat                                                 
9     e(j2,kk+nmat)=e(j2,kk+nmat)-                                    
     1e(i2,kk+nmat)*e(j2,i2)/e(i2,i2)                                 
      e(j2,i2)=0.0D0
 8    CONTINUE                                                        
 7    CONTINUE                                                        
      DO 10 ii=1,nmat                                                 
      DO 12  jj=1,nmat                                                
12    iA(ii,jj)=e(ii,jj+nmat)/e(ii,ii)                                
10    CONTINUE                                                        
      deallocate(E)
      RETURN                                                          
      END