Program T_GaussFMM

! Test Decomp + Solve z modulu LinAlgFMM

use std_type
use VstupVystup		! Verse >= 2.2
	! aliasy pro procedury decomp a solve z LinAlgFMM
use LinAlgFMM, FMMdecomp=>decomp, FMMsolve=>solve

implicit none

! Zadani presnosti pro realna cisla
! pro pohodlnou zamenu presnosti: FP=SP (single prec.), FP=DP (double prec.)
! procedury z modulu FMM je mozne volat s SP i DP
integer,parameter :: FP = DP

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

IF (VstupVystupver<"2.2") then
	print *,"Vyzaduje VstupVystup verse >= 2.2"
	goto 99
end if

print *,"Presnost: realna cisla ",FP,"byte",char(10)

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(A,FmtRetez)

!..............................................................................
call FMMDecomp(A,IndxPvt,JeReg,condA)
!..............................................................................

print *,char(10)//"Matice A po Decomp:"
FmtRetez = "(f12.3)"
call pis_matici(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)//"Odhad condA=",condA
print *,char(10)//"IndxPvt:"
FmtRetez = "(i3)"
call pis_vektor(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 FMMsolve(A,b,IndxPvt)
!..............................................................................
    
    print *, "Reseni je :"
    call pis_vektor(b,"(es15.6e3)")
end if

99 print *,"Pro ukonceni zadej ENTER "
read *


end program T_GaussFMM
