Program T_svdFMM

! Test procedury svd z LinAlgFMM
! Procedura svd z svdFMM pouzivava JINY ALGORITMUS (Golub and Reinsch)
! nez procedura svd z modulu LinEqu
! Porovnejte ortogonalni matice U,V z obou procedur (mohou byt odlisne)
! a singularni cisla (musi byt v mezich zaokrouhlovacich chyb shodna 
! Uprava pro VV ver 2.2


use std_type
use VstupVystup
! zavede alias FMMsvd pro proceduru svd z modulu LinAlgFMM	
use LinAlgFMM,FMMsvd=>svd

implicit none

! pro snadnou zmenu presnosti realnych cisel (FP = SP | DP)
integer,parameter :: FP = DP

integer :: radku,sloupcu
real(FP),dimension(:,:),allocatable :: A,U,V,sA
character(len=20) :: FmtRetez
real(FP),dimension(:),allocatable :: w
integer :: ierr

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 radku = ",radku)
call cti_cislo("Pocet sloupcu = ",sloupcu)
allocate(A(radku,sloupcu))
Allocate(U(radku,radku))
Allocate(V(sloupcu,sloupcu))
allocate(w(sloupcu))

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

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

! ............................................................................
call FMMsvd	(A,w,.true.,U,.true.,V,ierr)
!.............................................................................


print *,char(10)//"Singularni cisla z Svd:"
FmtRetez = "(es12.4e2)"
call pis_vektor(w,FmtRetez,.true.)

print "(a,e12.3)",char(10)//"Cislo podminenosti S_max/S_min=",maxval(w)/minval(w)

print *,char(10)//"Matice U z Svd:"
FmtRetez = "(f12.4)"
call pis_matici(U,FmtRetez)

print *,char(10)//"Matice V z Svd:"
FmtRetez = "(f12.4)"
call pis_matici(V,FmtRetez)

Print *,char(10),"Kontrola ortogonality U: U_tran*U=I ?"
call pis_matici(MatMul(transpose(U),U),FmtRetez)

Print *,char(10),"Kontrola ortogonality V: V_tran*V=I ?"
call pis_matici(MatMul(transpose(V),V),FmtRetez)


99 print *,char(10),"Pro ukonceni zadej ENTER "
read *

end program T_svdFMM
