"""
Fortran module DecompositionOps: May 2017 (dj)
Containing decompositions for tensors.
**Authors**
* D. Jaschke
* M. L. Wall
**Details**
The following subroutines / functions are defined in the
module as interfaces.
+------------------------+
| Procedure |
+========================+
| qr |
+------------------------+
| rq |
+------------------------+
| split |
+------------------------+
"""
[docs]def eigsvd_tensor():
"""
fortran-subroutine - November 2016 (dj)
Split a tensor into two tensors using a SVD.
**Arguments**
Ut : TYPE(tensor), inout
Contains on exit the left unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Lam : TYPE(tensor), inout
Contains on exit the singular values. Intent in to keep dimension
etc.
Vt : TYPE(tensor), inout
Contains on exit the right unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Tens : TYPE(tensor), inout
Contains the tensor to be decomposed. Destroyed on exit.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor.
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor.
multlr : INTEGER, OPTIONAL, in
Multiply singular values to the left (< 0) or right (> 0). Default
to no multiplication (0), but no valid option for eigsvd.
trunc : REAL, OPTIONAL, in
Keep infidelity below trunc (infidelity is sum of squared discarded
singular values for MPS). Default does not truncate.
ncut : INTEGER, OPTIONAL, in
Maximal bond dimension / number of singular values. Default is
keeping all singular values.
err : REAL, OPTIONAL, out
Sum over the square of the discarded singular value for MPS.
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor before returning it. Default to
no permutation.
renorm : CHARACTER, OPTIONAL, in
Renormalization of the singular values.
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def eigsvd_tensorc():
"""
fortran-subroutine - November 2016 (dj)
Split a tensor into two tensors using a SVD.
**Arguments**
Ut : TYPE(tensorc), inout
Contains on exit the left unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Lam : TYPE(tensor), inout
Contains on exit the singular values. Intent in to keep dimension
etc.
Vt : TYPE(tensorc), inout
Contains on exit the right unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Tens : TYPE(tensorc), inout
Contains the tensor to be decomposed. Destroyed on exit.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor.
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor.
multlr : INTEGER, OPTIONAL, in
Multiply singular values to the left (< 0) or right (> 0). Default
to no multiplication (0), but no valid option for eigsvd.
trunc : REAL, OPTIONAL, in
Keep infidelity below trunc (infidelity is sum of squared discarded
singular values for MPS). Default does not truncate.
ncut : INTEGER, OPTIONAL, in
Maximal bond dimension / number of singular values. Default is
keeping all singular values.
err : REAL, OPTIONAL, out
Sum over the square of the discarded singular value for MPS.
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor before returning it. Default to
no permutation.
renorm : CHARACTER, OPTIONAL, in
Renormalization of the singular values.
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def get_chi_mps():
"""
fortran-subroutine - May 2017 (dj)
Calculate the new bond dimension based on the truncation criteria.
**Arguments**
Lam : TYPE(tensor), in
Contains the singular values or singular values squared.
caller : CHARACTER, in
Called from eigendecomposition ('E'). Singular values are squared.
Called from singular values decomposition ('S').
chi : INTEGER, out
On exit the new bond dimension according to the truncation
criteria.
trunc : REAL, OPTIONAL, in
Specifies the allowed truncation in the norm as soft cut-off.
Default to zero.
ncut : INTEGER, OPTIONAL, in
Specifies the maximal number of singular values as hard cut off.
err : REAL, OPTIONAL, out
On exit the truncation error.
renorm : REAL, OPTIONAL, out
On exit the dot-product of the singular values with
themselves, which can be used for renormalizing them.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
Tens : TYPE(tensor), inout
Permute indices in this tensor.
idxl : INTEGER(\*), in
Contains indices which are permuted to the front.
idxr : INTEGER(\*), in
Contains indices which are permuted to the back.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
Tens : TYPE(tensorc), inout
Permute indices in this tensor.
idxl : INTEGER(\*), in
Contains indices which are permuted to the front.
idxr : INTEGER(\*), in
Contains indices which are permuted to the back.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def qr_tensor():
"""
fortran-subroutine - November 2016 (dj)
QR decomposition of a tensor.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(tensor), inout
R matrix of the QR decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (Q).
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (R).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor Q before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor R before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def qr_tensorc():
"""
fortran-subroutine - November 2016 (dj)
QR decomposition of a tensor.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(tensorc), inout
R matrix of the QR decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (Q).
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (R).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor Q before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor R before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def qr_qtensor():
"""
fortran-subroutine - May 2017 (dj)
QR decomposition of a tensor with symmetries.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(qtensor), inout
R matrix of the QR decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (Q).
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (R).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor Q before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor R before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def qr_qtensorc():
"""
fortran-subroutine - May 2017 (dj)
QR decomposition of a tensor with symmetries.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(qtensorc), inout
R matrix of the QR decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (Q).
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (R).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor Q before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor R before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def rq_tensor():
"""
fortran-subroutine - November 2016 (dj)
RQ decomposition of a tensor.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(tensor), inout
R matrix of the RQ decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (R)
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (Q).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor R before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor Q before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def rq_tensorc():
"""
fortran-subroutine - November 2016 (dj)
RQ decomposition of a tensor.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(tensorc), inout
R matrix of the RQ decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (R)
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (Q).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor R before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor Q before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def rq_qtensor():
"""
fortran-subroutine - May 2017 (dj)
RQ decomposition of a tensor with symmetries.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(qtensor), inout
R matrix of the RQ decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (R)
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (Q).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor R before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor Q before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def rq_qtensorc():
"""
fortran-subroutine - May 2017 (dj)
RQ decomposition of a tensor with symmetries.
**Arguments**
Tens : TYPE(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.
Rtens : TYPE(qtensorc), inout
R matrix of the RQ decomposition in tensor form.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor (R)
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor (Q).
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor R before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor Q before returning it. Default to
no permutation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def split_tensor():
"""
fortran-subroutine - November 2016 (dj)
Split a tensor into two tensors using a SVD or SVD via eigenvalue
decomposition.
**Arguments**
Ut : TYPE(tensor), inout
Contains on exit the left unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Lam : TYPE(tensor), inout
Contains on exit the singular values. Intent in to keep dimension
etc.
Vt : TYPE(tensor), inout
Contains on exit the right unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Tens : TYPE(tensor), inout
Contains the tensor to be decomposed. Destroyed on exit.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor.
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor.
multlr : INTEGER, OPTIONAL, in
Multiply singular values to the left (< 0) or right (> 0). Default
to no multiplication (0).
trunc : REAL, OPTIONAL, in
Keep infidelity below trunc (infidelity is sum of squared discarded
singular values for MPS). Default does not truncate.
ncut : INTEGER, OPTIONAL, in
Maximal bond dimension / number of singular values. Default is
keeping all singular values.
err : REAL, OPTIONAL, out
Sum over the square of the discarded singular value for MPS.
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor before returning it. Default to
no permutation.
renorm : CHARACTER, OPTIONAL, in
Renormalization of the singular values.
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
method : CHARACTER, OPTIONAL, in
Choose method: 'Y' GESDD method; 'N' use GESVD; 'E' eigsvd. If
not specified, ...
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def split_tensorc():
"""
fortran-subroutine - November 2016 (dj)
Split a tensor into two tensors using a SVD or SVD via eigenvalue
decomposition.
**Arguments**
Ut : TYPE(tensorc), inout
Contains on exit the left unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Lam : TYPE(tensor), inout
Contains on exit the singular values. Intent in to keep dimension
etc.
Vt : TYPE(tensorc), inout
Contains on exit the right unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Tens : TYPE(tensorc), inout
Contains the tensor to be decomposed. Destroyed on exit.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor.
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor.
multlr : INTEGER, OPTIONAL, in
Multiply singular values to the left (< 0) or right (> 0). Default
to no multiplication (0).
trunc : REAL, OPTIONAL, in
Keep infidelity below trunc (infidelity is sum of squared discarded
singular values for MPS). Default does not truncate.
ncut : INTEGER, OPTIONAL, in
Maximal bond dimension / number of singular values. Default is
keeping all singular values.
err : REAL, OPTIONAL, out
Sum over the square of the discarded singular value for MPS.
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor before returning it. Default to
no permutation.
renorm : CHARACTER, OPTIONAL, in
Renormalization of the singular values.
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
method : CHARACTER, OPTIONAL, in
Choose method: 'Y' GESDD method; 'N' use GESVD; 'E' eigsvd. If
not specified, ...
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def svd_tensor():
"""
fortran-subroutine - November 2016 (dj)
Split a tensor into two tensors using a SVD.
**Arguments**
Ut : TYPE(tensor), inout
Contains on exit the left unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Lam : TYPE(tensor), inout
Contains on exit the singular values. Intent in to keep dimension
etc.
Vt : TYPE(tensor), inout
Contains on exit the right unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Tens : TYPE(tensor), inout
Contains the tensor to be decomposed. Destroyed on exit.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor.
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor.
multlr : INTEGER, OPTIONAL, in
Multiply singular values to the left (< 0) or right (> 0). Default
to no multiplication (0).
trunc : REAL, OPTIONAL, in
Keep infidelity below trunc (infidelity is sum of squared discarded
singular values for MPS). Default does not truncate.
ncut : INTEGER, OPTIONAL, in
Maximal bond dimension / number of singular values. Default is
keeping all singular values.
err : REAL, OPTIONAL, out
Sum over the square of the discarded singular value for MPS.
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor before returning it. Default to
no permutation.
renorm : CHARACTER, OPTIONAL, in
Renormalization of the singular values.
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
gesdd : CHARACTER, OPTIONAL, in
Choose SVD method, 'Y' use GESDD method (default); 'N' use GESVD
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def svd_tensorc():
"""
fortran-subroutine - November 2016 (dj)
Split a tensor into two tensors using a SVD.
**Arguments**
Ut : TYPE(tensorc), inout
Contains on exit the left unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Lam : TYPE(tensor), inout
Contains on exit the singular values. Intent in to keep dimension
etc.
Vt : TYPE(tensorc), inout
Contains on exit the right unitary tensor unless a multiplication
is specified. Intent in to keep dimension etc.
Tens : TYPE(tensorc), inout
Contains the tensor to be decomposed. Destroyed on exit.
idxl : INTEGER(\*), in
Indices of the legs going into the left tensor.
idxr : INTEGER(\*), in
Indices of the legs going into the right tensor.
multlr : INTEGER, OPTIONAL, in
Multiply singular values to the left (< 0) or right (> 0). Default
to no multiplication (0).
trunc : REAL, OPTIONAL, in
Keep infidelity below trunc (infidelity is sum of squared discarded
singular values for MPS). Default does not truncate.
ncut : INTEGER, OPTIONAL, in
Maximal bond dimension / number of singular values. Default is
keeping all singular values.
err : REAL, OPTIONAL, out
Sum over the square of the discarded singular value for MPS.
lpermout : INTEGER(\*), OPTIONAL, in
Permutation of the left tensor before returning it. Default to
no permutation.
rpermout : INTEGER(\*), OPTIONAL, in
Permutation of the right tensor before returning it. Default to
no permutation.
renorm : CHARACTER, OPTIONAL, in
Renormalization of the singular values.
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
gesdd : CHARACTER, OPTIONAL, in
Choose SVD method, 'Y' use GESDD method (default); 'N' use GESVD
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return