FUNCTION IDENTITY(n)
   implicit none
   integer n,i
   real(8), dimension(n,n) ::identity

   identity=0.0
   do i=1,n
      identity(i,i)=1.0d0
   end do
END FUNCTION IDENTITY



SUBROUTINE INVERSE(A,matrix_ill)
   implicit none
   real(8), dimension(:,:)                         ::A
   real(8), dimension(size(A,dim=1),size(A,dim=2)) ::inv,L,U
   real(8), dimension(size(A,dim=1))               ::S1,S2
   real(8)  eps
   integer  n,i,c,k
   logical  matrix_ill

   interface
      function IDENTITY(n)
         integer n
         real(8), dimension(n,n):: IDENTITY
      end function IDENTITY
   end interface

   matrix_ill=.false.
   n=size(A,dim=1)
   if(n.ne.size(A,dim=2)) then
      write(*,*)'matrix is not square - cannot invert'
      write(*,*)'PROGRAM STOPPED'
      stop
   end if
   eps=epsilon(0.0d0)

!  do LU factorisation of matrix A
!  this is a simple LU decomposition which always assumes
!  the diagonal elements are non-zero.

   L=0.0d0
   U=0.0d0
   L(1:n,1)=A(1:n,1)
   if(abs(L(1,1)).lt.eps) then
      write(*,*)'matrix ill conditioned',abs(L(1,1))
      matrix_ill=.true.
      goto 90
   end if
   U(1,1:n)=A(1,1:n)/L(1,1)

   do c=2,N
      S1=0.0d0
      S2=0.0d0
      do k=1,c-1
         S1(:)=S1(:)+L(:,k)*U(k,c)
         S2(:)=S2(:)+L(c,k)*U(k,:)
      end do
      L(:,c)=A(:,c)-S1(:)
      if(abs(L(c,c)).lt.eps) then
         write(*,*)'matrix ill conditioned',abs(L(c,c))
         matrix_ill=.true.
         goto 90
      end if
      U(c,:)=(A(c,:)-S2(:))/L(c,c)
   end do

!  solve for inverse in    A*inverse=I
   inv=identity(N)
   inv(1,:)=inv(1,:)/L(1,1)
   do i=2,N
      s1(:)=0.0d0
      do k=1,i-1
         s1(:)=s1(:)+L(i,k)*inv(k,:)
      end do
      inv(i,:)=(inv(i,:)-s1(:))/L(i,i)
   end do

   do i=N-1,1,-1
      s1(:)=0.0d0
      do k=i+1,N
         s1(:)=s1(:)+U(i,k)*inv(k,:)
      end do
      inv(i,:)=inv(i,:)-s1(:)
   end do
   A=inv
90 return
END SUBROUTINE INVERSE



FUNCTION MTX_CROSS_PRODUCT(A,B)
   implicit none
   real(8), dimension(:,:)              ::A
   real(8), dimension(:,:)              ::B
   real(8), dimension(size(A,dim=1)*size(B,dim=1),size(      &
            A,dim=2)*size(B,dim=2))     ::mtx_cross_product
   integer i,j,dB1,dB2

   dB1=size(B,dim=1)
   dB2=size(B,dim=2)
   do i=1,size(A,dim=1)
      do j=1,size(A,dim=2)
         mtx_cross_product((i-1)*dB1+1:i*dB1,                &
                           (j-1)*dB2+1:j*dB2)=A(i,j)*B
      end do
   end do
END FUNCTION MTX_CROSS_PRODUCT



FUNCTION NORM2(x)
   implicit none
   real(8), dimension(:) ::x
   real(8)  norm2

   norm2=sqrt(dot_product(x,x))
END FUNCTION NORM2



FUNCTION OUTER_PRODUCT(x,y)
   implicit none
   real(8),  dimension(:)                     ::x
   real(8),  dimension(:)                     ::y
   real(8),  dimension(size(x),1)             ::x1
   real(8),  dimension(1,size(y))             ::y1
   real(8),  dimension(size(x),size(y))       ::outer_product

   x1(:,1)=x
   y1(1,:)=y
   outer_product=matmul(x1,y1)
END FUNCTION OUTER_PRODUCT