DecompositionOps

Fortran module DecompositionOps: May 2017 (dj)

Containing decompositions for tensors.

Authors

• 4. Jaschke

• 13. 12. Wall

Details

The following subroutines / functions are defined in the module as interfaces.

Procedure
qr
rq
split

DecompositionOps_f90.eigsvd_tensor()

fortran-subroutine - November 2016 (dj) Split a tensor into two tensors using a SVD.

Arguments

UtTYPE(tensor), inout

Contains on exit the left unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

LamTYPE(tensor), inout

Contains on exit the singular values. Intent in to keep dimension etc.

VtTYPE(tensor), inout

Contains on exit the right unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

TensTYPE(tensor), inout

Contains the tensor to be decomposed. Destroyed on exit.

idxlINTEGER(*), in

Indices of the legs going into the left tensor.

idxrINTEGER(*), in

Indices of the legs going into the right tensor.

multlrINTEGER, OPTIONAL, in

Multiply singular values to the left (< 0) or right (> 0). Default to no multiplication (0), but no valid option for eigsvd.

truncREAL, OPTIONAL, in

Keep infidelity below trunc (infidelity is sum of squared discarded singular values for MPS). Default does not truncate.

ncutINTEGER, OPTIONAL, in

Maximal bond dimension / number of singular values. Default is keeping all singular values.

errREAL, OPTIONAL, out

Sum over the square of the discarded singular value for MPS.

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor before returning it. Default to no permutation.

renormCHARACTER, OPTIONAL, in

Renormalization of the singular values. ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

Source Code

