!  compile with
!   gfortran -o CalcTaylorMatrices  CalcTaylorMatrices.f90
!   gfortran -O0 -g -fbounds-check -o CalcTaylorMatrices CalcTaylorMatrices.f90

module var

  real, allocatable :: b( : , : ) , n( : , : )  !  directions, planes
  real, allocatable :: Ematrix( : , : , : )  !  matrices
  real, allocatable :: TransposeMatrix( : , : , : )  !  matrices for stress

  integer :: NumberSystems

  logical :: verbose = .TRUE.
  
end module var

module var2

  INTEGER , allocatable :: combext( : , : )
  integer :: counter

end module var2

!!$module operations
!!$use types, only: dp
!!$implicit none
!!$private
!!$public scalar
!!$
!!$contains

!  --------------------
  real function scalar( a , b )
    implicit none
    real, intent(in) :: a(3), b(3)
    !real, intent(out) :: scalar
    integer :: i
    scalar = 0.
    do 100, i = 1 , 3
       scalar = scalar + a(i) * b(i)
100 end do
!!$  return
  end function scalar

!!$end module

!!$!.............
!!$module integrals
!!$use types, only: dp
!!$implicit none
!!$private
!!$public simpson
!!$
!!$contains
!!$
!!$real(dp) function simpson(f, a, b) result(s)
!!$real(dp), intent(in) :: a, b
!!$interface
!!$    real(dp) function func(x)
!!$    use types, only: dp
!!$    implicit none
!!$    real(dp), intent(in) :: x
!!$    end function
!!$end interface
!!$procedure(func) :: f
!!$s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
!!$end function
!!$
!!$end module integrals

