qp_plugins_scemama/devel/svdwf/buildpsi_SVDit.irp.f

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