Program T_GaussFMM

! Test Decomp + Solve z LinAlgFMM

use std_type
use VstupVystup
use LinalgFMM

implicit none

integer :: dimA,i
real(kind=SP) :: condA,det
real(Kind=SP),dimension(:,:),allocatable :: A
integer,dimension(:),allocatable :: IndxPvt
Integer :: JeReg
character(len=20) :: FmtRetez
real(kind=SP),dimension(:),allocatable :: b

call cti_cislo("Pocet rovnic = ",dimA)
allocate(A(dimA,dimA))
allocate(IndxPvt(dimA))

print *, "Zadej matici A po radcich"
call cti_matici(dimA,dimA,A)

print *
Print *,"Zadana matice A:"
FmtRetez = "(f10.3)"
call pis_matici(dimA,dimA,A,FmtRetez)

call Decomp(A,IndxPvt,JeReg,condA)

print *,char(10)//"Matice A po Decomp:"
FmtRetez = "(f12.3)"
call pis_matici(dimA,dimA,A,FmtRetez)
if (JeReg==0) then
	print *,char(10)//"Matice je regularni"
else
	print *,char(10)//"Matice je singularni"
end if
print "(a,es12.3e3)",char(10)//"condA=",condA
print *,char(10)//"IndxPvt:"
FmtRetez = "(i3)"
call pis_vektor(dimA,IndxPvt,FmtRetez,.true.)
det = 1
do i = 1,dimA
    det = det*A(i,i)
end do    
print "(a,es12.3e3)",char(10)//"Determinant=",IndxPvt(dimA)*det

if (JeReg==0) THEN
    allocate(b(dimA))
    print *, char(10)//"Zadej vektor prave strany: "
    call cti_vektor(dimA,b)
    
    call Solve(A,b,IndxPvt)
    
    print *, "Reseni je :"
    call pis_vektor(dimA,b,"(es15.6e3)")
end if

print *,"Pro ukonceni zadej ENTER "
read *


end program T_GaussFMM