!.............
program CalcTaylorMatrices

  use var
  use var2
  
  implicit none
  
  ! local variables
  integer :: i , j , k , ierror
  integer :: icount, jcount
  integer, allocatable :: MarkedAsDuplicate( : ) , MarkedAsSingular( : )
  
  real :: E1(5,5) , E2(5,5) , difference , det
  real :: epsilon = 1.e-5

  real :: Bcoeff,Ccoeff,Fcoeff,Gcoeff,Hcoeff
  real :: BH_stress_list( 6 , 1000 )
  real :: sqr6 = sqrt(6.)

  !  preamble
  print*, 'This program is written by AD Rollett, last edited 6 March 2016'
  print*, 'It is intended to compute all the available Taylor matrices'
  print*, 'See lecture notes on Multiple Slip (27-750) for details'

  ! code
  call SetUpFCC  ! systems for FCC

  allocate( combext( 5 , 2000 ) , STAT = ierror )  !  only needs sets of 5
  if ( ierror /= 0 ) stop 'problem allocating combext'
  call combinations ( 5 , 12 )  !  computes 5C12
  !  num. combos in counter; combext contains the list
  if(verbose) print * , '   number of combinations = ' , counter

  allocate( Ematrix( 5 , 5 , counter ) , STAT = ierror )
  if ( ierror /= 0 ) stop 'problem allocating Ematrix'
  allocate( TransposeMatrix( 5 , 5 , counter ) , STAT = ierror )
  if ( ierror /= 0 ) stop 'problem allocating TransposeMatrix'

  icount = 0 ; jcount = 0
  if( verbose ) print*,'starting up matrix calculation'
  do i = 1 , counter
     icount = icount + 1
     do j = 1 , 5  !  which slip system, the column in E
        !do k = i , 5  !  which row in E
        !  see slide 47 in L11-PolyXtal_plast-Aniso3-25Feb16.pptx
        Ematrix( 1 , j , icount ) = b( 2 , combext( j , icount ) ) * n( 2 , combext( j , icount ) )
        TransposeMatrix( j , 1 , icount ) = Ematrix( 1 , j , icount )
        Ematrix( 2 , j , icount ) = b( 3 , combext( j , icount ) ) * n( 3 , combext( j , icount ) )
        TransposeMatrix( j , 2 , icount ) = Ematrix( 2 , j , icount )
        Ematrix( 3 , j , icount ) = b( 2 , combext( j , icount ) ) * n( 3 , combext( j , icount ) ) &
                                 +  b( 3 , combext( j , icount ) ) * n( 2 , combext( j , icount ) )
        TransposeMatrix( j , 3 , icount ) = Ematrix( 3 , j , icount )
        Ematrix( 4 , j , icount ) = b( 1 , combext( j , icount ) ) * n( 3 , combext( j , icount ) ) &
                                 +  b( 3 , combext( j , icount ) ) * n( 1 , combext( j , icount ) )
        TransposeMatrix( j , 4 , icount ) = Ematrix( 4 , j , icount )
        Ematrix( 5 , j , icount ) = b( 2 , combext( j , icount ) ) * n( 1 , combext( j , icount ) ) &
                                 +  b( 1 , combext( j , icount ) ) * n( 2 , combext( j , icount ) )
        TransposeMatrix( j , 5 , icount ) = Ematrix( 5 , j , icount )
     end do
  end do
  if(verbose) print * , '   number of E matrices = ' , counter
  if(verbose) print * , ' now to go through and eliminate singular cases'
  allocate( MarkedAsDuplicate(icount) , STAT = ierror )
  if ( ierror /= 0 ) stop 'problem allocating MarkedAsDuplicate'
  MarkedAsDuplicate = 0
  allocate( MarkedAsSingular(icount) , STAT = ierror )
  if ( ierror /= 0 ) stop 'problem allocating MarkedAsSingular'
  MarkedAsSingular = 0
  
  do i = 1 , counter-1
     E1 = Ematrix( : , : , i )
     do j = i+1 , counter
        E2 = Ematrix( : , : , j )
        call DifferenceMatrix( 5, 5 , E1 , E2 , difference )
        if ( difference < epsilon ) then
           print*,'Found redundant case for i,j = ', i , j
           MarkedAsDuplicate = 99
        end if
     end do
  end do

  open ( unit = 4, file = 'Stress-matrix-list.txt', status = 'unknown' )
  !write(4,*) ' number  A  B  C  F  G  H'
  jcount = 0
  do i = 1 , counter
     MarkedAsSingular(i) = 0
     E1 = TransposeMatrix( : , : , i )
     call matinv5( E1 , E2 , det )  !  here E2 is the inverted matrix
     if ( abs(det) < epsilon ) then
        print*,'Found singular case for i = ', i, ' determinant = ', det
        MarkedAsSingular(i) = 99
        jcount = jcount + 1
     else
        !  matrix is good, compute stresses
        Ccoeff = -1. * ( E2(1,1) + E2(1,2) + E2(1,3) + E2(1,4) + E2(1,5) ) / sqr6
        Bcoeff = ( E2(2,1) + E2(2,2) + E2(2,3) + E2(2,4) + E2(2,5) )  / sqr6
        Fcoeff = ( E2(3,1) + E2(3,2) + E2(3,3) + E2(3,4) + E2(3,5) ) / sqr6
        Gcoeff = ( E2(4,1) + E2(4,2) + E2(4,3) + E2(4,4) + E2(4,5) ) / sqr6
        Hcoeff = ( E2(5,1) + E2(5,2) + E2(5,3) + E2(5,4) + E2(5,5) ) / sqr6
        print*,' for # ', i, ' stress values are ', Bcoeff,Ccoeff,Fcoeff,Gcoeff,Hcoeff
        BH_stress_list( 2 , i ) = Bcoeff
        BH_stress_list( 3 , i ) = Ccoeff
        BH_stress_list( 4 , i ) = Fcoeff
        BH_stress_list( 5 , i ) = Gcoeff
        BH_stress_list( 6 , i ) = Hcoeff
        BH_stress_list( 1 , i ) = (-1. * BH_stress_list( 2 , i )) - BH_stress_list( 3 , i )
        write(4,"(/'matrix # ',i6 / 5(2x,f9.4) )" ) i , E1( 1 , : )
        write(4,"(5(2x,f9.4) )" ) E1( 2 , : )
        write(4,"(5(2x,f9.4) )" ) E1( 3 , : )
        write(4,"(5(2x,f9.4) )" ) E1( 4 , : )
        write(4,"(5(2x,f9.4) )" ) E1( 5 , : )
     end if
  end do
  close( 4 )
  print * , 'Found ', jcount, ' singular cases'
  print * , 'Total of  ', counter-jcount, ' non-singular E matrices'

  open ( unit = 4, file = 'Bishop-Hill-list.txt', status = 'unknown' )
  write(4,*) ' number  A  B  C  F  G  H'
  do i = 1 , counter
     if ( MarkedAsSingular(i) == 0 ) then
        write(4,"(i4,1x,6f8.2)" ) i , BH_stress_list( : , i )
     end if
  end do
  close ( 4 )
  
  call exit (  0 )
  
end program CalcTaylorMatrices

! .......
subroutine DifferenceMatrix( m, n , A1 , A2 , difference )
  ! stuff to compute difference between 2 matrices
  ! http://www.fortran90.org/src/best-practices.html

  integer, intent(in) :: m, n
  real, intent(in) :: A1( m , n ) , A2( m , n )
  real, intent( out ) :: difference
  integer :: i , j
  difference = 0.
  do i = 1 , m
     do j = 1 , n
        difference = difference + abs( A1(i,j) - A2(i,j) )
     end do
  end do
    
end subroutine DifferenceMatrix

! .......
subroutine SetUpFCC

  use var