show / hide f90 code

  subroutine eigsvd_tensor(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, ncut, &
                                err, lpermout, rpermout, renorm, errst)
      type(tensor), intent(inout) :: Ut, Vt
      type(tensor), intent(inout) :: Lam
      type(tensor), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(in), optional :: multlr
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      character, intent(in), optional :: renorm
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! number of singular values kept
      integer :: nkeep

      ! number of legs going to the left / right
      integer :: nl, nr

      ! size of the matrix to be decomposed
      integer :: a_row, a_col

      ! first singular value kept
      integer :: first

      ! renormalization due to truncated singular values
      real(KIND=rKind) :: renormval

      ! all singular values
      type(tensor) :: Allsing

      ! temporary unitary matrix before truncation
      real(KIND=rKind), dimension(:, :), allocatable :: utmp

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('eigsvd_tensor: permute_split failed', &
      !              errst=errst)) return

      ! Run the SVD as eigenvalue decomposition (1st part)
      ! ---------------------------------------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      select case(multlr)
          case(1:)
              allocate(utmp(a_row, a_row))
              call create(allsing, [a_row])
              call eigsvd_lu_step1_real(utmp, Allsing%elem, Tens%elem, &
                                               a_row, a_col, errst=errst)
              !if(prop_error('eigsvd_tensor: eigsvd_lu_step1_'//&
              !              'real failed.', 'DecompositionOps_include.f90:159', &
              !              errst=errst)) return
          case(:-1)
              allocate(utmp(a_col, a_col))
              call create(Allsing, [a_col])
              call eigsvd_ru_step1_real(utmp, Allsing%elem, Tens%elem, &
                                               a_row, a_col, errst=errst)
              !if(prop_error('eigsvd_tensor: eigsvd_ru_step1_'//&
              !              'real failed.', 'DecompositionOps_include.f90:167', &
              !              errst=errst)) return
          case default
              errst = raise_error('eigsvd_tensor: multlr == 0!', &
                                  99, 'DecompositionOps_include.f90:171', errst=errst)
              return
      end select

      ! Decide on truncation
      ! --------------------

      call get_chi_mps(Allsing, 'E', nkeep, trunc, ncut, err, &
                       renormval, errst=errst)
      !if(prop_error('eigsvd_tensor: get_chi_mps failed.', &
      !              'DecompositionOps_include.f90:181', errst=errst)) return

      if(present(renorm)) then
          select case(renorm)
          case('M')
              renormval = 1.0_rKind / sqrt(renormval)
          case default
              renormval = 1.0_rKind
          end select
      else
          renormval = 1.0_rKind
      end if

      call create(Lam, [nkeep])
      first = size(Allsing%elem) - nkeep + 1
      Lam%elem(1:nkeep) = Allsing%elem(first:)
      call destroy(Allsing)

      ! Run the SVD as eigenvalue decomposition (2nd part)
      ! ---------------------------------------

      select case(multlr)
          case(1:)
              call create(Vt, [nkeep, Tens%dl(nl + 1:)])
              call eigsvd_lu_step2_real(utmp, Vt%elem, Tens%elem, a_row, &
                                               a_col, nkeep, renorm=renormval, &
                                               errst=errst)
              !if(prop_error('eigsvd_tensor: eigsvd_lu_step2_'//&
              !              'real failed.', 'DecompositionOps_include.f90:209', &
              !              errst=errst)) return

              call create(Ut, [Tens%dl(:nl), nkeep])
              Ut%elem = reshape(utmp(:, first:), [a_row * nkeep])
          case(:-1)
              call create(Ut, [Tens%dl(:nl), nkeep])
              call eigsvd_ru_step2_real(utmp, Ut%elem, Tens%elem, a_row, &
                                               a_col, nkeep, renorm=renormval, &
                                               errst=errst)
              !if(prop_error('eigsvd_tensor: eigsvd_ru_step2_'//&
              !              'real failed.', 'DecompositionOps_include.f90:220', &
              !              errst=errst)) return

              call create(Vt, [nkeep, Tens%dl(nl + 1:)])
              Vt%elem = reshape(utmp(first:, :), [nkeep * a_col])
      end select

      ! Permutation on the output
      ! -------------------------

      if(present(lpermout)) then
          call transposed(Ut, lpermout, doperm=.true., errst=errst)
          !if(prop_error('eigsvd_tensor: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Vt, rpermout, doperm=.true., errst=errst)
          !if(prop_error('eigsvd_tensor: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine eigsvd_tensor

DecompositionOps_f90.eigsvd_tensorc()

fortran-subroutine - November 2016 (dj) Split a tensor into two tensors using a SVD.

Arguments

UtTYPE(tensorc), inout

Contains on exit the left unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

LamTYPE(tensor), inout

Contains on exit the singular values. Intent in to keep dimension etc.

VtTYPE(tensorc), inout

Contains on exit the right unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

TensTYPE(tensorc), inout

Contains the tensor to be decomposed. Destroyed on exit.

idxlINTEGER(*), in

Indices of the legs going into the left tensor.

idxrINTEGER(*), in

Indices of the legs going into the right tensor.

multlrINTEGER, OPTIONAL, in

Multiply singular values to the left (< 0) or right (> 0). Default to no multiplication (0), but no valid option for eigsvd.

truncREAL, OPTIONAL, in

Keep infidelity below trunc (infidelity is sum of squared discarded singular values for MPS). Default does not truncate.

ncutINTEGER, OPTIONAL, in

Maximal bond dimension / number of singular values. Default is keeping all singular values.

errREAL, OPTIONAL, out

Sum over the square of the discarded singular value for MPS.

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor before returning it. Default to no permutation.

renormCHARACTER, OPTIONAL, in

Renormalization of the singular values. ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

Source Code

show / hide f90 code

  subroutine eigsvd_tensorc(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, ncut, &
                                err, lpermout, rpermout, renorm, errst)
      type(tensorc), intent(inout) :: Ut, Vt
      type(tensor), intent(inout) :: Lam
      type(tensorc), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(in), optional :: multlr
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      character, intent(in), optional :: renorm
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! number of singular values kept
      integer :: nkeep

      ! number of legs going to the left / right
      integer :: nl, nr

      ! size of the matrix to be decomposed
      integer :: a_row, a_col

      ! first singular value kept
      integer :: first

      ! renormalization due to truncated singular values
      real(KIND=rKind) :: renormval

      ! all singular values
      type(tensor) :: Allsing

      ! temporary unitary matrix before truncation
      complex(KIND=rKind), dimension(:, :), allocatable :: utmp

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('eigsvd_tensorc: permute_split failed', &
      !              errst=errst)) return

      ! Run the SVD as eigenvalue decomposition (1st part)
      ! ---------------------------------------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      select case(multlr)
          case(1:)
              allocate(utmp(a_row, a_row))
              call create(allsing, [a_row])
              call eigsvd_lu_step1_complex(utmp, Allsing%elem, Tens%elem, &
                                               a_row, a_col, errst=errst)
              !if(prop_error('eigsvd_tensorc: eigsvd_lu_step1_'//&
              !              'complex failed.', 'DecompositionOps_include.f90:159', &
              !              errst=errst)) return
          case(:-1)
              allocate(utmp(a_col, a_col))
              call create(Allsing, [a_col])
              call eigsvd_ru_step1_complex(utmp, Allsing%elem, Tens%elem, &
                                               a_row, a_col, errst=errst)
              !if(prop_error('eigsvd_tensorc: eigsvd_ru_step1_'//&
              !              'complex failed.', 'DecompositionOps_include.f90:167', &
              !              errst=errst)) return
          case default
              errst = raise_error('eigsvd_tensorc: multlr == 0!', &
                                  99, 'DecompositionOps_include.f90:171', errst=errst)
              return
      end select

      ! Decide on truncation
      ! --------------------

      call get_chi_mps(Allsing, 'E', nkeep, trunc, ncut, err, &
                       renormval, errst=errst)
      !if(prop_error('eigsvd_tensorc: get_chi_mps failed.', &
      !              'DecompositionOps_include.f90:181', errst=errst)) return

      if(present(renorm)) then
          select case(renorm)
          case('M')
              renormval = 1.0_rKind / sqrt(renormval)
          case default
              renormval = 1.0_rKind
          end select
      else
          renormval = 1.0_rKind
      end if

      call create(Lam, [nkeep])
      first = size(Allsing%elem) - nkeep + 1
      Lam%elem(1:nkeep) = Allsing%elem(first:)
      call destroy(Allsing)

      ! Run the SVD as eigenvalue decomposition (2nd part)
      ! ---------------------------------------

      select case(multlr)
          case(1:)
              call create(Vt, [nkeep, Tens%dl(nl + 1:)])
              call eigsvd_lu_step2_complex(utmp, Vt%elem, Tens%elem, a_row, &
                                               a_col, nkeep, renorm=renormval, &
                                               errst=errst)
              !if(prop_error('eigsvd_tensorc: eigsvd_lu_step2_'//&
              !              'complex failed.', 'DecompositionOps_include.f90:209', &
              !              errst=errst)) return

              call create(Ut, [Tens%dl(:nl), nkeep])
              Ut%elem = reshape(utmp(:, first:), [a_row * nkeep])
          case(:-1)
              call create(Ut, [Tens%dl(:nl), nkeep])
              call eigsvd_ru_step2_complex(utmp, Ut%elem, Tens%elem, a_row, &
                                               a_col, nkeep, renorm=renormval, &
                                               errst=errst)
              !if(prop_error('eigsvd_tensorc: eigsvd_ru_step2_'//&
              !              'complex failed.', 'DecompositionOps_include.f90:220', &
              !              errst=errst)) return

              call create(Vt, [nkeep, Tens%dl(nl + 1:)])
              Vt%elem = reshape(utmp(first:, :), [nkeep * a_col])
      end select

      ! Permutation on the output
      ! -------------------------

      if(present(lpermout)) then
          call transposed(Ut, lpermout, doperm=.true., errst=errst)
          !if(prop_error('eigsvd_tensorc: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Vt, rpermout, doperm=.true., errst=errst)
          !if(prop_error('eigsvd_tensorc: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine eigsvd_tensorc

DecompositionOps_f90.get_chi_mps()

fortran-subroutine - May 2017 (dj) Calculate the new bond dimension based on the truncation criteria.

Arguments

LamTYPE(tensor), in

Contains the singular values or singular values squared.

callerCHARACTER, in

Called from eigendecomposition (‘E’). Singular values are squared. Called from singular values decomposition (‘S’).

chiINTEGER, out

On exit the new bond dimension according to the truncation criteria.

truncREAL, OPTIONAL, in

Specifies the allowed truncation in the norm as soft cut-off. Default to zero.

ncutINTEGER, OPTIONAL, in

Specifies the maximal number of singular values as hard cut off.

errREAL, OPTIONAL, out

On exit the truncation error.

renormREAL, OPTIONAL, out

On exit the dot-product of the singular values with themselves, which can be used for renormalizing them.

Source Code

show / hide f90 code

  subroutine get_chi_mps(Lam, caller, chi, trunc, ncut, err, renorm, errst)
      use BasicOps, only : numzero

      type(tensor), intent(in) :: Lam
      character, intent(in) :: caller
      integer, intent(out) :: chi
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      real(KIND=rKind), intent(out), optional :: renorm
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii, upper

      ! Duplettes for optional arguments
      real(KIND=rKind) :: trunc_

      ! norms of values kept and truncated
      real(KIND=rKind) :: normk, normt, tmp

      ! complete norm of all Lam
      real(KIND=rKind) :: cnorm

      ! BLAS scalar product
      real(KIND=rKind), external :: ddot

      !if(present(errst)) errst = 0

      trunc_ = numzero
      if(present(trunc)) trunc_ = trunc

      chi = Lam%dim
      if(present(ncut)) chi = min(ncut, chi)

      if(caller == 'E') then
          ! Call from eigenvalue SVD: values ascending, squared
          ! ...................................................

          normk = sum(Lam%elem(Lam%dim - chi + 1:Lam%dim))
          normt = sum(Lam%elem(1:Lam%dim - chi))
          cnorm = normk + normt

          if(normt / cnorm < trunc_) then
              ! Shrink bond dimension

              upper = Lam%dim - chi + 1
              do ii = upper, (Lam%dim -1)
                  !tmp = Lam%elem(ii)

                  if((normt + Lam%elem(ii)) / cnorm < trunc_) then
                      ! Continue shrinking
                      chi = chi - 1
                      normt = normt + Lam%elem(ii)
                      normk = normk - Lam%elem(ii)
                  else
                      ! Leave loop
                      exit
                  end if
              end do
          end if

      elseif(caller == 'S') then
          ! Call from SVD: values descending, not squared
          ! .............................................

          normk = ddot(chi, Lam%elem, 1, Lam%elem, 1)
          if(chi < Lam%dim) then
              normt = ddot(Lam%dim - chi, Lam%elem(chi + 1:Lam%dim), 1, &
                           Lam%elem(chi + 1:Lam%dim), 1)
          else
              normt = 0.0_rKind
          end if
          cnorm = normk + normt

          ! The local truncation error is
          ! 1 - sqrt(\sum_{i=1}^{\chi} \lambda_i^2) / cnorm
          ! <= \sum_{i = \chi + 1}^{chimax} \lambda_i^2 / cnorm^2
          if(normt / cnorm < trunc_) then
              ! Shrink bond dimension

              upper = chi
              do ii = upper, 2, -1
                  tmp = Lam%elem(ii)**2

                  if((normt + tmp) / cnorm < trunc_) then
                      ! Continue shrinking
                      chi = ii - 1
                      normt = normt + tmp
                      normk = normk - tmp
                  else
                      ! Leave loop
                      exit
                  end if
              end do

          end if

      else
          stop 'Bad argument.'
      end if

      if(present(renorm)) renorm = normk
      if(present(err)) err = normt / cnorm

  end subroutine get_chi_mps

DecompositionOps_f90.permute_split_tensor()

fortran-subroutine - May 2017 (dj) Perform permutation to bring indices in idxl to the left and indices in idxr to the right.

Arguments

TensTYPE(tensor), inout

Permute indices in this tensor.

idxlINTEGER(*), in

Contains indices which are permuted to the front.

idxrINTEGER(*), in

Contains indices which are permuted to the back.

Source Code

show / hide f90 code

  subroutine permute_split_tensor(Tens, idxl, idxr, errst)
      type(tensor), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! number of the left/right indices
      integer :: nl, nr

      ! array for permutation
      integer, dimension(:), allocatable :: perm

      nl = size(idxl)
      nr = size(idxr)

      allocate(perm(nl + nr))
      perm(:nl) = idxl
      perm(nl + 1:) = idxr

      call transposed(Tens, perm=perm, doperm=.true., errst=errst)
      !if(prop_error('permute_split_tensor: transpose failed', &
      !              errst=errst)) return
      deallocate(perm)

  end subroutine permute_split_tensor

DecompositionOps_f90.permute_split_tensorc()

fortran-subroutine - May 2017 (dj) Perform permutation to bring indices in idxl to the left and indices in idxr to the right.

Arguments

TensTYPE(tensorc), inout

Permute indices in this tensor.

idxlINTEGER(*), in

Contains indices which are permuted to the front.

idxrINTEGER(*), in

Contains indices which are permuted to the back.

Source Code

show / hide f90 code

  subroutine permute_split_tensorc(Tens, idxl, idxr, errst)
      type(tensorc), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! number of the left/right indices
      integer :: nl, nr

      ! array for permutation
      integer, dimension(:), allocatable :: perm

      nl = size(idxl)
      nr = size(idxr)

      allocate(perm(nl + nr))
      perm(:nl) = idxl
      perm(nl + 1:) = idxr

      call transposed(Tens, perm=perm, doperm=.true., errst=errst)
      !if(prop_error('permute_split_tensorc: transpose failed', &
      !              errst=errst)) return
      deallocate(perm)

  end subroutine permute_split_tensorc

DecompositionOps_f90.qr_tensor()

fortran-subroutine - November 2016 (dj) QR decomposition of a tensor.

Arguments

TensTYPE(tensor), inout

Decompose this tensor in a QR decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(tensor), inout

R matrix of the QR decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (Q).

idxrINTEGER(*), in

Indices of the legs going into the right tensor (R).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor Q before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor R before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine qr_tensor(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(tensor), intent(inout) :: Tens, Rtens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! dimensions of columns and rows
      integer :: a_row, a_col

      ! number of left / right legs
      integer :: nl, nr

      ! dimension of tensor for R-tensor
      integer, dimension(:), allocatable :: rdl

      ! dimension of tensor for trapezoidal case
      integer, dimension(:), allocatable :: qdl

      ! temporary tensor for matrix Q
      type(tensor) :: Qtens

      !if(present(errst)) errst = 0

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('qr_tensor: permute_split failed', &
      !              errst=errst)) return

      ! Run the QR
      ! ----------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      allocate(rdl(nr + 1))
      rdl(1) = min(a_row, a_col)
      rdl(2:) = Tens%dl(nl + 1:)
      call create(Rtens, rdl)
      deallocate(rdl)

      if(a_row >= a_col) then
          ! Regular case: Rtens is upper triangular
          call qr_inplace_real(Tens%elem, a_row, a_col, a_row, &
                                      Rtens%elem, errst=errst)
          !if(prop_error('qr_tensor: qr (1) failed', &
          !              errst=errst)) return

          if(nr > 1) then
              ! Re-use rdl for left tensor
              allocate(rdl(nl + 1))
              rdl(:nl) = Tens%dl(:nl)
              rdl(nl + 1) = a_col

              deallocate(Tens%dl)
              allocate(Tens%dl(nl + 1))
              Tens%dl = rdl

              ! Maximal dimension
              rdl(:nl) = Tens%mdl(:nl)
              rdl(nl + 1) = product(Tens%mdl(nl + 1:))

              deallocate(Tens%mdl)
              allocate(Tens%mdl(nl + 1))
              Tens%mdl = rdl

              ! Order
              deallocate(Tens%idx)
              allocate(Tens%idx(nl + 1))
              Tens%idx = [(ii, ii = 1, nl + 1)]

              ! Rank
              Tens%rank = nl + 1

              deallocate(rdl)
          end if
      else
          ! Rtens is trapezoidal
          allocate(qdl(nl + 1))
          qdl(:nl) = Tens%dl(:nl)
          qdl(nl + 1) = min(a_row, a_col)
          call create(Qtens, qdl)
          deallocate(qdl)

          call qr_trapez_real(Tens%elem, Qtens%elem, a_row, a_col, &
                                     a_row, Rtens%elem)
          !if(prop_error('qr_tensor: qr (1) failed', &
          !              errst=errst)) return

          call destroy(Tens)
          call pointto(Tens, Qtens)
          !call destroy(Qtens)
      end if

      ! Permuation on the output
      ! ------------------------

      if(present(lpermout)) then
          call transposed(Tens, lpermout, doperm=.true., errst=errst)
          !if(prop_error('qr_tensor: tranpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Rtens, rpermout, doperm=.true., errst=errst)
          !if(prop_error('qr_tensor: tranpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine qr_tensor

DecompositionOps_f90.qr_tensorc()

fortran-subroutine - November 2016 (dj) QR decomposition of a tensor.

Arguments

TensTYPE(tensorc), inout

Decompose this tensor in a QR decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(tensorc), inout

R matrix of the QR decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (Q).

idxrINTEGER(*), in

Indices of the legs going into the right tensor (R).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor Q before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor R before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine qr_tensorc(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(tensorc), intent(inout) :: Tens, Rtens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! dimensions of columns and rows
      integer :: a_row, a_col

      ! number of left / right legs
      integer :: nl, nr

      ! dimension of tensor for R-tensor
      integer, dimension(:), allocatable :: rdl

      ! dimension of tensor for trapezoidal case
      integer, dimension(:), allocatable :: qdl

      ! temporary tensor for matrix Q
      type(tensorc) :: Qtens

      !if(present(errst)) errst = 0

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('qr_tensorc: permute_split failed', &
      !              errst=errst)) return

      ! Run the QR
      ! ----------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      allocate(rdl(nr + 1))
      rdl(1) = min(a_row, a_col)
      rdl(2:) = Tens%dl(nl + 1:)
      call create(Rtens, rdl)
      deallocate(rdl)

      if(a_row >= a_col) then
          ! Regular case: Rtens is upper triangular
          call qr_inplace_complex(Tens%elem, a_row, a_col, a_row, &
                                      Rtens%elem, errst=errst)
          !if(prop_error('qr_tensorc: qr (1) failed', &
          !              errst=errst)) return

          if(nr > 1) then
              ! Re-use rdl for left tensor
              allocate(rdl(nl + 1))
              rdl(:nl) = Tens%dl(:nl)
              rdl(nl + 1) = a_col

              deallocate(Tens%dl)
              allocate(Tens%dl(nl + 1))
              Tens%dl = rdl

              ! Maximal dimension
              rdl(:nl) = Tens%mdl(:nl)
              rdl(nl + 1) = product(Tens%mdl(nl + 1:))

              deallocate(Tens%mdl)
              allocate(Tens%mdl(nl + 1))
              Tens%mdl = rdl

              ! Order
              deallocate(Tens%idx)
              allocate(Tens%idx(nl + 1))
              Tens%idx = [(ii, ii = 1, nl + 1)]

              ! Rank
              Tens%rank = nl + 1

              deallocate(rdl)
          end if
      else
          ! Rtens is trapezoidal
          allocate(qdl(nl + 1))
          qdl(:nl) = Tens%dl(:nl)
          qdl(nl + 1) = min(a_row, a_col)
          call create(Qtens, qdl)
          deallocate(qdl)

          call qr_trapez_complex(Tens%elem, Qtens%elem, a_row, a_col, &
                                     a_row, Rtens%elem)
          !if(prop_error('qr_tensorc: qr (1) failed', &
          !              errst=errst)) return

          call destroy(Tens)
          call pointto(Tens, Qtens)
          !call destroy(Qtens)
      end if

      ! Permuation on the output
      ! ------------------------

      if(present(lpermout)) then
          call transposed(Tens, lpermout, doperm=.true., errst=errst)
          !if(prop_error('qr_tensorc: tranpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Rtens, rpermout, doperm=.true., errst=errst)
          !if(prop_error('qr_tensorc: tranpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine qr_tensorc

DecompositionOps_f90.qr_qtensor()

fortran-subroutine - May 2017 (dj) QR decomposition of a tensor with symmetries.

Arguments

TensTYPE(qtensor), inout

Decompose this tensor in a QR decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(qtensor), inout

R matrix of the QR decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (Q).

idxrINTEGER(*), in

Indices of the legs going into the right tensor (R).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor Q before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor R before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine qr_qtensor(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(qtensor), intent(inout) :: Tens, RTens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! number of quantum numbers at splitting
      integer :: nunique

      ! number of left/right indices
      integer :: nl, nr

      ! Temporary tensor for Q
      type(qtensor) :: Qtmp

      !
      type(vector_int), dimension(:), allocatable :: Lmap, Rmap

      !
      type(vector_int), dimension(:), allocatable :: Rowcut, Colcut

      !
      integer, dimension(:), allocatable :: nts_row, nts_col, deghash, idxhash

      !
      type(tensorlist) :: Rmats, Qmats

      ! Build the blocks with the same quantum number
      ! ---------------------------------------------

      call block(Tens, idxl, idxr, Qmats, nunique, Lmap, Rmap, &
                 Rowcut, Colcut, nts_row, nts_col, deghash, idxhash, &
                 errst=errst)
      !if(prop_error('split_qtensor: block failed.', &
      !              errst=errst)) return

      ! Get arrays for results of QR
      call create(Rmats, nunique)

      ! Decompose all blocks via QR
      ! ---------------------------

      do ii = 1, nunique
          call qr(Qmats%Li(ii), Rmats%Li(ii), [1], [2], errst=errst)
          !if(prop_error('qr_qtensor: qr failed.', &
          !             errst=errst)) return
      end do

      ! Regain tensors from matrices
      ! ----------------------------

      nl = size(idxl)
      nr = size(idxr)

      ! -99 will not be referenced
      call block2tensor_left(Qtmp, Tens, Qmats, nl, nunique, -99, Lmap, &
                             Rowcut, nts_row, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensor: block2tensor_left failed.', &
      !             errst=errst)) return

      ! -99 will not be referenced
      call block2tensor_right(Rtens, Tens, Rmats, nl, nr, nunique, -99, Rmap, &
                              Colcut, nts_col, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensor: block2tensor_right failed.', &
      !             errst=errst)) return

      call destroy(Tens)
      call pointto(Tens, Qtmp)!!, errst=errst)
      !!if(prop_error('rq_qtensor: copy Tens <- Qtmp failed.', &
      !!             errst=errst)) return
      !call destroy(Qtmp)

      deallocate(idxhash, deghash)

      if(present(lpermout)) then
          call transposed(Tens, lpermout, errst=errst)
          !if(prop_error('rq_qtensor: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Rtens, rpermout, errst=errst)
          !if(prop_error('rq_qtensor: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine qr_qtensor

DecompositionOps_f90.qr_qtensorc()

fortran-subroutine - May 2017 (dj) QR decomposition of a tensor with symmetries.

Arguments

TensTYPE(qtensorc), inout

Decompose this tensor in a QR decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(qtensorc), inout

R matrix of the QR decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (Q).

idxrINTEGER(*), in

Indices of the legs going into the right tensor (R).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor Q before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor R before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine qr_qtensorc(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(qtensorc), intent(inout) :: Tens, RTens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! number of quantum numbers at splitting
      integer :: nunique

      ! number of left/right indices
      integer :: nl, nr

      ! Temporary tensor for Q
      type(qtensorc) :: Qtmp

      !
      type(vector_int), dimension(:), allocatable :: Lmap, Rmap

      !
      type(vector_int), dimension(:), allocatable :: Rowcut, Colcut

      !
      integer, dimension(:), allocatable :: nts_row, nts_col, deghash, idxhash

      !
      type(tensorlistc) :: Rmats, Qmats

      ! Build the blocks with the same quantum number
      ! ---------------------------------------------

      call block(Tens, idxl, idxr, Qmats, nunique, Lmap, Rmap, &
                 Rowcut, Colcut, nts_row, nts_col, deghash, idxhash, &
                 errst=errst)
      !if(prop_error('split_qtensorc: block failed.', &
      !              errst=errst)) return

      ! Get arrays for results of QR
      call create(Rmats, nunique)

      ! Decompose all blocks via QR
      ! ---------------------------

      do ii = 1, nunique
          call qr(Qmats%Li(ii), Rmats%Li(ii), [1], [2], errst=errst)
          !if(prop_error('qr_qtensorc: qr failed.', &
          !             errst=errst)) return
      end do

      ! Regain tensors from matrices
      ! ----------------------------

      nl = size(idxl)
      nr = size(idxr)

      ! -99 will not be referenced
      call block2tensor_left(Qtmp, Tens, Qmats, nl, nunique, -99, Lmap, &
                             Rowcut, nts_row, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensorc: block2tensor_left failed.', &
      !             errst=errst)) return

      ! -99 will not be referenced
      call block2tensor_right(Rtens, Tens, Rmats, nl, nr, nunique, -99, Rmap, &
                              Colcut, nts_col, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensorc: block2tensor_right failed.', &
      !             errst=errst)) return

      call destroy(Tens)
      call pointto(Tens, Qtmp)!!, errst=errst)
      !!if(prop_error('rq_qtensorc: copy Tens <- Qtmp failed.', &
      !!             errst=errst)) return
      !call destroy(Qtmp)

      deallocate(idxhash, deghash)

      if(present(lpermout)) then
          call transposed(Tens, lpermout, errst=errst)
          !if(prop_error('rq_qtensorc: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Rtens, rpermout, errst=errst)
          !if(prop_error('rq_qtensorc: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine qr_qtensorc

DecompositionOps_f90.rq_tensor()

fortran-subroutine - November 2016 (dj) RQ decomposition of a tensor.

Arguments

TensTYPE(tensor), inout

Decompose this tensor in a RQ decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(tensor), inout

R matrix of the RQ decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (R)

idxrINTEGER(*), in

Indices of the legs going into the right tensor (Q).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor R before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor Q before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine rq_tensor(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(tensor), intent(inout) :: Tens, Rtens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! dimensions of columns and rows
      integer :: a_row, a_col

      ! number of left / right legs
      integer :: nl, nr

      ! dimension of tensor for R-tensor
      integer, dimension(:), allocatable :: rdl

      ! dimension of tensor for trapezoidal case
      integer, dimension(:), allocatable :: qdl

      ! temporary tensor for matrix Q
      type(tensor) :: Qtens

      !if(present(errst)) errst = 0

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('rq_tensor: permute_split failed', &
      !              'DecompositionOps_include.f90:868', errst=errst)) return

      ! Run the RQ
      ! ----------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      allocate(rdl(nl + 1))
      rdl(:nl) = Tens%dl(:nl)
      rdl(nl + 1) = min(a_row, a_col)
      call create(Rtens, rdl)
      deallocate(rdl)

      if(a_row <= a_col) then
          ! Regular case: Rtens is upper triangular
          call rq_inplace_real(Tens%elem, a_row, a_col, a_row, &
                                      Rtens%elem, errst=errst)
          !if(prop_error('rq_tensor: rq (1) failed', &
          !              'DecompositionOps_include.f90:889', errst=errst)) return

          if(nl > 1) then
              ! Re-use rdl for unitary tensor dimensions
              allocate(rdl(nr + 1))

              rdl(1) = a_row
              rdl(2:) = Tens%dl(nl + 1:)

              deallocate(Tens%dl)
              allocate(Tens%dl(nr + 1))
              Tens%dl = rdl

              ! Maximal dimension
              rdl(1) = product(Tens%mdl(:nl))
              rdl(2:) = Tens%mdl(nl + 1:)

              deallocate(Tens%mdl)
              allocate(Tens%mdl(nr + 1))
              Tens%mdl = rdl

              ! Order
              deallocate(Tens%idx)
              allocate(Tens%idx(nr + 1))
              Tens%idx = [(ii, ii = 1, nr + 1)]

              ! Rank
              Tens%rank = nr + 1

              deallocate(rdl)
          end if
      else
          ! Rtens is trapezoidal
          allocate(qdl(nr + 1))
          qdl(1) = min(a_row, a_col)
          qdl(2:) = Tens%dl(nl + 1:)
          call create(Qtens, qdl)
          deallocate(qdl)

          call rq_trapez_real(Tens%elem, Qtens%elem, a_row, a_col, &
                                     a_row, Rtens%elem)
          !if(prop_error('rq_tensor: rq (2) failed', &
          !              'DecompositionOps_include.f90:931', errst=errst)) return

          call destroy(Tens)
          call pointto(Tens, Qtens)
          !call destroy(Qtens)
      end if

      ! Permutation on the output
      ! -------------------------

      if(present(lpermout)) then
          call transposed(Rtens, lpermout, doperm=.true., errst=errst)
          !if(prop_error('rq_tensor: tranpose (l) failed', &
          !              'DecompositionOps_include.f90:944', errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Tens, rpermout, doperm=.true., errst=errst)
          !if(prop_error('rq_tensor: tranpose (r) failed', &
          !              'DecompositionOps_include.f90:950', errst=errst)) return
      end if

  end subroutine rq_tensor

DecompositionOps_f90.rq_tensorc()

fortran-subroutine - November 2016 (dj) RQ decomposition of a tensor.

Arguments

TensTYPE(tensorc), inout

Decompose this tensor in a RQ decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(tensorc), inout

R matrix of the RQ decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (R)

idxrINTEGER(*), in

Indices of the legs going into the right tensor (Q).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor R before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor Q before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine rq_tensorc(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(tensorc), intent(inout) :: Tens, Rtens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! dimensions of columns and rows
      integer :: a_row, a_col

      ! number of left / right legs
      integer :: nl, nr

      ! dimension of tensor for R-tensor
      integer, dimension(:), allocatable :: rdl

      ! dimension of tensor for trapezoidal case
      integer, dimension(:), allocatable :: qdl

      ! temporary tensor for matrix Q
      type(tensorc) :: Qtens

      !if(present(errst)) errst = 0

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('rq_tensorc: permute_split failed', &
      !              'DecompositionOps_include.f90:868', errst=errst)) return

      ! Run the RQ
      ! ----------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      allocate(rdl(nl + 1))
      rdl(:nl) = Tens%dl(:nl)
      rdl(nl + 1) = min(a_row, a_col)
      call create(Rtens, rdl)
      deallocate(rdl)

      if(a_row <= a_col) then
          ! Regular case: Rtens is upper triangular
          call rq_inplace_complex(Tens%elem, a_row, a_col, a_row, &
                                      Rtens%elem, errst=errst)
          !if(prop_error('rq_tensorc: rq (1) failed', &
          !              'DecompositionOps_include.f90:889', errst=errst)) return

          if(nl > 1) then
              ! Re-use rdl for unitary tensor dimensions
              allocate(rdl(nr + 1))

              rdl(1) = a_row
              rdl(2:) = Tens%dl(nl + 1:)

              deallocate(Tens%dl)
              allocate(Tens%dl(nr + 1))
              Tens%dl = rdl

              ! Maximal dimension
              rdl(1) = product(Tens%mdl(:nl))
              rdl(2:) = Tens%mdl(nl + 1:)

              deallocate(Tens%mdl)
              allocate(Tens%mdl(nr + 1))
              Tens%mdl = rdl

              ! Order
              deallocate(Tens%idx)
              allocate(Tens%idx(nr + 1))
              Tens%idx = [(ii, ii = 1, nr + 1)]

              ! Rank
              Tens%rank = nr + 1

              deallocate(rdl)
          end if
      else
          ! Rtens is trapezoidal
          allocate(qdl(nr + 1))
          qdl(1) = min(a_row, a_col)
          qdl(2:) = Tens%dl(nl + 1:)
          call create(Qtens, qdl)
          deallocate(qdl)

          call rq_trapez_complex(Tens%elem, Qtens%elem, a_row, a_col, &
                                     a_row, Rtens%elem)
          !if(prop_error('rq_tensorc: rq (2) failed', &
          !              'DecompositionOps_include.f90:931', errst=errst)) return

          call destroy(Tens)
          call pointto(Tens, Qtens)
          !call destroy(Qtens)
      end if

      ! Permutation on the output
      ! -------------------------

      if(present(lpermout)) then
          call transposed(Rtens, lpermout, doperm=.true., errst=errst)
          !if(prop_error('rq_tensorc: tranpose (l) failed', &
          !              'DecompositionOps_include.f90:944', errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Tens, rpermout, doperm=.true., errst=errst)
          !if(prop_error('rq_tensorc: tranpose (r) failed', &
          !              'DecompositionOps_include.f90:950', errst=errst)) return
      end if

  end subroutine rq_tensorc

DecompositionOps_f90.rq_qtensor()

fortran-subroutine - May 2017 (dj) RQ decomposition of a tensor with symmetries.

Arguments

TensTYPE(qtensor), inout

Decompose this tensor in a RQ decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(qtensor), inout

R matrix of the RQ decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (R)

idxrINTEGER(*), in

Indices of the legs going into the right tensor (Q).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor R before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor Q before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine rq_qtensor(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(qtensor), intent(inout) :: Tens, Rtens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! number of quantum numbers at splitting
      integer :: nunique

      ! number of left/right indices
      integer :: nl, nr

      ! Temporary tensor for Q
      type(qtensor) :: Qtmp

      !
      type(vector_int), dimension(:), allocatable :: Lmap, Rmap

      !
      type(vector_int), dimension(:), allocatable :: Rowcut, Colcut

      !
      integer, dimension(:), allocatable :: nts_row, nts_col, deghash, idxhash

      !
      type(tensorlist) :: Rmats, Qmats

      ! Build the blocks with the same quantum number
      ! ---------------------------------------------

      call block(Tens, idxl, idxr, Qmats, nunique, Lmap, Rmap, &
                 Rowcut, Colcut, nts_row, nts_col, deghash, idxhash, &
                 errst=errst)
      !if(prop_error('split_qtensor: block failed.', &
      !              errst=errst)) return

      ! Get arrays for results of RQ
      call create(Rmats, nunique)

      ! Decompose all blocks via RQ
      ! ---------------------------

      do ii = 1, nunique
          call rq(Qmats%Li(ii), Rmats%Li(ii), [1], [2], errst=errst)
          !if(prop_error('rq_qtensor: rq failed.', &
          !             errst=errst)) return
      end do

      ! Regain tensors from matrices
      ! ----------------------------

      nl = size(idxl)
      nr = size(idxr)

      ! -99 will not be referenced
      call block2tensor_left(Rtens, Tens, Rmats, nl, nunique, -99, Lmap, &
                             Rowcut, nts_row, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensor: block2tensor_left failed.', &
      !             errst=errst)) return

      ! -99 will not be referenced
      call block2tensor_right(Qtmp, Tens, Qmats, nl, nr, nunique, -99, Rmap, &
                              Colcut, nts_col, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensor: block2tensor_right failed.', &
      !             errst=errst)) return

      call destroy(Tens)
      call pointto(Tens, Qtmp)!, errst=errst)
      !!if(prop_error('rq_qtensor: copy Tens <- Qtmp failed.', &
      !!             errst=errst)) return
      !call destroy(Qtmp)

      deallocate(idxhash, deghash)

      if(present(lpermout)) then
          call transposed(Rtens, lpermout, errst=errst)
          !if(prop_error('rq_qtensor: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Tens, rpermout, errst=errst)
          !if(prop_error('rq_qtensor: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine rq_qtensor

DecompositionOps_f90.rq_qtensorc()

fortran-subroutine - May 2017 (dj) RQ decomposition of a tensor with symmetries.

Arguments

TensTYPE(qtensorc), inout

Decompose this tensor in a RQ decomposition according to the indices specified. The Q matrix is stored on exit here in tensor form.

RtensTYPE(qtensorc), inout

R matrix of the RQ decomposition in tensor form.

idxlINTEGER(*), in

Indices of the legs going into the left tensor (R)

idxrINTEGER(*), in

Indices of the legs going into the right tensor (Q).

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor R before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor Q before returning it. Default to no permutation.

Source Code

show / hide f90 code

  subroutine rq_qtensorc(Tens, Rtens, idxl, idxr, lpermout, rpermout, errst)
      type(qtensorc), intent(inout) :: Tens, Rtens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! number of quantum numbers at splitting
      integer :: nunique

      ! number of left/right indices
      integer :: nl, nr

      ! Temporary tensor for Q
      type(qtensorc) :: Qtmp

      !
      type(vector_int), dimension(:), allocatable :: Lmap, Rmap

      !
      type(vector_int), dimension(:), allocatable :: Rowcut, Colcut

      !
      integer, dimension(:), allocatable :: nts_row, nts_col, deghash, idxhash

      !
      type(tensorlistc) :: Rmats, Qmats

      ! Build the blocks with the same quantum number
      ! ---------------------------------------------

      call block(Tens, idxl, idxr, Qmats, nunique, Lmap, Rmap, &
                 Rowcut, Colcut, nts_row, nts_col, deghash, idxhash, &
                 errst=errst)
      !if(prop_error('split_qtensorc: block failed.', &
      !              errst=errst)) return

      ! Get arrays for results of RQ
      call create(Rmats, nunique)

      ! Decompose all blocks via RQ
      ! ---------------------------

      do ii = 1, nunique
          call rq(Qmats%Li(ii), Rmats%Li(ii), [1], [2], errst=errst)
          !if(prop_error('rq_qtensorc: rq failed.', &
          !             errst=errst)) return
      end do

      ! Regain tensors from matrices
      ! ----------------------------

      nl = size(idxl)
      nr = size(idxr)

      ! -99 will not be referenced
      call block2tensor_left(Rtens, Tens, Rmats, nl, nunique, -99, Lmap, &
                             Rowcut, nts_row, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensorc: block2tensor_left failed.', &
      !             errst=errst)) return

      ! -99 will not be referenced
      call block2tensor_right(Qtmp, Tens, Qmats, nl, nr, nunique, -99, Rmap, &
                              Colcut, nts_col, deghash, idxhash, errst=errst)
      !if(prop_error('rq_qtensorc: block2tensor_right failed.', &
      !             errst=errst)) return

      call destroy(Tens)
      call pointto(Tens, Qtmp)!, errst=errst)
      !!if(prop_error('rq_qtensorc: copy Tens <- Qtmp failed.', &
      !!             errst=errst)) return
      !call destroy(Qtmp)

      deallocate(idxhash, deghash)

      if(present(lpermout)) then
          call transposed(Rtens, lpermout, errst=errst)
          !if(prop_error('rq_qtensorc: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Tens, rpermout, errst=errst)
          !if(prop_error('rq_qtensorc: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine rq_qtensorc

DecompositionOps_f90.split_tensor()

fortran-subroutine - November 2016 (dj) Split a tensor into two tensors using a SVD or SVD via eigenvalue decomposition.

Arguments

UtTYPE(tensor), inout

Contains on exit the left unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

LamTYPE(tensor), inout

Contains on exit the singular values. Intent in to keep dimension etc.

VtTYPE(tensor), inout

Contains on exit the right unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

TensTYPE(tensor), inout

Contains the tensor to be decomposed. Destroyed on exit.

idxlINTEGER(*), in

Indices of the legs going into the left tensor.

idxrINTEGER(*), in

Indices of the legs going into the right tensor.

multlrINTEGER, OPTIONAL, in

Multiply singular values to the left (< 0) or right (> 0). Default to no multiplication (0).

truncREAL, OPTIONAL, in

Keep infidelity below trunc (infidelity is sum of squared discarded singular values for MPS). Default does not truncate.

ncutINTEGER, OPTIONAL, in

Maximal bond dimension / number of singular values. Default is keeping all singular values.

errREAL, OPTIONAL, out

Sum over the square of the discarded singular value for MPS.

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor before returning it. Default to no permutation.

renormCHARACTER, OPTIONAL, in

Renormalization of the singular values. ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

methodCHARACTER, OPTIONAL, in

Choose method: ‘Y’ GESDD method; ‘N’ use GESVD; ‘E’ eigsvd. If not specified, …

Source Code

show / hide f90 code

  subroutine split_tensor(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, ncut, &
                               err, lpermout, rpermout, renorm, method, errst)
      type(tensor), intent(inout) :: Ut, Vt
      type(tensor), intent(inout) :: Lam
      type(tensor), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(in), optional :: multlr
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      character, intent(in), optional :: renorm, method
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! flag for SVD vs EIGSVD
      character :: switch

      if(present(method)) then
          if(method == 'Y' .or. method == 'N') then
              switch = 'S'
          else
              switch = 'E'
              !if(multlr == 0) then
              !    errst = raise_error('split_tensor:'//&
              !                        'multlr=0 + eigsvd', 99, &
              !                        errst=errst)
              !    return
              !end if
          end if
      else
          if(.not. present(multlr)) then
              switch = 'S'
          elseif(multlr == 0) then
              switch = 'S'
          else
              switch = 'E'
          end if
      end if

      select case(switch)
          case('S')
              call svd(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, &
                       ncut, err, lpermout, rpermout, renorm, method, errst)
              !if(prop_error('split_tensor: svd failed.', &
              !              errst=errst)) return
          case('E')
              call eigsvd(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, &
                          ncut, err, lpermout, rpermout, renorm, errst)
              !if(prop_error('split_tensor: eigsvd failed.', &
              !              errst=errst)) return
      end select

  end subroutine split_tensor

DecompositionOps_f90.split_tensorc()

fortran-subroutine - November 2016 (dj) Split a tensor into two tensors using a SVD or SVD via eigenvalue decomposition.

Arguments

UtTYPE(tensorc), inout

Contains on exit the left unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

LamTYPE(tensor), inout

Contains on exit the singular values. Intent in to keep dimension etc.

VtTYPE(tensorc), inout

Contains on exit the right unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

TensTYPE(tensorc), inout

Contains the tensor to be decomposed. Destroyed on exit.

idxlINTEGER(*), in

Indices of the legs going into the left tensor.

idxrINTEGER(*), in

Indices of the legs going into the right tensor.

multlrINTEGER, OPTIONAL, in

Multiply singular values to the left (< 0) or right (> 0). Default to no multiplication (0).

truncREAL, OPTIONAL, in

Keep infidelity below trunc (infidelity is sum of squared discarded singular values for MPS). Default does not truncate.

ncutINTEGER, OPTIONAL, in

Maximal bond dimension / number of singular values. Default is keeping all singular values.

errREAL, OPTIONAL, out

Sum over the square of the discarded singular value for MPS.

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor before returning it. Default to no permutation.

renormCHARACTER, OPTIONAL, in

Renormalization of the singular values. ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

methodCHARACTER, OPTIONAL, in

Choose method: ‘Y’ GESDD method; ‘N’ use GESVD; ‘E’ eigsvd. If not specified, …

Source Code

show / hide f90 code

  subroutine split_tensorc(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, ncut, &
                               err, lpermout, rpermout, renorm, method, errst)
      type(tensorc), intent(inout) :: Ut, Vt
      type(tensor), intent(inout) :: Lam
      type(tensorc), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(in), optional :: multlr
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      character, intent(in), optional :: renorm, method
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! flag for SVD vs EIGSVD
      character :: switch

      if(present(method)) then
          if(method == 'Y' .or. method == 'N') then
              switch = 'S'
          else
              switch = 'E'
              !if(multlr == 0) then
              !    errst = raise_error('split_tensorc:'//&
              !                        'multlr=0 + eigsvd', 99, &
              !                        errst=errst)
              !    return
              !end if
          end if
      else
          if(.not. present(multlr)) then
              switch = 'S'
          elseif(multlr == 0) then
              switch = 'S'
          else
              switch = 'E'
          end if
      end if

      select case(switch)
          case('S')
              call svd(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, &
                       ncut, err, lpermout, rpermout, renorm, method, errst)
              !if(prop_error('split_tensorc: svd failed.', &
              !              errst=errst)) return
          case('E')
              call eigsvd(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, &
                          ncut, err, lpermout, rpermout, renorm, errst)
              !if(prop_error('split_tensorc: eigsvd failed.', &
              !              errst=errst)) return
      end select

  end subroutine split_tensorc

DecompositionOps_f90.svd_tensor()

fortran-subroutine - November 2016 (dj) Split a tensor into two tensors using a SVD.

Arguments

UtTYPE(tensor), inout

Contains on exit the left unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

LamTYPE(tensor), inout

Contains on exit the singular values. Intent in to keep dimension etc.

VtTYPE(tensor), inout

Contains on exit the right unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

TensTYPE(tensor), inout

Contains the tensor to be decomposed. Destroyed on exit.

idxlINTEGER(*), in

Indices of the legs going into the left tensor.

idxrINTEGER(*), in

Indices of the legs going into the right tensor.

multlrINTEGER, OPTIONAL, in

Multiply singular values to the left (< 0) or right (> 0). Default to no multiplication (0).

truncREAL, OPTIONAL, in

Keep infidelity below trunc (infidelity is sum of squared discarded singular values for MPS). Default does not truncate.

ncutINTEGER, OPTIONAL, in

Maximal bond dimension / number of singular values. Default is keeping all singular values.

errREAL, OPTIONAL, out

Sum over the square of the discarded singular value for MPS.

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor before returning it. Default to no permutation.

renormCHARACTER, OPTIONAL, in

Renormalization of the singular values. ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

gesddCHARACTER, OPTIONAL, in

Choose SVD method, ‘Y’ use GESDD method (default); ‘N’ use GESVD

Source Code

show / hide f90 code

  subroutine svd_tensor(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, ncut, &
                             err, lpermout, rpermout, renorm, gesdd, errst)
      type(tensor), intent(inout) :: Ut, Vt
      type(tensor), intent(inout) :: Lam
      type(tensor), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(in), optional :: multlr
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      character, intent(in), optional :: renorm, gesdd
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! LAPACK info flag
      integer :: info

      ! number of columns and rows
      integer :: a_row, a_col

      ! number of singular values kept
      integer :: nkeep

      ! number of left/right legs 
      integer :: nl, nr

      ! new dimension for target tensors
      integer, dimension(:), allocatable :: newdl

      ! duplicate flag for gesdd
      logical :: do_gesdd

      ! real saving dot product for renormalization
      real(KIND=rKind) :: renormval

      !if(present(errst)) errst = 0

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('eigsvd_tensor: permute_split failed', &
      !              errst=errst)) return

      ! Run the svd
      ! -----------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      ! Check which method should be used
      do_gesdd = .true.
      if(present(gesdd)) do_gesdd = (gesdd == 'Y')

      allocate(newdl(nl + 1))
      newdl(:nl) = Tens%dl(:nl)
      newdl(nl + 1) = min(a_row, a_col)
      call create(Ut, newdl)
      deallocate(newdl)

      allocate(newdl(nr + 1))
      newdl(1) = min(a_row, a_col)
      newdl(2:) = Tens%dl(nl + 1:)
      call create(Vt, newdl)
      deallocate(newdl)

      call create(Lam, [min(a_row, a_col)])

      if(do_gesdd) then
          if(a_row >= a_col) then
              ! mode='O' overwrites input matrix
              Ut%elem(:Tens%dim) = Tens%elem(:Tens%dim)
              call svdd_elem_real(Tens%elem, Lam%elem, Vt%elem, &
                                         Ut%elem, a_row, a_col, 'O', info)
          else
              ! mode='O' overwrites input matrix
              Vt%elem(:Tens%dim) = Tens%elem(:Tens%dim)
              call svdd_elem_real(Ut%elem, Lam%elem, Tens%elem, &
                                         Vt%elem, a_row, a_col, 'O', info)
          end if

          !if(info > 0) then
          !    ! Run again with more stable GESVD
          !    write(elog, *) 'Exit status svdd_elem', info
          !    write(elog, *) 'svd_tensor: GESDD failed, '//&
          !                   'switch to GESVD.'
          !    do_gesdd = .false.
          !end if
      end if

      if(.not. do_gesdd) then
          call svd_elem_real(Ut%elem, Lam%elem, Vt%elem, Tens%elem, &
                                    a_row, a_col, info)

          !if(info > 0) then
          !    write(elog, *) 'Exit status svd_elem', info
          !    errst = raise_error('svd_tensor: svd_elem failed', &
          !                        6, errst=errst)
          !    return
          !end if
      end if

      ! Decide on truncation
      ! --------------------

      call get_chi_mps(Lam, 'S', nkeep, trunc, ncut, err, renormval, &
                       errst=errst)

      if(nkeep < Lam%dl(1)) then
          ! Truncation in Lam
          Lam%dl(1) = nkeep
          Lam%dim = nkeep

          ! Truncation in Ut
          Ut%dl(Ut%rank) = nkeep
          Ut%dim = product(Ut%dl)

          ! Vt have to manipulate
          call truncate_right(Vt, nkeep)
      end if

      if(present(renorm)) then
          select case(renorm)
              case('N')
                  ! Nothing
              case('R')
                  call scale(1.0_rKind / renormval, Lam)
              case default
                  call scale(1.0_rKind / sqrt(renormval), Lam)
          end select
      end if

      ! Mulitplication to the left or right
      ! -----------------------------------

      if(present(multlr)) then
          if(multlr < 0) then
              call contr_diag(Ut, Lam, Ut%rank, errst=errst)
              !if(prop_error('svd_tensor: diag_contr (l) failed', &
              !              errst=errst)) return
          elseif(multlr > 0) then
              call contr_diag(Vt, Lam, 1, errst=errst)
              !if(prop_error('svd_tensor: diag_contr (r) failed', &
              !              errst=errst)) return
          end if
      end if

      ! Permutation on the output
      ! -------------------------

      if(present(lpermout)) then
          call transposed(Ut, lpermout, doperm=.true., errst=errst)
          !if(prop_error('svd_tensor: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Vt, rpermout, doperm=.true., errst=errst)
          !if(prop_error('svd_tensor: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine svd_tensor

DecompositionOps_f90.svd_tensorc()

fortran-subroutine - November 2016 (dj) Split a tensor into two tensors using a SVD.

Arguments

UtTYPE(tensorc), inout

Contains on exit the left unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

LamTYPE(tensor), inout

Contains on exit the singular values. Intent in to keep dimension etc.

VtTYPE(tensorc), inout

Contains on exit the right unitary tensor unless a multiplication is specified. Intent in to keep dimension etc.

TensTYPE(tensorc), inout

Contains the tensor to be decomposed. Destroyed on exit.

idxlINTEGER(*), in

Indices of the legs going into the left tensor.

idxrINTEGER(*), in

Indices of the legs going into the right tensor.

multlrINTEGER, OPTIONAL, in

Multiply singular values to the left (< 0) or right (> 0). Default to no multiplication (0).

truncREAL, OPTIONAL, in

Keep infidelity below trunc (infidelity is sum of squared discarded singular values for MPS). Default does not truncate.

ncutINTEGER, OPTIONAL, in

Maximal bond dimension / number of singular values. Default is keeping all singular values.

errREAL, OPTIONAL, out

Sum over the square of the discarded singular value for MPS.

lpermoutINTEGER(*), OPTIONAL, in

Permutation of the left tensor before returning it. Default to no permutation.

rpermoutINTEGER(*), OPTIONAL, in

Permutation of the right tensor before returning it. Default to no permutation.

renormCHARACTER, OPTIONAL, in

Renormalization of the singular values. ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

gesddCHARACTER, OPTIONAL, in

Choose SVD method, ‘Y’ use GESDD method (default); ‘N’ use GESVD

Source Code

show / hide f90 code

  subroutine svd_tensorc(Ut, Lam, Vt, Tens, idxl, idxr, multlr, trunc, ncut, &
                             err, lpermout, rpermout, renorm, gesdd, errst)
      type(tensorc), intent(inout) :: Ut, Vt
      type(tensor), intent(inout) :: Lam
      type(tensorc), intent(inout) :: Tens
      integer, dimension(:), intent(in) :: idxl, idxr
      integer, intent(in), optional :: multlr
      real(KIND=rKind), intent(in), optional :: trunc
      integer, intent(in), optional :: ncut
      real(KIND=rKind), intent(out), optional :: err
      integer, dimension(:), intent(in), optional :: lpermout, rpermout
      character, intent(in), optional :: renorm, gesdd
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! LAPACK info flag
      integer :: info

      ! number of columns and rows
      integer :: a_row, a_col

      ! number of singular values kept
      integer :: nkeep

      ! number of left/right legs 
      integer :: nl, nr

      ! new dimension for target tensors
      integer, dimension(:), allocatable :: newdl

      ! duplicate flag for gesdd
      logical :: do_gesdd

      ! real saving dot product for renormalization
      real(KIND=rKind) :: renormval

      !if(present(errst)) errst = 0

      ! Permute them into the right order
      call permute_split(Tens, idxl, idxr, errst=errst)
      !if(prop_error('eigsvd_tensorc: permute_split failed', &
      !              errst=errst)) return

      ! Run the svd
      ! -----------

      nl = size(idxl)
      nr = size(idxr)
      a_row = product(Tens%dl(:nl))
      a_col = product(Tens%dl(nl + 1:))

      ! Check which method should be used
      do_gesdd = .true.
      if(present(gesdd)) do_gesdd = (gesdd == 'Y')

      allocate(newdl(nl + 1))
      newdl(:nl) = Tens%dl(:nl)
      newdl(nl + 1) = min(a_row, a_col)
      call create(Ut, newdl)
      deallocate(newdl)

      allocate(newdl(nr + 1))
      newdl(1) = min(a_row, a_col)
      newdl(2:) = Tens%dl(nl + 1:)
      call create(Vt, newdl)
      deallocate(newdl)

      call create(Lam, [min(a_row, a_col)])

      if(do_gesdd) then
          if(a_row >= a_col) then
              ! mode='O' overwrites input matrix
              Ut%elem(:Tens%dim) = Tens%elem(:Tens%dim)
              call svdd_elem_complex(Tens%elem, Lam%elem, Vt%elem, &
                                         Ut%elem, a_row, a_col, 'O', info)
          else
              ! mode='O' overwrites input matrix
              Vt%elem(:Tens%dim) = Tens%elem(:Tens%dim)
              call svdd_elem_complex(Ut%elem, Lam%elem, Tens%elem, &
                                         Vt%elem, a_row, a_col, 'O', info)
          end if

          !if(info > 0) then
          !    ! Run again with more stable GESVD
          !    write(elog, *) 'Exit status svdd_elem', info
          !    write(elog, *) 'svd_tensorc: GESDD failed, '//&
          !                   'switch to GESVD.'
          !    do_gesdd = .false.
          !end if
      end if

      if(.not. do_gesdd) then
          call svd_elem_complex(Ut%elem, Lam%elem, Vt%elem, Tens%elem, &
                                    a_row, a_col, info)

          !if(info > 0) then
          !    write(elog, *) 'Exit status svd_elem', info
          !    errst = raise_error('svd_tensorc: svd_elem failed', &
          !                        6, errst=errst)
          !    return
          !end if
      end if

      ! Decide on truncation
      ! --------------------

      call get_chi_mps(Lam, 'S', nkeep, trunc, ncut, err, renormval, &
                       errst=errst)

      if(nkeep < Lam%dl(1)) then
          ! Truncation in Lam
          Lam%dl(1) = nkeep
          Lam%dim = nkeep

          ! Truncation in Ut
          Ut%dl(Ut%rank) = nkeep
          Ut%dim = product(Ut%dl)

          ! Vt have to manipulate
          call truncate_right(Vt, nkeep)
      end if

      if(present(renorm)) then
          select case(renorm)
              case('N')
                  ! Nothing
              case('R')
                  call scale(1.0_rKind / renormval, Lam)
              case default
                  call scale(1.0_rKind / sqrt(renormval), Lam)
          end select
      end if

      ! Mulitplication to the left or right
      ! -----------------------------------

      if(present(multlr)) then
          if(multlr < 0) then
              call contr_diag(Ut, Lam, Ut%rank, errst=errst)
              !if(prop_error('svd_tensorc: diag_contr (l) failed', &
              !              errst=errst)) return
          elseif(multlr > 0) then
              call contr_diag(Vt, Lam, 1, errst=errst)
              !if(prop_error('svd_tensorc: diag_contr (r) failed', &
              !              errst=errst)) return
          end if
      end if

      ! Permutation on the output
      ! -------------------------

      if(present(lpermout)) then
          call transposed(Ut, lpermout, doperm=.true., errst=errst)
          !if(prop_error('svd_tensorc: transpose (l) failed', &
          !              errst=errst)) return
      end if

      if(present(rpermout)) then
          call transposed(Vt, rpermout, doperm=.true., errst=errst)
          !if(prop_error('svd_tensorc: transpose (r) failed', &
          !              errst=errst)) return
      end if

  end subroutine svd_tensorc