program PorovnejSVD

! Porovnani vysledku z procedury svd v NM a svd z FMM
! Matice se zde nacita ze souboru Mat_A.dat
! 11.1.2008

use std_type
use LinEqu						! NM
use LinalgFMM, svdFMM=>svd		! FMM
use VstupVystup					! verze >= 2.2

implicit none

integer :: i,j,radku,sloupcu,iostav,ierr
real(DP),dimension(:,:),allocatable :: A,AF,U,UF,V,VF	! pro NM a FMM
real(DP),dimension(:),allocatable :: s,sF				! sing.cisla a 1D vstup
logical :: matUF, matVF					! pro FMM, true = matice UF/VF se pocita
character(len=40) :: vstupniA			! jmeno souboru se vstupni matici A

!call cti_cislo("Pocet radku matice A= ",radku)
!call cti_cislo("Pocet sloupcu matice A= ",sloupcu)
! Udaje jsou v prvnim radku vstupniho souboru, viz Mat_A.txt

call cti_text('Vstupni soubor = ',vstupniA)
open(unit=20,file=vstupniA,iostat=iostav,position='rewind',&
	action='read',status='old') 
if (iostav/=0) then
	print *,"IOSTAT= ",iostav
	goto 99
end if
read(unit=20,fmt=*)radku,sloupcu
!read(unit=20,fmt=*)sloupcu
print *,"matice ",radku,"x",sloupcu


allocate(A(radku,sloupcu))			! vstup pro NM, prepise ji U
allocate(AF(radku,sloupcu))			! vstup pro FMM, zachova se
allocate(U(radku,sloupcu))			! = prepsana matice A
allocate(UF(radku,radku))			! U pro FMM a matUF=true
allocate(V(sloupcu,sloupcu))		! V
allocate(VF(sloupcu,sloupcu))		! V pro FMM a matVF=true
allocate(s(sloupcu))				! vektor singularnich cisel z NM
allocate(sF(sloupcu))				! vektor singularnich cisel z FMM

do i=1,radku  		!cteni vstupu po radcich
	do j=1,sloupcu
		read(unit=20,fmt=*) AF(i,j)
	end do
end do
close(unit=20)

A=AF
print "(a,es12.3)","Minimalni prvek A = ",minval(A)
print "(a,es12.3)","Maximalni prvek A = ",maxval(A),char(10)

call svd(radku,sloupcu,A,V,s)

call svdFMM(AF,sF,.true.,UF,.true.,VF,ierr)
								! singularni cisla z svdFMM nejsou usporadana
call ShellSort(sF,.false.)		! Usporadat sestupne
print *,"ierr z procedury svdFMM =",ierr,char(10)


print *,"Singularni cisla  z NM, FMM a rozdily (NM-FMM)"//char(10)
print "(3a15)","NM ","FMM ","(NM-FMM)"//char(10)
do i=1,sloupcu
	print "(3es15.4)",s(i),sF(i),s(i)-sF(i)
end do


99	print *,char(10),"Konec = ENTER"
read *


!*****************************************************************************
contains
!*****************************************************************************


subroutine ShellSort(pole,vzestupne)
!
! Shelluv algoritmus pro usporadani pole
! Pouzije se pro usporadani singularnich cisel z svdFMM
! VSTUP:
!	pole	 	neusporadane 1D pole
!	vzestupne	OPTIONAL, .true./.false. usporada vzestupne(default)/sestupne
! VYSTUP:
!	pole		usporadane pole
!
	real(DP),dimension(:),intent(inout) :: pole
	logical,intent(in),optional :: vzestupne
	
	logical :: sVzestupne, flag
	integer :: ibeg,iend,j,k,ma,mb
	real(DP) :: w
	
	sVzestupne = .true.
	if(present(Vzestupne)) sVzestupne=Vzestupne
	
	ibeg = lbound(pole,1)
	iend = ubound(pole,1)
	ma = iend-ibeg+1
	mb = ma
	do
		mb = mb/2
		do
			flag = .true.
			do j=ibeg,(ibeg+ma-mb-1)

				k = j+mb
				if (sVzestupne) then
					if (pole(j)>pole(k)) then
						w = pole(j)
						pole(j) = pole(k)
						pole(k) = w
						flag = .false.
					end if
				else
					if (pole(j)<pole(k)) then
						w = pole(j)
						pole(j) = pole(k)
						pole(k) = w
						flag = .false.
					end if
				end if
			end do
			if (flag) exit
		end do
		if(mb<=1) exit
	end do
end subroutine ShellSort


end program PorovnejSVD
