scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-23 04:43:38 +01:00
Code Issues Releases Wiki Activity
a264f8bb6f
qp_plugins_scemama/devel/svdwf/buildpsi_SVDit.irp.f
Abdallah Ammar 2fd2fcac5d svd save
2021-07-28 17:19:18 +02:00

678 lines
22 KiB
Fortran
Raw Blame History

program buildpsi_SVDit
implicit none
BEGIN_DOC
! perturbative approach to build psi_postsvd
END_DOC
read_wf = .True.
TOUCH read_wf
PROVIDE N_int
call run()
end
subroutine run
USE OMP_LIB
implicit none
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: degree, i_state
double precision :: h12
integer :: i, j, k, l, ii, jj
double precision :: norm_psi
double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)
double precision :: err0, err_tmp, e_tmp
double precision :: E0, E0pt2, ept2, E0_old, tol_energy
double precision :: ctmp, htmp
double precision, allocatable :: H0(:,:), Hdiag(:), Hkl(:,:)
double precision, allocatable :: coeff_psi_selected(:), coeff_psi_toselect(:)
integer :: n_FSVD, n_selected, n_toselect, it_svd, it_svd_max
integer, allocatable :: numalpha_selected(:), numbeta_selected(:)
integer, allocatable :: numalpha_toselect(:), numbeta_toselect(:)
integer, allocatable :: numalpha_tmp(:), numbeta_tmp(:)
integer(kind=8) :: W_tbeg, W_tend, W_tbeg_it, W_tend_it, W_ir
real(kind=8) :: W_tot_time, W_tot_time_it
real(kind=8) :: CPU_tbeg, CPU_tend, CPU_tbeg_it, CPU_tend_it
real(kind=8) :: CPU_tot_time, CPU_tot_time_it
real(kind=8) :: speedup, speedup_it
integer :: nb_taches
!$OMP PARALLEL
nb_taches = OMP_GET_NUM_THREADS()
!$OMP END PARALLEL
call CPU_TIME(CPU_tbeg)
call SYSTEM_CLOCK(COUNT=W_tbeg, COUNT_RATE=W_ir)
i_state = 1
det1(:,1) = psi_det_alpha_unique(:,1)
det2(:,1) = psi_det_alpha_unique(:,1)
call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int)
det1(:,2) = psi_det_beta_unique(:,1)
det2(:,2) = psi_det_beta_unique(:,1)
call get_excitation_degree(det1,det2,degree,N_int)
call i_H_j(det1, det2, N_int, h12)
! ---------------------------------------------------------------------------------------
! construct the initial CISD matrix
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
print *, ' CISD matrix:', n_det_alpha_unique,'x',n_det_beta_unique
print *, ' N det :', N_det
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
norm_psi = 0.d0
do k = 1, N_det
norm_psi = norm_psi + psi_bilinear_matrix_values(k,i_state) &
* psi_bilinear_matrix_values(k,i_state)
enddo
print *, ' initial norm = ', norm_psi
allocate( Aref(n_det_alpha_unique,n_det_beta_unique) )
Aref(:,:) = 0.d0
do k = 1, N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Aref(i,j) = psi_bilinear_matrix_values(k,i_state)
enddo
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! perform a Full SVD
allocate( Uref(n_det_alpha_unique,n_det_beta_unique) )
allocate( Dref(n_det_beta_unique) )
allocate( Vref(n_det_beta_unique,n_det_beta_unique) )
allocate( Vtref(n_det_beta_unique,n_det_beta_unique) )
call svd_s(Aref, size(Aref,1), Uref, size(Uref,1), Dref, Vtref, size(Vtref,1) &
, n_det_alpha_unique, n_det_beta_unique)
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
print *, ' --- First SVD: ok --- '
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
do l = 1, n_det_beta_unique
do i = 1, n_det_beta_unique
Vref(i,l) = Vtref(l,i)
enddo
enddo
deallocate( Vtref )
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! numerote vectors
! Full rank
n_FSVD = n_det_beta_unique*n_det_beta_unique
print*, ' Full psi space rank = ', n_FSVD
! Truncated rank
n_selected = 20
print*, ' initial psi space rank = ', n_selected
allocate( numalpha_selected(n_selected) , numbeta_selected(n_selected) )
do i = 1, n_selected
numalpha_selected(i) = i
numbeta_selected (i) = i
enddo
! check SVD error
err0 = 0.d0
do j = 1, n_det_beta_unique
do i = 1, n_det_alpha_unique
err_tmp = 0.d0
do l = 1, n_selected
ii = numalpha_selected(l)
jj = numbeta_selected (l)
err_tmp = err_tmp + Dref(l) * Uref(i,ii) * Vref(j,jj)
enddo
err_tmp = Aref(i,j) - err_tmp
err0 += err_tmp * err_tmp
enddo
enddo
print *, ' SVD err (%) = ', 100.d0 * dsqrt(err0/norm_psi)
deallocate( Aref )
! perturbative space rank
l = 3
k = 0
if( l.eq.1 ) then
n_toselect = 2*n_selected * ( n_det_beta_unique - n_selected )
allocate( numalpha_toselect(n_toselect) , numbeta_toselect(n_toselect) )
! nondiagonal blocs
do i = 1, n_selected
do j = n_selected+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
do j = 1, n_selected
do i = n_selected+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
elseif( l.eq.2 ) then
n_toselect = n_FSVD - n_selected*n_selected
allocate( numalpha_toselect(n_toselect) , numbeta_toselect(n_toselect) )
! nondiagonal blocs
do i = 1, n_selected
do j = n_selected+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
do j = 1, n_selected
do i = n_selected+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
! diagonal bloc
do i = n_selected+1, n_det_beta_unique
do j = n_selected+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
elseif( l.eq.3 ) then
n_toselect = n_FSVD - n_selected
allocate( numalpha_toselect(n_toselect) , numbeta_toselect(n_toselect) )
do i = 1, n_det_beta_unique
do j = 1, n_det_beta_unique
if( (i.eq.j).and.(i.le.n_selected)) then
cycle
else
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
endif
enddo
enddo
elseif( l.eq.4 ) then
n_toselect = n_FSVD - n_selected - (n_det_beta_unique-n_selected)**2
allocate( numalpha_toselect(n_toselect) , numbeta_toselect(n_toselect) )
! nondiagonal blocs
do i = 1, n_selected
do j = i+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
do j = 1, n_selected
do i = j+1, n_det_beta_unique
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
endif
if( k.ne.n_toselect ) then
print*, ' error in numeroting '
stop
endif
print*, ' perturbative psi space rank = ', n_toselect
! ---------------------------------------------------------------------------------------
!________________________________________________________________________________________________________
!
! loop over SVD iterations
!________________________________________________________________________________________________________
E0_old = 0.d0
tol_energy = 1.d0
it_svd = 0
it_svd_max = 100
do while( ( it_svd .lt. it_svd_max) .and. ( tol_energy .gt. 1d-8 ) )
call CPU_TIME(CPU_tbeg_it)
call SYSTEM_CLOCK(COUNT=W_tbeg_it, COUNT_RATE=W_ir)
it_svd = it_svd + 1
print*, '+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +'
print*, ' '
print*, ' iteration', it_svd
! ---------------------------------------------------------------------------------------
! calculate the energy
allocate( coeff_psi_selected(n_selected) )
! normalise | psi0 >
norm_psi = 0.d0
do i = 1, n_selected
norm_psi += Dref(i) * Dref(i)
enddo
norm_psi = 1.d0 / dsqrt(norm_psi)
do i = 1, n_selected
coeff_psi_selected(i) = Dref(i) * norm_psi
enddo
! H0(i,j) = < u_i v_j | H | u_i v_j >
print *, ''
print *, ''
print *, ''
print *, '-- Compute H --'
allocate( H0(n_selected,n_selected) )
call const_psiHpsi(n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0)
! avant SVD
! E0 = < psi_0 | H | psi_0 > / < psi_0 | psi_0 >
E0 = 0.d0
do i = 1, n_selected
ii = numalpha_selected(i)
htmp = 0.d0
do j = 1, n_selected
jj = numalpha_selected(j)
htmp = htmp + coeff_psi_selected(j) * H0(jj,ii)
enddo
E0 = E0 + htmp * coeff_psi_selected(i)
enddo
E0 = E0 + nuclear_repulsion
print *,' E0 (avant SVD) =', E0
deallocate( H0 )
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! nondiagonal elements
print *, ' --- Perturbation --- '
allocate( Hdiag(n_toselect), Hkl(n_selected,n_toselect) )
call const_Hdiag_Hkl(n_selected, n_toselect, Uref, Vref &
, numalpha_selected, numbeta_selected, numalpha_toselect, numbeta_toselect, Hdiag, Hkl)
! evaluate the coefficients for all the vectors
allocate( coeff_psi_toselect(n_toselect) )
ept2 = 0.d0
do ii = 1, n_toselect
ctmp = 0.d0
do l = 1, n_selected
ctmp += coeff_psi_selected(l) * Hkl(l,ii)
enddo
coeff_psi_toselect(ii) = ctmp / ( E0 - (Hdiag(ii)+nuclear_repulsion) )
ept2 += ctmp * ctmp / ( E0 - (Hdiag(ii)+nuclear_repulsion) )
enddo
E0pt2 = E0 + ept2
deallocate( Hdiag, Hkl)
print *, ' perturb energy = ', E0pt2, ept2
print*, ' delta E0 = ', E0pt2 - E0_old
tol_energy = 100.d0 * dabs(E0pt2-E0_old)/dabs(E0pt2)
E0_old = E0pt2
! normalize the new psi and perform a new SVD
norm_psi = 0.d0
do l = 1, n_toselect
norm_psi = norm_psi + coeff_psi_toselect(l)*coeff_psi_toselect(l)
enddo
norm_psi = norm_psi + 1.d0
norm_psi = 1.d0 / dsqrt(norm_psi)
do i = 1, n_toselect
coeff_psi_toselect(i) = coeff_psi_toselect(i) * norm_psi
enddo
do i = 1, n_selected
coeff_psi_selected(i) = coeff_psi_selected(i) * norm_psi
enddo
print *, ' --- SVD --- '
call perform_newSVD(n_selected, n_toselect, numalpha_selected, numbeta_selected &
, numalpha_toselect, numbeta_toselect, coeff_psi_selected, coeff_psi_toselect &
, Uref, Vref, Dref )
! ---------------------------------------------------------------------------------------
deallocate( coeff_psi_toselect )
deallocate( coeff_psi_selected )
write(55,'(i5,4x,4(f22.15,2x))') it_svd, E0, E0pt2
call CPU_TIME(CPU_tend_it)
call SYSTEM_CLOCK(COUNT=W_tend_it, COUNT_RATE=W_ir)
CPU_tot_time_it = CPU_tend_it - CPU_tbeg_it
W_tot_time_it = real(W_tend_it-W_tbeg_it, kind=8) / real(W_ir, kind=8)
speedup_it = CPU_tot_time_it / W_tot_time_it
print '(//, 3X, "elapsed time = ", 1PE10.3, " min.", /, &
& 3X, "CPU time = ", 1PE10.3, " min.", /, &
& 3X, "speed up = ", 1PE10.3,//)', W_tot_time_it/60.d0, CPU_tot_time_it/60.d0, speedup_it
!print*, '+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +'
!print*, ' '
end do
!________________________________________________________________________________________________________
!________________________________________________________________________________________________________
deallocate( Uref, Vref, Dref )
call CPU_TIME(CPU_tend)
call SYSTEM_CLOCK(COUNT=W_tend, COUNT_RATE=W_ir)
CPU_tot_time = CPU_tend - CPU_tbeg
W_tot_time = real(W_tend - W_tbeg, kind=8) / real(W_ir, kind=8)
speedup = CPU_tot_time / W_tot_time
print *,' ___________________________________________________________________'
print '(//,3X,"Execution avec ",i2," threads")',nb_taches
print *,' ___________________________________________________________________'
print '(//, 3X, "elapsed time = ", 1PE10.3, " min.", /, &
& 3X, "CPU time = ", 1PE10.3, " min.", /, &
& 3X, "speed up = ", 1PE10.3 ,//)', W_tot_time/60.d0, CPU_tot_time/60.d0, speedup
end
subroutine const_psiHpsi(n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0)
USE OMP_LIB
implicit none
integer, intent(in) :: n_selected
integer, intent(in) :: numalpha_selected(n_selected), numbeta_selected(n_selected)
double precision, intent(in) :: Uref(n_det_alpha_unique,n_det_beta_unique)
double precision, intent(in) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(out) :: H0(n_selected,n_selected)
integer(bit_kind) :: det1(N_int,2)
integer(bit_kind) :: det2(N_int,2)
integer :: degree
integer :: i, j, k, l
integer :: ii0, jj0, ii, jj, n, m, np, mp
double precision :: h12, x
H0(:,:) = 0.d0
do i = 1, n_det_alpha_unique
det1(:,1) = psi_det_alpha_unique(:,i)
do k = 1, n_det_alpha_unique
det2(:,1) = psi_det_alpha_unique(:,k)
call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int)
if (degree .gt. 2) then
cycle
endif
do j = 1, n_det_beta_unique
det1(:,2) = psi_det_beta_unique(:,j)
do l = 1, n_det_beta_unique
det2(:,2) = psi_det_beta_unique(:,l)
call get_excitation_degree(det1,det2,degree,N_int)
if (degree .gt. 2) then
cycle
endif
! !!!
call i_H_j(det1, det2, N_int, h12)
! !!!
! ~~~ H0 ~~~
do n = 1, n_selected
ii0 = numalpha_selected(n)
jj0 = numbeta_selected (n)
x = Uref(k,ii0) * Vref(l,jj0) * h12
do m = 1, n_selected
ii = numalpha_selected(m)
jj = numbeta_selected (m)
H0(m,n) += Uref(i,ii) * Vref(j,jj) * x
enddo
enddo
! ~~~ ~~~~~~ ~~~
enddo
enddo
enddo
enddo
end subroutine const_psiHpsi
subroutine const_Hdiag_Hkl(n_selected, n_toselect, Uref, Vref &
, numalpha_selected, numbeta_selected, numalpha_toselect, numbeta_toselect, Hdiag, Hkl)
USE OMP_LIB
implicit none
integer, intent(in) :: n_selected,n_toselect
integer, intent(in) :: numalpha_selected(n_selected), numbeta_selected(n_selected)
integer, intent(in) :: numalpha_toselect(n_toselect), numbeta_toselect(n_toselect)
double precision, intent(in) :: Uref(n_det_alpha_unique,n_det_beta_unique)
double precision, intent(in) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(out) :: Hdiag(n_toselect), Hkl(n_selected,n_toselect)
integer(bit_kind) :: det1(N_int,2)
integer(bit_kind) :: det2(N_int,2)
integer :: degree
integer :: i, j, k, l
integer :: ii0, jj0, ii, jj, n, m, np, mp
double precision :: h12, y
double precision, allocatable :: Hdiag_tmp(:), Hkl_tmp(:,:)
Hdiag(:) = 0.d0
Hkl(:,:) = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(n,ii0,jj0,y,m,ii,jj,i,j,k,l,h12,det1,det2,Hdiag_tmp,Hkl_tmp,degree) &
!$OMP SHARED(n_det_alpha_unique,n_det_beta_unique,psi_det_alpha_unique,psi_det_beta_unique, &
!$OMP N_int,n_selected,n_toselect,Uref,Vref,numalpha_toselect,numbeta_toselect, &
!$OMP numalpha_selected, numbeta_selected,Hkl,Hdiag )
allocate( Hdiag_tmp(n_toselect), Hkl_tmp(n_selected,n_toselect) )
Hdiag_tmp(:) = 0.d0
Hkl_tmp(:,:) = 0.d0
!$OMP DO COLLAPSE(2) SCHEDULE(DYNAMIC,8)
do i = 1, n_det_alpha_unique
do k = 1, n_det_alpha_unique
det1(:,1) = psi_det_alpha_unique(:,i)
det2(:,1) = psi_det_alpha_unique(:,k)
! !!!
! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int)
if (degree .gt. 2) then
cycle
endif
! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
! !!!
do j = 1, n_det_beta_unique
det1(:,2) = psi_det_beta_unique(:,j)
do l = 1, n_det_beta_unique
det2(:,2) = psi_det_beta_unique(:,l)
! !!!
! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
call get_excitation_degree(det1,det2,degree,N_int)
if (degree .gt. 2) then
cycle
endif
! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
! !!!
call i_H_j(det1, det2, N_int, h12)
! ~ ~ ~ H ~ ~ ~
do n = 1, n_toselect
ii0 = numalpha_toselect(n)
jj0 = numbeta_toselect (n)
y = Uref(k,ii0) * Vref(l,jj0) * h12
! Hdiag
Hdiag_tmp(n) += Uref(i,ii0) * Vref(j,jj0) * y
do m = 1, n_selected
ii = numalpha_selected(m)
jj = numbeta_selected (m)
! Hkl
Hkl_tmp(m,n) += Uref(i,ii) * Vref(j,jj) * y
enddo
enddo
! ~ ~ ~ ! ! ! ~ ~ ~
enddo
enddo
! !!!
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
do n = 1, n_toselect
Hdiag(n) += Hdiag_tmp(n)
do m = 1, n_selected
Hkl(m,n) += Hkl_tmp(m,n)
enddo
enddo
!$OMP END CRITICAL
deallocate( Hdiag_tmp,Hkl_tmp )
!$OMP END PARALLEL
end subroutine const_Hdiag_Hkl
subroutine perform_newSVD(n_selected, n_toselect, numalpha_selected, numbeta_selected &
, numalpha_toselect, numbeta_toselect, coeff_psi_selected, coeff_psi_toselect &
, Uref, Vref, Dref )
USE OMP_LIB
implicit none
integer, intent(in) :: n_selected, n_toselect
integer, intent(in) :: numalpha_selected(n_selected), numbeta_selected(n_selected)
integer, intent(in) :: numalpha_toselect(n_toselect), numbeta_toselect(n_toselect)
double precision, intent(in) :: coeff_psi_selected(n_selected), coeff_psi_toselect(n_toselect)
double precision, intent(inout) :: Uref(n_det_alpha_unique,n_det_beta_unique)
double precision, intent(inout) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(inout) :: Dref(n_det_beta_unique)
integer :: mm, nn, i, j, ii0, ii, l, jj
double precision :: err0, err_norm, err_tmp
double precision, allocatable :: S_mat(:,:), SxVt(:,:)
double precision, allocatable :: U_svd(:,:), V_svd(:,:)
double precision, allocatable :: U_newsvd(:,:), V_newsvd(:,:), Vt_newsvd(:,:), D_newsvd(:), A_newsvd(:,:)
mm = n_det_alpha_unique
nn = n_det_beta_unique
allocate( U_svd(n_det_alpha_unique,n_det_beta_unique) )
allocate( V_svd(n_det_beta_unique ,n_det_beta_unique) )
allocate( S_mat(n_det_beta_unique ,n_det_beta_unique) )
U_svd(:,:) = Uref(:,:)
V_svd(:,:) = Vref(:,:)
S_mat(:,:) = 0.d0
do l = 1, n_selected
ii = numalpha_selected(l)
jj = numbeta_selected (l)
S_mat(ii,jj) = coeff_psi_selected(l)
enddo
do l = 1, n_toselect
ii = numalpha_toselect(l)
jj = numbeta_toselect (l)
S_mat(ii,jj) = coeff_psi_toselect(l)
enddo
! construct the new matrix: U_svd x S_mat x transpose(V_svd)
! (NaxNb) (NbxNb) transpose(NbxNb)
! first compute S_mat x transpose(V_svd)
allocate( SxVt(nn,nn) )
call dgemm( 'N', 'T', nn, nn, nn, 1.d0 &
, S_mat , size(S_mat,1) &
, V_svd , size(V_svd,1) &
, 0.d0, SxVt, size(SxVt ,1) )
! then compute U_svd x SxVt
allocate( A_newsvd(mm,nn) )
call dgemm( 'N', 'N', mm, nn, nn, 1.d0 &
, U_svd , size(U_svd ,1) &
, SxVt , size(SxVt ,1) &
, 0.d0, A_newsvd, size(A_newsvd,1) )
deallocate( SxVt )
! perform new SVD
allocate( U_newsvd(mm,nn), Vt_newsvd(nn,nn), D_newsvd(nn) )
call svd_s( A_newsvd, size(A_newsvd,1), U_newsvd, size(U_newsvd,1), D_newsvd &
, Vt_newsvd, size(Vt_newsvd,1), mm, nn)
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
print *, ' +++ new SVD is performed +++ '
allocate( V_newsvd(nn,nn) )
do l = 1, nn
do j = 1, nn
V_newsvd(j,l) = Vt_newsvd(l,j)
enddo
enddo
! check SVD error
err0 = 0.d0
err_norm = 0.d0
do j = 1, nn
do i = 1, mm
err_tmp = 0.d0
do l = 1, mm
err_tmp = err_tmp + D_newsvd(l) * U_newsvd(i,l) * V_newsvd(j,l)
enddo
err_tmp = A_newsvd(i,j) - err_tmp
err0 += err_tmp * err_tmp
err_norm += A_newsvd(i,j) * A_newsvd(i,j)
enddo
enddo
print *, ' SVD err (%) = ', 100.d0 * dsqrt(err0/err_norm)
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
do l = 1, nn
Dref(l) = D_newsvd(l)
Uref(:,l) = U_newsvd(:,l)
Vref(:,l) = V_newsvd(:,l)
enddo
deallocate( U_newsvd, V_newsvd, Vt_newsvd, D_newsvd, A_newsvd )
return
end subroutine perform_newSVD
Reference in New Issue View Git Blame Copy Permalink