!!$  use operations, only: scalar
  implicit none
  real :: transfer(3), transfer2(3), dot
  real :: scalar

  integer :: i , j , ierror

  if ( verbose ) print * , 'setting up fcc systems'
  
  NumberSystems = 12
  allocate( b( 3 , NumberSystems ) , STAT = ierror )
  allocate( n( 3 , NumberSystems ) , STAT = ierror )
  b( 1 , 1 ) = -1 ; b( 2 , 1 ) = 1 ; b( 3 , 1 ) = 0
  b( 1 , 2 ) = 1  ; b( 2 , 2 ) = 0 ; b( 3 , 2 ) = 1
  b( 1 , 3 ) = 0  ; b( 2 , 3 ) = 1 ; b( 3 , 3 ) = 1
  n( 1 , 1 ) = 1  ; n( 2 , 1 ) = 1 ; n( 3 , 1 ) = -1
  n( 1 , 2 ) = 1  ; n( 2 , 2 ) = 1 ; n( 3 , 2 ) = -1
  n( 1 , 3 ) = 1  ; n( 2 , 3 ) = 1 ; n( 3 , 3 ) = -1
  !
  b( 1 , 4 ) = -1 ; b( 2 , 4 ) = 1  ; b( 3 , 4 ) = 0
  b( 1 , 5 ) = -1 ; b( 2 , 5 ) = 0  ; b( 3 , 5 ) = 1
  b( 1 , 6 ) = 0  ; b( 2 , 6 ) = -1 ; b( 3 , 6 ) = 1
  n( 1 , 4 ) = 1  ; n( 2 , 4 ) = 1  ; n( 3 , 4 ) = 1
  n( 1 , 5 ) = 1  ; n( 2 , 5 ) = 1  ; n( 3 , 5 ) = 1
  n( 1 , 6 ) = 1  ; n( 2 , 6 ) = 1  ; n( 3 , 6 ) = 1
  !
  b( 1 , 7 ) = 1 ; b( 2 , 7 ) = 1  ; b( 3 , 7 ) = 0
  b( 1 , 8 ) = 1 ; b( 2 , 8 ) = 0  ; b( 3 , 8 ) = -1
  b( 1 , 9 ) = 0 ; b( 2 , 9 ) = 1  ; b( 3 , 9 ) = 1
  n( 1 , 7 ) = 1 ; n( 2 , 7 ) = -1 ; n( 3 , 7 ) = 1
  n( 1 , 8 ) = 1 ; n( 2 , 8 ) = -1 ; n( 3 , 8 ) = 1
  n( 1 , 9 ) = 1 ; n( 2 , 9 ) = -1 ; n( 3 , 9 ) = 1
  !
  b( 1 , 10 ) = 1 ; b( 2 , 10 ) = 1  ; b( 3 , 10 ) = 0
  b( 1 , 11 ) = 1 ; b( 2 , 11 ) = 0  ; b( 3 , 11 ) = 1
  b( 1 , 12 ) = 0 ; b( 2 , 12 ) = -1 ; b( 3 , 12 ) = 1
  n( 1 , 10 ) = 1 ; n( 2 , 10 ) = -1 ; n( 3 , 10 ) = -1
  n( 1 , 11 ) = 1 ; n( 2 , 11 ) = -1 ; n( 3 , 11 ) = -1
  n( 1 , 12 ) = 1 ; n( 2 , 12 ) = -1 ; n( 3 , 12 ) = -1

  do j = 1 , NumberSystems
     transfer = b( : , j )
     call vecnorm ( transfer )
     b( : , j ) = transfer
     transfer = n( : , j )
     call vecnorm ( transfer )
     n( : , j ) = transfer
  end do

  if (verbose ) then
     do i = 1 , NumberSystems
        transfer = b( : , i )
        transfer2 = n( : , i )
        !call scalar(dot, transfer,transfer2)
        print*,' system 3 ',i, ' b , n dot_product', b(:,i), n(:,i), scalar(transfer,transfer2)
     end do
  end if
  
end subroutine SetUpFCC

! .......
subroutine vecnorm( vec )
  implicit none
  real :: vec(3),rnorm
  
  integer :: i

  rnorm = 0.
  do 10, i = 1 , 3
     rnorm = rnorm + vec(i)**2
10 end do
  if ( rnorm <= 0.0 ) return
  rnorm = sqrt(rnorm)
  do 20, i = 1 , 3
     vec(i) = vec(i) / rnorm
20 end do
  return
end subroutine vecnorm

!................
subroutine combinations ( r , nn )

  use var2 

  !INTEGER :: r=4,nn=32
  INTEGER :: r , nn
  INTEGER , allocatable :: comb(:)
  !integer :: counter

  allocate( comb( r ) )
  counter = 0
  CALL iterate(1,nn-r+1,1)

  deallocate ( comb ) 

CONTAINS

  RECURSIVE SUBROUTINE iterate(s,e,j)
    INTEGER , INTENT(IN) :: s,e
    integer :: k
    DO i = s , e
       comb( j ) = i
       IF( j < r ) CALL iterate(i+1,e+1,j+1)
       IF( j == r ) then
          !WRITE(*,*) comb
          counter = counter + 1
          do k = 1 , r
             combext( k , counter ) = comb( k )
          end do
       end if
    END do
  END SUBROUTINE iterate

END SUBROUTINE combinations

