program test_matinv
use fvn_linear
implicit none

integer, parameter :: n=3
complex(kind=dp_kind),dimension(n,n) :: m1,m2,res
real(kind=dp_kind),dimension(n,n) :: rtmp,itmp
character(len=80) :: fmcmplx
integer :: i

fmcmplx='(3("(",f8.5,",",f8.5,")  "))'
! initialize pseudo random generator
!call init_random_seed()
! fill real and imaginary part
call random_number(rtmp)
call random_number(itmp)
! create the complex matrix (fvn_i is defined in the fvn module)
m1=rtmp+fvn_i*itmp
write(*,*) "Matrix M1"
do i=1,n
    write(*,fmcmplx) m1(i,:)
end do

call random_number(rtmp)
call random_number(itmp)
m2=rtmp+fvn_i*itmp
write(*,*)
write(*,*) "Matrix M2"
do i=1,n
    write(*,fmcmplx) m2(i,:)
end do


write(*,*)
write(*,*) "M1.x.M2"
res=m1.x.m2
call write_res()
write(*,*)
write(*,*) "M1.x.M2 standard"
res=matmul(m1,m2)
call write_res()

write(*,*)
write(*,*) ".i.M1"
res=.i.m1
call write_res()
write(*,*)
write(*,*) ".i.M1 standard"
call fvn_matinv(3,m1,res)
call write_res()

write(*,*)
write(*,*) "M1.ix.M2"
res=m1.ix.m2
call write_res()
write(*,*)
write(*,*) "M1.ix.M2 standard"
call fvn_matinv(3,m1,res)
res=matmul(res,m2)
call write_res()

write(*,*)
write(*,*) "M1.xi.M2"
res=m1.xi.m2
call write_res()
write(*,*)
write(*,*) "M1.xi.M2 standard"
res=m1.xi.m2
call fvn_matinv(3,m2,res)
res=matmul(m1,res)
call write_res()

write(*,*)
write(*,*) ".t.M1"
res=.t.m1
call write_res()
write(*,*)
write(*,*) ".t.M1 standard"
res=transpose(m1)
call write_res()

write(*,*)
write(*,*) "M1.tx.M2"
res=m1.tx.m2
call write_res()
write(*,*)
write(*,*) "M1.tx.M2 standard"
res=matmul(transpose(m1),m2)
call write_res()

write(*,*)
write(*,*) "M1.xt.M2"
res=m1.xt.m2
call write_res()
write(*,*)
write(*,*) "M1.xt.M2 standard"
res=matmul(m1,transpose(m2))
call write_res()

write(*,*)
write(*,*) ".h.M1"
res=.h.m1
call write_res()
write(*,*)
write(*,*) ".h.M1 standard"
res=transpose(conjg(m1))
call write_res()

write(*,*)
write(*,*) "M1.hx.M2"
res=m1.hx.m2
call write_res()
write(*,*)
write(*,*) "M1.hx.M2 standard"
res=matmul(transpose(conjg(m1)),m2)
call write_res()

write(*,*)
write(*,*) "M1.xh.M2"
res=m1.xh.m2
call write_res()
write(*,*)
write(*,*) "M1.xh.M2 standard"
res=matmul(m1,transpose(conjg(m2)))
call write_res()

contains
subroutine write_res()
do i=1,n
    write(*,fmcmplx) res(i,:)
end do
end subroutine
end program