!...............................
subroutine matinv5( mat , mati , det )

  implicit none

  real mat(5,5) , mati(5,5) , minv1(5) , det
  real start1 , start2
  !	double precision denom
  !	double precision mat(4,4), mati(4,4)
  double precision a1,b1,c1,d1,e1,a2,b2,c2,d2,e2
  double precision a3,b3,c3,d3,e3,a4,b4,c4,d4,e4
  double precision a5,b5,c5,d5,e5

  a1=mat(1,1)
  b1=mat(2,1)
  c1=mat(3,1)
  d1=mat(4,1)
  e1=mat(5,1)
  a2=mat(1,2)
  b2=mat(2,2)
  c2=mat(3,2)
  d2=mat(4,2)
  e2=mat(5,2)
  a3=mat(1,3)
  b3=mat(2,3)
  c3=mat(3,3)
  d3=mat(4,3)
  e3=mat(5,3)
  a4=mat(1,4)
  b4=mat(2,4)
  c4=mat(3,4)
  d4=mat(4,4)
  e4=mat(5,4)
  a5=mat(1,5)
  b5=mat(2,5)
  c5=mat(3,5)
  d5=mat(4,5)
  e5=mat(5,5)

  !	call cpu_time(start1)

  !	do i=1,steps

  det=a5*(b4*(c3*d2*e1-c2*d3*e1-c3*d1*e2+c1*d3  &
       *e2+c2*d1*e3-c1*d2*e3)+b3*(-c4*d2*e1+c2*d4  &
       *e1+c4*d1*e2-c1*d4*e2-c2*d1*e4+c1*d2*e4)+b2  &
       *(c4*d3*e1-c3*d4*e1-c4*d1*e3+c1*d4*e3+c3*d1  &
       *e4-c1*d3*e4)+b1*(-c4*d3*e2+c3*d4*e2+c4*d2  &
       *e3-c2*d4*e3-c3*d2*e4+c2*d3*e4))+a4*(b5*(-c3  &
       *d2*e1+c2*d3*e1+c3*d1*e2-c1*d3*e2-c2*d1*e3  &
       +c1*d2*e3)+b3*(c5*d2*e1-c2*d5*e1-c5*d1*e2  &
       +c1*d5*e2+c2*d1*e5-c1*d2*e5)+b2*(-c5*d3*e1  &
       +c3*d5*e1+c5*d1*e3-c1*d5*e3-c3*d1*e5+c1*d3  &
       *e5)+b1*(c5*d3*e2-c3*d5*e2-c5*d2*e3+c2*d5  &
       *e3+c3*d2*e5-c2*d3*e5))+a3*(b5*(c4*d2*e1-c2  &
       *d4*e1-c4*d1*e2+c1*d4*e2+c2*d1*e4-c1*d2*e4)  &
       +b4*(-c5*d2*e1+c2*d5*e1+c5*d1*e2-c1*d5*e2  &
       -c2*d1*e5+c1*d2*e5)+b2*(c5*d4*e1-c4*d5*e1  &
       -c5*d1*e4+c1*d5*e4+c4*d1*e5-c1*d4*e5)  &
       +b1*(-c5*d4*e2+c4*d5*e2+c5*d2*e4-c2*d5*e4  &
       -c4*d2*e5+c2*d4*e5))+a2*(b5*(-c4*d3*e1+c3*d4  &
       *e1+c4*d1*e3-c1*d4*e3-c3*d1*e4+c1*d3*e4)  &
       +b4*(c5*d3*e1-c3*d5*e1-c5*d1*e3+c1*d5*e3+c3  &
       *d1*e5-c1*d3*e5)+b3*(-c5*d4*e1+c4*d5*e1+c5  &
       *d1*e4-c1*d5*e4-c4*d1*e5+c1*d4*e5)+b1*(c5  &
       *d4*e3-c4*d5*e3-c5*d3*e4+c3*d5*e4+c4*d3*e5  &
       -c3*d4*e5))+a1*(b5*(c4*d3*e2-c3*d4*e2-c4*d2  &
       *e3+c2*d4*e3+c3*d2*e4-c2*d3*e4)+b4*(-c5*d3  &
       *e2+c3*d5*e2+c5*d2*e3-c2*d5*e3-c3*d2*e5+c2  &
       *d3*e5)+b3*(c5*d4*e2-c4*d5*e2-c5*d2*e4+c2  &
       *d5*e4+c4*d2*e5-c2*d4*e5)+b2*(-c5*d4*e3+c4  &
       *d5*e3+c5*d3*e4-c3*d5*e4-c4*d3*e5+c3*d4*e5))

  if(det.eq.0) return

  mati(1,1)=(b5*c4*d3*e2-b4*c5*d3*e2-b5*c3  &
       *d4*e2+b3*c5*d4*e2+b4*c3*d5*e2-b3*c4*d5  &
       *e2-b5*c4*d2*e3+b4*c5*d2*e3+b5*c2*d4*e3  &
       -b2*c5*d4*e3-b4*c2*d5*e3+b2*c4*d5*e3+b5  &
       *c3*d2*e4-b3*c5*d2*e4-b5*c2*d3*e4+b2*c5  &
       *d3*e4+b3*c2*d5*e4-b2*c3*d5*e4-b4*c3*d2  &
       *e5+b3*c4*d2*e5+b4*c2*d3*e5-b2*c4*d3*e5  &
       -b3*c2*d4*e5+b2*c3*d4*e5)/det
  mati(1,2)=(-a5*c4*d3*e2+a4*c5*d3*e2+a5*c3  &
       *d4*e2-a3*c5*d4*e2-a4*c3*d5*e2+a3*c4*d5  &
       *e2+a5*c4*d2*e3-a4*c5*d2*e3-a5*c2*d4*e3  &
       +a2*c5*d4*e3+a4*c2*d5*e3-a2*c4*d5*e3-a5  &
       *c3*d2*e4+a3*c5*d2*e4+a5*c2*d3*e4-a2*c5  &
       *d3*e4-a3*c2*d5*e4+a2*c3*d5*e4+a4*c3*d2  &
       *e5-a3*c4*d2*e5-a4*c2*d3*e5+a2*c4*d3*e5  &
       +a3*c2*d4*e5-a2*c3*d4*e5)/det
  mati(1,3)=(a5*b4*d3*e2-a4*b5*d3*e2-a5*b3  &
       *d4*e2+a3*b5*d4*e2+a4*b3*d5*e2-a3*b4  &
       *d5*e2-a5*b4*d2*e3+a4*b5*d2*e3+a5*b2*d4  &
       *e3-a2*b5*d4*e3-a4*b2*d5*e3+a2*b4*d5*e3  &
       +a5*b3*d2*e4-a3*b5*d2*e4-a5*b2*d3*e4+a2  &
       *b5*d3*e4+a3*b2*d5*e4-a2*b3*d5*e4-a4*b3  &
       *d2*e5+a3*b4*d2*e5+a4*b2*d3*e5-a2*b4*d3  &
       *e5-a3*b2*d4*e5+a2*b3*d4*e5)/det
  mati(1,4)=(-a5*b4*c3*e2+a4*b5*c3*e2+a5*b3  &
       *c4*e2-a3*b5*c4*e2-a4*b3*c5*e2+a3*b4*c5  &
       *e2+a5*b4*c2*e3-a4*b5*c2*e3-a5*b2*c4*e3  &
       +a2*b5*c4*e3+a4*b2*c5*e3-a2*b4*c5*e3-a5  &
       *b3*c2*e4+a3*b5*c2*e4+a5*b2*c3*e4-a2*b5  &
       *c3*e4-a3*b2*c5*e4+a2*b3*c5*e4+a4*b3*c2  &
       *e5-a3*b4*c2*e5-a4*b2*c3*e5+a2*b4*c3*e5  &
       +a3*b2*c4*e5-a2*b3*c4*e5)/det
  mati(1,5)=(a5*b4*c3*d2-a4*b5*c3*d2-a5*b3  &
       *c4*d2+a3*b5*c4*d2+a4*b3*c5*d2-a3*b4*c5  &
       *d2-a5*b4*c2*d3+a4*b5*c2*d3+a5*b2*c4*d3  &
       -a2*b5*c4*d3-a4*b2*c5*d3+a2*b4*c5*d3+a5  &
       *b3*c2*d4-a3*b5*c2*d4-a5*b2*c3*d4+a2*b5  &
       *c3*d4+a3*b2*c5*d4-a2*b3*c5*d4-a4*b3*c2  &
       *d5+a3*b4*c2*d5+a4*b2*c3*d5-a2*b4*c3*d5  &
       -a3*b2*c4*d5+a2*b3*c4*d5)/det

  mati(2,1)=(-b5*c4*d3*e1+b4*c5*d3*e1+b5  &
       *c3*d4*e1-b3*c5*d4*e1-b4*c3*d5*e1+b3  &
       *c4*d5*e1+b5*c4*d1*e3-b4*c5*d1*e3-b5  &
       *c1*d4*e3+b1*c5*d4*e3+b4*c1*d5*e3-b1  &
       *c4*d5*e3-b5*c3*d1*e4+b3*c5*d1*e4+b5  &
       *c1*d3*e4-b1*c5*d3*e4-b3*c1*d5*e4+b1  &
       *c3*d5*e4+b4*c3*d1*e5-b3*c4*d1*e5-b4  &
       *c1*d3*e5+b1*c4*d3*e5+b3*c1*d4*e5-b1  &
       *c3*d4*e5)/det
  mati(2,2)=(a5*c4*d3*e1-a4*c5*d3*e1-a5  &
       *c3*d4*e1+a3*c5*d4*e1+a4*c3*d5*e1-a3  &
       *c4*d5*e1-a5*c4*d1*e3+a4*c5*d1*e3+a5  &
       *c1*d4*e3-a1*c5*d4*e3-a4*c1*d5*e3+a1  &
       *c4*d5*e3+a5*c3*d1*e4-a3*c5*d1*e4-a5  &
       *c1*d3*e4+a1*c5*d3*e4+a3*c1*d5*e4-a1  &
       *c3*d5*e4-a4*c3*d1*e5+a3*c4*d1*e5+a4  &
       *c1*d3*e5-a1*c4*d3*e5-a3*c1*d4*e5+a1  &
       *c3*d4*e5)/det
  mati(2,3)=(-a5*b4*d3*e1+a4*b5*d3*e1+a5  &
       *b3*d4*e1-a3*b5*d4*e1-a4*b3*d5*e1+a3  &
       *b4*d5*e1+a5*b4*d1*e3-a4*b5*d1*e3-a5  &
       *b1*d4*e3+a1*b5*d4*e3+a4*b1*d5*e3-a1  &
       *b4*d5*e3-a5*b3*d1*e4+a3*b5*d1*e4+a5  &
       *b1*d3*e4-a1*b5*d3*e4-a3*b1*d5*e4+a1  &
       *b3*d5*e4+a4*b3*d1*e5-a3*b4*d1*e5-a4  &
       *b1*d3*e5+a1*b4*d3*e5+a3*b1*d4*e5-a1  &
       *b3*d4*e5)/det
  mati(2,4)=(a5*b4*c3*e1-a4*b5*c3*e1-a5  &
       *b3*c4*e1+a3*b5*c4*e1+a4*b3*c5*e1-a3  &
       *b4*c5*e1-a5*b4*c1*e3+a4*b5*c1*e3+a5  &
       *b1*c4*e3-a1*b5*c4*e3-a4*b1*c5*e3+a1  &
       *b4*c5*e3+a5*b3*c1*e4-a3*b5*c1*e4-a5  &
       *b1*c3*e4+a1*b5*c3*e4+a3*b1*c5*e4-a1  &
       *b3*c5*e4-a4*b3*c1*e5+a3*b4*c1*e5+a4  &
       *b1*c3*e5-a1*b4*c3*e5-a3*b1*c4*e5+a1  &
       *b3*c4*e5)/det
  mati(2,5)=(-a5*b4*c3*d1+a4*b5*c3*d1+a5  &
       *b3*c4*d1-a3*b5*c4*d1-a4*b3*c5*d1+a3  &
       *b4*c5*d1+a5*b4*c1*d3-a4*b5*c1*d3-a5  &
       *b1*c4*d3+a1*b5*c4*d3+a4*b1*c5*d3-a1  &
       *b4*c5*d3-a5*b3*c1*d4+a3*b5*c1*d4+a5  &
       *b1*c3*d4-a1*b5*c3*d4-a3*b1*c5*d4+a1  &
       *b3*c5*d4+a4*b3*c1*d5-a3*b4*c1*d5-a4  &
       *b1*c3*d5+a1*b4*c3*d5+a3*b1*c4*d5-a1  &
       *b3*c4*d5)/det

  mati(3,1)=(b5*c4*d2*e1-b4*c5*d2*e1-b5  &
       *c2*d4*e1+b2*c5*d4*e1+b4*c2*d5*e1-b2  &
       *c4*d5*e1-b5*c4*d1*e2+b4*c5*d1*e2+b5  &
       *c1*d4*e2-b1*c5*d4*e2-b4*c1*d5*e2+b1  &
       *c4*d5*e2+b5*c2*d1*e4-b2*c5*d1*e4-b5  &
       *c1*d2*e4+b1*c5*d2*e4+b2*c1*d5*e4-b1  &
       *c2*d5*e4-b4*c2*d1*e5+b2*c4*d1*e5+b4  &
       *c1*d2*e5-b1*c4*d2*e5-b2*c1*d4*e5+b1  &
       *c2*d4*e5)/det
  mati(3,2)=(-a5*c4*d2*e1+a4*c5*d2*e1+a5  &
       *c2*d4*e1-a2*c5*d4*e1-a4*c2*d5*e1+a2  &
       *c4*d5*e1+a5*c4*d1*e2-a4*c5*d1*e2-a5  &
       *c1*d4*e2+a1*c5*d4*e2+a4*c1*d5*e2-a1  &
       *c4*d5*e2-a5*c2*d1*e4+a2*c5*d1*e4+a5  &
       *c1*d2*e4-a1*c5*d2*e4-a2*c1*d5*e4+a1  &
       *c2*d5*e4+a4*c2*d1*e5-a2*c4*d1*e5-a4  &
       *c1*d2*e5+a1*c4*d2*e5+a2*c1*d4*e5-a1  &
       *c2*d4*e5)/det
  mati(3,3)=(a5*b4*d2*e1-a4*b5*d2*e1-a5  &
       *b2*d4*e1+a2*b5*d4*e1+a4*b2*d5*e1-a2  &
       *b4*d5*e1-a5*b4*d1*e2+a4*b5*d1*e2+a5  &
       *b1*d4*e2-a1*b5*d4*e2-a4*b1*d5*e2+a1  &
       *b4*d5*e2+a5*b2*d1*e4-a2*b5*d1*e4-a5  &
       *b1*d2*e4+a1*b5*d2*e4+a2*b1*d5*e4-a1  &
       *b2*d5*e4-a4*b2*d1*e5+a2*b4*d1*e5+a4  &
       *b1*d2*e5-a1*b4*d2*e5-a2*b1*d4*e5+a1  &
       *b2*d4*e5)/det
  mati(3,4)=(-a5*b4*c2*e1+a4*b5*c2*e1+a5  &
       *b2*c4*e1-a2*b5*c4*e1-a4*b2*c5*e1+a2  &
       *b4*c5*e1+a5*b4*c1*e2-a4*b5*c1*e2-a5  &
       *b1*c4*e2+a1*b5*c4*e2+a4*b1*c5*e2-a1  &
       *b4*c5*e2-a5*b2*c1*e4+a2*b5*c1*e4+a5  &
       *b1*c2*e4-a1*b5*c2*e4-a2*b1*c5*e4+a1  &
       *b2*c5*e4+a4*b2*c1*e5-a2*b4*c1*e5-a4  &
       *b1*c2*e5+a1*b4*c2*e5+a2*b1*c4*e5-a1  &
       *b2*c4*e5)/det
  mati(3,5)=(a5*b4*c2*d1-a4*b5*c2*d1-a5  &
       *b2*c4*d1+a2*b5*c4*d1+a4*b2*c5*d1-a2  &
       *b4*c5*d1-a5*b4*c1*d2+a4*b5*c1*d2+a5  &
       *b1*c4*d2-a1*b5*c4*d2-a4*b1*c5*d2+a1  &
       *b4*c5*d2+a5*b2*c1*d4-a2*b5*c1*d4-a5  &
       *b1*c2*d4+a1*b5*c2*d4+a2*b1*c5*d4-a1  &
       *b2*c5*d4-a4*b2*c1*d5+a2*b4*c1*d5+a4  &
       *b1*c2*d5-a1*b4*c2*d5-a2*b1*c4*d5+a1  &
       *b2*c4*d5)/det

  mati(4,1)=(-b5*c3*d2*e1+b3*c5*d2*e1+b5  &
       *c2*d3*e1-b2*c5*d3*e1-b3*c2*d5*e1+b2  &
       *c3*d5*e1+b5*c3*d1*e2-b3*c5*d1*e2-b5  &
       *c1*d3*e2+b1*c5*d3*e2+b3*c1*d5*e2-b1  &
       *c3*d5*e2-b5*c2*d1*e3+b2*c5*d1*e3+b5  &
       *c1*d2*e3-b1*c5*d2*e3-b2*c1*d5*e3+b1  &
       *c2*d5*e3+b3*c2*d1*e5-b2*c3*d1*e5-b3  &
       *c1*d2*e5+b1*c3*d2*e5+b2*c1*d3*e5-b1  &
       *c2*d3*e5)/det
  mati(4,2)=(a5*c3*d2*e1-a3*c5*d2*e1-a5  &
       *c2*d3*e1+a2*c5*d3*e1+a3*c2*d5*e1-a2  &
       *c3*d5*e1-a5*c3*d1*e2+a3*c5*d1*e2+a5  &
       *c1*d3*e2-a1*c5*d3*e2-a3*c1*d5*e2+a1  &
       *c3*d5*e2+a5*c2*d1*e3-a2*c5*d1*e3-a5  &
       *c1*d2*e3+a1*c5*d2*e3+a2*c1*d5*e3-a1  &
       *c2*d5*e3-a3*c2*d1*e5+a2*c3*d1*e5+a3  &
       *c1*d2*e5-a1*c3*d2*e5-a2*c1*d3*e5+a1  &
       *c2*d3*e5)/det
  mati(4,3)=(-a5*b3*d2*e1+a3*b5*d2*e1+a5  &
       *b2*d3*e1-a2*b5*d3*e1-a3*b2*d5*e1+a2  &
       *b3*d5*e1+a5*b3*d1*e2-a3*b5*d1*e2-a5  &
       *b1*d3*e2+a1*b5*d3*e2+a3*b1*d5*e2-a1  &
       *b3*d5*e2-a5*b2*d1*e3+a2*b5*d1*e3+a5  &
       *b1*d2*e3-a1*b5*d2*e3-a2*b1*d5*e3+a1  &
       *b2*d5*e3+a3*b2*d1*e5-a2*b3*d1*e5-a3  &
       *b1*d2*e5+a1*b3*d2*e5+a2*b1*d3*e5-a1  &
       *b2*d3*e5)/det
  mati(4,4)=(a5*b3*c2*e1-a3*b5*c2*e1-a5  &
       *b2*c3*e1+a2*b5*c3*e1+a3*b2*c5*e1-a2  &
       *b3*c5*e1-a5*b3*c1*e2+a3*b5*c1*e2+a5  &
       *b1*c3*e2-a1*b5*c3*e2-a3*b1*c5*e2+a1  &
       *b3*c5*e2+a5*b2*c1*e3-a2*b5*c1*e3-a5  &
       *b1*c2*e3+a1*b5*c2*e3+a2*b1*c5*e3-a1  &
       *b2*c5*e3-a3*b2*c1*e5+a2*b3*c1*e5+a3  &
       *b1*c2*e5-a1*b3*c2*e5-a2*b1*c3*e5+a1  &
       *b2*c3*e5)/det
  mati(4,5)=(-a5*b3*c2*d1+a3*b5*c2*d1+a5  &
       *b2*c3*d1-a2*b5*c3*d1-a3*b2*c5*d1+a2  &
       *b3*c5*d1+a5*b3*c1*d2-a3*b5*c1*d2-a5  &
       *b1*c3*d2+a1*b5*c3*d2+a3*b1*c5*d2-a1  &
       *b3*c5*d2-a5*b2*c1*d3+a2*b5*c1*d3+a5  &
       *b1*c2*d3-a1*b5*c2*d3-a2*b1*c5*d3+a1  &
       *b2*c5*d3+a3*b2*c1*d5-a2*b3*c1*d5-a3  &
       *b1*c2*d5+a1*b3*c2*d5+a2*b1*c3*d5-a1  &
       *b2*c3*d5)/det

  mati(5,1)=(b4*c3*d2*e1-b3*c4*d2*e1-b4  &
       *c2*d3*e1+b2*c4*d3*e1+b3*c2*d4*e1-b2  &
       *c3*d4*e1-b4*c3*d1*e2+b3*c4*d1*e2+b4  &
       *c1*d3*e2-b1*c4*d3*e2-b3*c1*d4*e2+b1  &
       *c3*d4*e2+b4*c2*d1*e3-b2*c4*d1*e3-b4  &
       *c1*d2*e3+b1*c4*d2*e3+b2*c1*d4*e3-b1  &
       *c2*d4*e3-b3*c2*d1*e4+b2*c3*d1*e4+b3  &
       *c1*d2*e4-b1*c3*d2*e4-b2*c1*d3*e4+b1  &
       *c2*d3*e4)/det
  mati(5,2)=(-a4*c3*d2*e1+a3*c4*d2*e1+a4  &
       *c2*d3*e1-a2*c4*d3*e1-a3*c2*d4*e1+a2  &
       *c3*d4*e1+a4*c3*d1*e2-a3*c4*d1*e2-a4  &
       *c1*d3*e2+a1*c4*d3*e2+a3*c1*d4*e2-a1  &
       *c3*d4*e2-a4*c2*d1*e3+a2*c4*d1*e3+a4  &
       *c1*d2*e3-a1*c4*d2*e3-a2*c1*d4*e3+a1  &
       *c2*d4*e3+a3*c2*d1*e4-a2*c3*d1*e4-a3  &
       *c1*d2*e4+a1*c3*d2*e4+a2*c1*d3*e4-a1  &
       *c2*d3*e4)/det
  mati(5,3)=(a4*b3*d2*e1-a3*b4*d2*e1-a4  &
       *b2*d3*e1+a2*b4*d3*e1+a3*b2*d4*e1-a2  &
       *b3*d4*e1-a4*b3*d1*e2+a3*b4*d1*e2+a4  &
       *b1*d3*e2-a1*b4*d3*e2-a3*b1*d4*e2+a1  &
       *b3*d4*e2+a4*b2*d1*e3-a2*b4*d1*e3-a4  &
       *b1*d2*e3+a1*b4*d2*e3+a2*b1*d4*e3-a1  &
       *b2*d4*e3-a3*b2*d1*e4+a2*b3*d1*e4+a3  &
       *b1*d2*e4-a1*b3*d2*e4-a2*b1*d3*e4+a1  &
       *b2*d3*e4)/det
  mati(5,4)=(-a4*b3*c2*e1+a3*b4*c2*e1+a4  &
       *b2*c3*e1-a2*b4*c3*e1-a3*b2*c4*e1+a2  &
       *b3*c4*e1+a4*b3*c1*e2-a3*b4*c1*e2-a4  &
       *b1*c3*e2+a1*b4*c3*e2+a3*b1*c4*e2-a1  &
       *b3*c4*e2-a4*b2*c1*e3+a2*b4*c1*e3+a4  &
       *b1*c2*e3-a1*b4*c2*e3-a2*b1*c4*e3+a1  &
       *b2*c4*e3+a3*b2*c1*e4-a2*b3*c1*e4-a3  &
       *b1*c2*e4+a1*b3*c2*e4+a2*b1*c3*e4-a1  &
       *b2*c3*e4)/det
  mati(5,5)=(a4*b3*c2*d1-a3*b4*c2*d1-a4  &
       *b2*c3*d1+a2*b4*c3*d1+a3*b2*c4*d1-a2  &
       *b3*c4*d1-a4*b3*c1*d2+a3*b4*c1*d2+a4  &
       *b1*c3*d2-a1*b4*c3*d2-a3*b1*c4*d2+a1  &
       *b3*c4*d2+a4*b2*c1*d3-a2*b4*c1*d3-a4  &
       *b1*c2*d3+a1*b4*c2*d3+a2*b1*c4*d3-a1  &
       *b2*c4*d3-a3*b2*c1*d4+a2*b3*c1*d4+a3  &
       *b1*c2*d4-a1*b3*c2*d4-a2*b1*c3*d4+a1  &
       *b2*c3*d4)/det

  return
end subroutine matinv5
