"""
Fortran module LinearAlgebra: May 2017 (dj)
Containing basic linear algebra operators defined for arrays.
**Authors**
* D. Jaschke
* M. L. Wall
**Details**
The following subroutines / functions are defined in the
linear algebra module as interfaces.
+-------------------------+
| Procedure |
+=========================+
| diag_times_mat |
+-------------------------+
| mat_times_diag |
+-------------------------+
| diag_contr_tensor |
+-------------------------+
| eigd |
+-------------------------+
| eigsvd |
+-------------------------+
| expm |
+-------------------------+
| expmh |
+-------------------------+
| kron |
+-------------------------+
| trace_rho_x_mat |
+-------------------------+
| tracematmul |
+-------------------------+
Various QR, RQ, and SVD are defined in addition.
"""
[docs]def diag_times_mat_real_real():
"""
fortran-subroutine - November 2016 (dj)
Multiplication of diagonal matrix with matrix.
**Arguments**
mat : real(nrow, ncol), inout
Matrix to be multiplied with diagonal matrix in vector form.
diag : real(nrow), in
Diagonal matrix in vector form.
nrow : INTEGER, in
Number of rows in matrix and length of vector for diagonal matrix.
ncol : INTEGER, in
Number of columns in matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine diag_times_mat_real_real(mat, diag, nrow, ncol, &
errst)
integer, intent(in) :: nrow, ncol
real(KIND=rKind), dimension(nrow, ncol), intent(inout) :: mat
real(KIND=rKind), dimension(nrow), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
do ii = 1, ncol
mat(:, ii) = diag(:) * mat(:, ii)
end do
end subroutine diag_times_mat_real_real
"""
return
[docs]def diag_times_mat_complex_real():
"""
fortran-subroutine - November 2016 (dj)
Multiplication of diagonal matrix with matrix.
**Arguments**
mat : complex(nrow, ncol), inout
Matrix to be multiplied with diagonal matrix in vector form.
diag : real(nrow), in
Diagonal matrix in vector form.
nrow : INTEGER, in
Number of rows in matrix and length of vector for diagonal matrix.
ncol : INTEGER, in
Number of columns in matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine diag_times_mat_complex_real(mat, diag, nrow, ncol, &
errst)
integer, intent(in) :: nrow, ncol
complex(KIND=rKind), dimension(nrow, ncol), intent(inout) :: mat
real(KIND=rKind), dimension(nrow), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
do ii = 1, ncol
mat(:, ii) = diag(:) * mat(:, ii)
end do
end subroutine diag_times_mat_complex_real
"""
return
[docs]def diag_times_mat_complex_complex():
"""
fortran-subroutine - November 2016 (dj)
Multiplication of diagonal matrix with matrix.
**Arguments**
mat : complex(nrow, ncol), inout
Matrix to be multiplied with diagonal matrix in vector form.
diag : complex(nrow), in
Diagonal matrix in vector form.
nrow : INTEGER, in
Number of rows in matrix and length of vector for diagonal matrix.
ncol : INTEGER, in
Number of columns in matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine diag_times_mat_complex_complex(mat, diag, nrow, ncol, &
errst)
integer, intent(in) :: nrow, ncol
complex(KIND=rKind), dimension(nrow, ncol), intent(inout) :: mat
complex(KIND=rKind), dimension(nrow), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
do ii = 1, ncol
mat(:, ii) = diag(:) * mat(:, ii)
end do
end subroutine diag_times_mat_complex_complex
"""
return
[docs]def mat_times_diag_real_real():
"""
fortran-subroutine - November 2016 (dj)
Multiplication of matrix with diagonal matrix.
**Arguments**
mat : real(nrow, ncol), inout
Matrix to be multiplied with diagonal matrix in vector form.
diag : real(nrow), in
Diagonal matrix in vector form.
nrow : INTEGER, in
Number of rows in matrix.
ncol : INTEGER, in
Number of columns in matrix and length of vector for diagonal matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine mat_times_diag_real_real(mat, diag, nrow, ncol, &
errst)
integer, intent(in) :: nrow, ncol
real(KIND=rKind), dimension(nrow, ncol), intent(inout) :: mat
real(KIND=rKind), dimension(ncol), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
do ii = 1, ncol
mat(:, ii) = diag(ii) * mat(:, ii)
end do
end subroutine mat_times_diag_real_real
"""
return
[docs]def mat_times_diag_complex_real():
"""
fortran-subroutine - November 2016 (dj)
Multiplication of matrix with diagonal matrix.
**Arguments**
mat : complex(nrow, ncol), inout
Matrix to be multiplied with diagonal matrix in vector form.
diag : real(nrow), in
Diagonal matrix in vector form.
nrow : INTEGER, in
Number of rows in matrix.
ncol : INTEGER, in
Number of columns in matrix and length of vector for diagonal matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine mat_times_diag_complex_real(mat, diag, nrow, ncol, &
errst)
integer, intent(in) :: nrow, ncol
complex(KIND=rKind), dimension(nrow, ncol), intent(inout) :: mat
real(KIND=rKind), dimension(ncol), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
do ii = 1, ncol
mat(:, ii) = diag(ii) * mat(:, ii)
end do
end subroutine mat_times_diag_complex_real
"""
return
[docs]def mat_times_diag_complex_complex():
"""
fortran-subroutine - November 2016 (dj)
Multiplication of matrix with diagonal matrix.
**Arguments**
mat : complex(nrow, ncol), inout
Matrix to be multiplied with diagonal matrix in vector form.
diag : complex(nrow), in
Diagonal matrix in vector form.
nrow : INTEGER, in
Number of rows in matrix.
ncol : INTEGER, in
Number of columns in matrix and length of vector for diagonal matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine mat_times_diag_complex_complex(mat, diag, nrow, ncol, &
errst)
integer, intent(in) :: nrow, ncol
complex(KIND=rKind), dimension(nrow, ncol), intent(inout) :: mat
complex(KIND=rKind), dimension(ncol), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
do ii = 1, ncol
mat(:, ii) = diag(ii) * mat(:, ii)
end do
end subroutine mat_times_diag_complex_complex
"""
return
[docs]def diag_contr_tensor_real_real():
"""
fortran-subroutine - November 2016 (dj)
Contract a rank-three tensor with a diagonal matrix. Together with the
diagonal matrix times matrix this covers all cases for contraction one
index with a diagonal matrix.
**Arguments**
tens : real(n1, n2, n3), inout
Rank three tensor to be contracted with diagonal matrix in vector
form over second index.
diag : real(n2), in
Diagonal matrix in vector form.
n1 : INTEGER, in
Dimension of the first index of the tensor.
n2 : INTEGER, in
Dimension of the second index of the tensor and the diagonal matrix.
n3 : INTEGER, in
Dimension of the third index of the tensor.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine diag_contr_tensor_real_real(tens, diag, n1, n2, n3, &
errst)
integer, intent(in) :: n1, n2, n3
real(KIND=rKind), dimension(n1, n2, n3), intent(inout) :: tens
real(KIND=rKind), dimension(n2), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
!if(present(errst)) errst = 0
do ii = 1, n2
do jj = 1, n3
tens(:, ii, jj) = diag(ii) * tens(:, ii, jj)
end do
end do
end subroutine diag_contr_tensor_real_real
"""
return
[docs]def diag_contr_tensor_complex_real():
"""
fortran-subroutine - November 2016 (dj)
Contract a rank-three tensor with a diagonal matrix. Together with the
diagonal matrix times matrix this covers all cases for contraction one
index with a diagonal matrix.
**Arguments**
tens : complex(n1, n2, n3), inout
Rank three tensor to be contracted with diagonal matrix in vector
form over second index.
diag : real(n2), in
Diagonal matrix in vector form.
n1 : INTEGER, in
Dimension of the first index of the tensor.
n2 : INTEGER, in
Dimension of the second index of the tensor and the diagonal matrix.
n3 : INTEGER, in
Dimension of the third index of the tensor.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine diag_contr_tensor_complex_real(tens, diag, n1, n2, n3, &
errst)
integer, intent(in) :: n1, n2, n3
complex(KIND=rKind), dimension(n1, n2, n3), intent(inout) :: tens
real(KIND=rKind), dimension(n2), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
!if(present(errst)) errst = 0
do ii = 1, n2
do jj = 1, n3
tens(:, ii, jj) = diag(ii) * tens(:, ii, jj)
end do
end do
end subroutine diag_contr_tensor_complex_real
"""
return
[docs]def diag_contr_tensor_complex_complex():
"""
fortran-subroutine - November 2016 (dj)
Contract a rank-three tensor with a diagonal matrix. Together with the
diagonal matrix times matrix this covers all cases for contraction one
index with a diagonal matrix.
**Arguments**
tens : complex(n1, n2, n3), inout
Rank three tensor to be contracted with diagonal matrix in vector
form over second index.
diag : complex(n2), in
Diagonal matrix in vector form.
n1 : INTEGER, in
Dimension of the first index of the tensor.
n2 : INTEGER, in
Dimension of the second index of the tensor and the diagonal matrix.
n3 : INTEGER, in
Dimension of the third index of the tensor.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine diag_contr_tensor_complex_complex(tens, diag, n1, n2, n3, &
errst)
integer, intent(in) :: n1, n2, n3
complex(KIND=rKind), dimension(n1, n2, n3), intent(inout) :: tens
complex(KIND=rKind), dimension(n2), intent(in) :: diag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
!if(present(errst)) errst = 0
do ii = 1, n2
do jj = 1, n3
tens(:, ii, jj) = diag(ii) * tens(:, ii, jj)
end do
end do
end subroutine diag_contr_tensor_complex_complex
"""
return
[docs]def eigd_symm_real():
"""
fortran-subroutine - April 2016 (dj)
Diagonalize symmetric mata using LAPACK's DSYEVD. Return the
eigenvectors in mata, eigenvalues in ev.
**Arguments**
mata : REAL(*,*), inout
Diagonalize this matrix. On exit overwritten with eigenvectors.
ev : REAL(*), out
Eigenvalues of diagonlization are stored in this vector.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_symm_real(mata, ev, dim, errst)
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(inout) :: mata
real(KIND=rKind), dimension(dim), intent(out) :: ev
integer, intent(out), optional :: errst
! Local variables
! ---------------
real(KIND=rKind) :: query(1)
integer :: wi(1), info
real(KIND=rKind), ALLOCATABLE :: work(:)
integer, allocatable :: iwork(:)
integer :: ll
!if(present(errst)) errst = 0
! workspace query
call dsyevd('V', 'U', dim, mata, dim, ev , query, -1, wi, -1, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: DSYEVD failed', &
! 6, errst=errst)
! return
!end if
ll = floor(query(1))
allocate(work(ll), iwork(wi(1)))
call dsyevd('V', 'U', dim, mata, dim, ev, work, ll, iwork, wi(1), info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: DSYEVD failed', &
! 6, errst=errst)
! return
!end if
deallocate(work, iwork)
end subroutine eigd_symm_real
"""
return
[docs]def eigd_symm_complex():
"""
fortran-subroutine - April 2016 (dj)
Diagonalize hermitian mata using LAPACK's ZHEEV. Return the
eigenvectors in mata, eigenvalues in ev.
**Arguments**
mata : COMPLEX(*,*), inout
Diagonalize this matrix. On exit overwritten with eigenvectors.
ev : REAL(*), out
Eigenvalues of diagonlization are stored in this vector.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_symm_complex(mata, ev, dim, errst)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(inout) :: mata
real(KIND=rKind), dimension(dim), intent(out) :: ev
integer, intent(out), optional :: errst
! Local variables
! ---------------
complex(KIND=rKind) :: query(1)
INTEGER :: info
COMPLEX(KIND=rKind), ALLOCATABLE :: work(:)
REAL(KIND=rKind), ALLOCATABLE :: rwork(:)
INTEGER :: ll
!if(present(errst)) errst = 0
! workspace query
allocate(rwork(3 * dim - 2))
call zheev('V', 'U', dim, mata, dim, ev, query, -1, rwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: ZHEEV failed', &
! 6, errst=errst)
! return
!end if
ll = floor(abs(query(1)))
allocate(work(ll))
call zheev('V', 'U', dim, mata, dim, ev, work, ll, rwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: ZHEEV failed', &
! 6, errst=errst)
! return
!end if
deallocate(work, rwork)
end subroutine eigd_symm_complex
"""
return
[docs]def eigd_real():
"""
fortran-subroutine - November 2014 (dj)
Diagonalization of a general real matrix
**Arguments**
mat : REAL(*,*), inout
on input containing the matrix to be diagonalized
on exit right eigenvectors of the matrix
eival : COMPLEX(*), out
containing on exit the (possibly complex) eigenvalues
dim : INTEGER, in
dimension of the eigenvalue problem.
matinv : REAL(*,*), OPTIONAL, out
containing on exit the inverse matrix of the right
eigenvectors.
**Details**
The original matrix is regained with mat * eival * matinv
(template defined in LinearAlgebra_template)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_real(mat, eival, dim, matinv, errst)
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
complex(KIND=rKind), dimension(dim), intent(out) :: eival
real(KIND=rKind), dimension(dim, dim), intent(out), optional :: matinv
integer, intent(out), optional :: errst
! for query
real(KIND=rKind) :: query(1)
! size of workspace array / info flag
integer :: lwork, info
! pivot elements in array
integer :: ipiv(size(mat, 1))
! real and imaginary part of the eigenvalues
real(KIND=rKind) :: wr(size(mat, 1)), wi(size(mat, 1))
! work space
real(KIND=rKind), dimension(:), allocatable :: work
! temporary matrix for eigenvectors
real(KIND=rKind) :: mtemp(size(mat, 1), size(mat, 2))
!if(present(errst)) errst = 0
! Eigenvalue decomposition
! ------------------------
call dgeev('N', 'V', dim, mat, dim, wr, wi, query, 1, &
mtemp, dim, query, -1, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: DGEEV failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call dgeev('N', 'V', dim, mat, dim, wr, wi, work, lwork, &
mtemp, dim, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: DGEEV failed', &
! 6, errst=errst)
! return
!end if
! Build complex eigenvalues and store eigenvectors in original matrix
eival = cmplx(wr, wi, KIND=rKind)
mat = mtemp
! Calculate the inverse if asked
! ------------------------------
if(present(matinv)) then
matinv = mat
call dgetrf(dim, dim, matinv, dim, ipiv, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: DGETRF failed', &
! 6, errst=errst)
! return
!end if
if(size(work) < dim) then
deallocate(work)
allocate(work(dim))
lwork = dim
end if
call dgetri(dim, matinv, dim, ipiv, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_symm_real: DGETRI failed', &
! 6, errst=errst)
! return
!end if
end if
deallocate(work)
end subroutine eigd_real
"""
return
[docs]def eigd_complex():
"""
fortran-subroutine - November 2014 (dj)
Diagonalization of a general complex matrix
**Arguments**
mat : COMPLEX(*,*), inout
on input containing the matrix to be diagonalized
on exit right eigenvectors of the matrix
eival : COMPLEX(*), out
containing on exit the (possibly complex) eigenvalues
dim : INTEGER, in
dimension of the eigenvalue problem.
matinv : COMPLEX(*,*), OPTIONAL, out
containing on exit the inverse matrix of the right
eigenvectors.
**Details**
The original matrix is regained with mat * eival * matinv
(template defined in LinearAlgebra_template)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_complex(mat, eival, dim, matinv, errst)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
complex(KIND=rKind), dimension(dim), intent(out) :: eival
complex(KIND=rKind), dimension(dim, dim), intent(out), optional :: matinv
integer, intent(out), optional :: errst
! for query
complex(KIND=rKind) :: query(1)
! size of workspace array / info flag
integer :: lwork, info
! pivot elements in array
integer :: ipiv(dim)
! work space
complex(KIND=rKind), dimension(:), allocatable :: work
real(KIND=rKind), dimension(2 * dim) :: rwork
! temporary matrix for eigenvectors
complex(KIND=rKind) :: mtemp(dim, dim)
!if(present(errst)) errst = 0
! Eigenvalue decomposition
! ------------------------
call zgeev('N', 'V', dim, mat, dim, eival, mtemp, dim, mtemp, dim, &
query, -1, rwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_complex: ZGEEV failed', &
! 6, 'LinearAlgebra_include.f90:545', errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call zgeev('N', 'V', dim, mat, dim, eival, mtemp, dim, mtemp, dim, &
work, lwork, rwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_complex: ZGEEV failed', &
! 6, 'LinearAlgebra_include.f90:555', errst=errst)
! return
!end if
! Store eigenvectors in original matrix
mat = mtemp
! Calculate the inverse if asked
! ------------------------------
if(present(matinv)) then
matinv = mat
call zgetrf(dim, dim, matinv, dim, ipiv, info)
!if(info /= 0) then
! errst = raise_error('eigd_complex: ZGETRF failed', &
! 6, 'LinearAlgebra_include.f90:570', errst=errst)
! return
!end if
if(size(work) < dim) then
deallocate(work)
allocate(work(dim))
lwork = dim
end if
call zgetri(dim, matinv, dim, ipiv, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('eigd_complex: ZGETRI failed', &
! 6, 'LinearAlgebra_include.f90:582', errst=errst)
! return
!end if
end if
deallocate(work)
end subroutine eigd_complex
"""
return
[docs]def eigd_tri_symm_real():
"""
fortran-subroutine - April 2016 (dj)
Diagonalize symmetric tridiagonal matrix with LAPACKs DSTEQR.
**Arguments**
mata : REAL(\*, \*), out
Contains on exit the eigenvectors of the symmetric tridiagonal
matrix.
diag : REAL(\*), inout
On entry the diagonal elements of the symmetric matrix, on exit the
eigenvalues in ascending order.
supdiag : REAL(\*), inout
Super-diagonal of the symmetric matrix. Destroyed on exit.
dim : INTEGER, in
Dimension of the square matrix
**Details**
(Template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_tri_symm_real(mata, diag, supdiag, dim, errst)
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(out) :: mata
real(KIND=rKind), dimension(dim), intent(inout) :: diag
real(KIND=rKind), dimension(dim - 1), intent(inout) :: supdiag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
! lapack info flag
integer :: info
! lapacj work array
real(KIND=rKind), dimension(2 * dim - 2) :: work
!if(present(errst)) errst = 0
! Initialize to identity matrix
mata = 0.0_rKind
do ii = 1, dim
mata(ii, ii) = done
end do
call DSTEQR('I', dim, diag, supdiag, mata, dim, work, info)
!if(info /= 0) then
! errst = raise_error('eigd_tri_symm_real: DSTEQR failed', &
! 6, 'LinearAlgebra_include.f90:664', errst=errst)
! return
!end if
end subroutine eigd_tri_symm_real
"""
return
[docs]def eigd_tri_symm_complex():
"""
fortran-subroutine - April 2016 (dj)
Diagonalize hermitian tridiagonal matrix using LAPACKS dense method
ZHEEV.
**Arguments**
mata : COMPLEX(\*, \*), out
Contains on exit the eigenvectors of the symmetric tridiagonal
matrix.
diag : REAL(\*), inout
On entry the diagonal elements of the hermitian matrix, on exit the
eigenvalues in ascending order.
supdiag : COMPLEX(\*), in
Super-diagonal of the hermitian matrix.
dim : INTEGER, in
Dimension of the square matrix
**Details**
(Template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_tri_symm_complex(mata, diag, supdiag, dim, errst)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mata
real(KIND=rKind), dimension(dim), intent(inout) :: diag
complex(KIND=rKind), dimension(dim - 1), intent(in) :: supdiag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii
!if(present(errst)) errst = 0
! Initialize matrix
mata = 0.0_rKind
do ii = 1, dim
mata(ii, ii) = diag(ii)
end do
! ZHEEV uses upper triangular
do ii = 1, (dim - 1)
mata(ii, ii+1) = supdiag(ii)
end do
call eigd(mata, diag, dim, errst=errst)
!if(prop_error('eigd_tri_symm_complex: eigd failed', &
! errst=errst)) return
end subroutine eigd_tri_symm_complex
"""
return
[docs]def eigd_symm_tri_first_real():
"""
fortran-subroutine - July 2017 (dj, updated)
Get the first eigenvalue and vector of a tridiagonal and real matrix.
**Arguments**
diags : REAL(\*), inout
The diagonal entries in the tridiagonal matrix.
odiags : REAL(\*), inout
The first off-diagonal entries in the triadiagonal matrix.
eigval : REAL, out
The first eigenvalue.
eigvec : REAL(\*), inout
The first eigenvector.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigd_symm_tri_first_real(diags, odiags, eigval, eigvec, errst)
real(KIND=rKind), dimension(:), intent(inout) :: diags, odiags
real(KIND=rKind), intent(out) :: eigval
real(KIND=rKind), dimension(:), intent(inout) :: eigvec
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Need a vector for eigenvalues
real(KIND=rKind) :: eigs(size(eigvec))
! Need a matrix for eigenvectors
real(KIND=rKind) :: zz(size(eigvec), 1)
! number of eigenvalues found
integer :: numeigs
! workspace for LAPACK
integer :: iwork(5 * size(eigvec))
real(KIND=rKind) :: work(5 * size(eigvec))
! Check on convergence of ii-th eigenvalue
integer :: ifail(SIZE(eigVec))
! Reporting errors from LAPACK
integer :: info
call dstevx('V', 'I', size(eigvec), diags, odiags, 0.0, 0.0, 1, 1, &
0.0_rKind, numeigs, eigs, zz, size(eigvec), work, iwork, &
ifail, info)
!if((info /= 0) .or. (ifail(1) /= 0)) then
! errst = raise_error('eigd_symm_first_real: DSTEVX failed.', &
! 99, errst=errst)
! return
!end if
eigval = eigs(1)
eigvec = zz(:, 1)
end subroutine eigd_symm_tri_first_real
"""
return
[docs]def eigsvd_lu_real():
"""
fortran-subroutine - Eigenvalue decomposition as singular value decomposition.
?? (dj, updated)
**Arguments**
um : real(a_row, a_row), out
Unitary matrix for decomposition.
sing : real(a_row), out
Singular values (squared?)
rm : real(a_row, a_col), out
Right matrix in decomposition, i.e., singular values times
right unitary matrix.
matin : real(a_row, a_col), in
Input matrix to be decomposed.
a_row : INTEGER, in
Number of row in the original matrix.
a_col : INTEGER, in
Number of columns in the original matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_lu_real(um, sing, rm, matin, a_row, a_col, errst)
integer, intent(in) :: a_row, a_col
real(KIND=rKIND), dimension(a_row, a_row), intent(out) :: um
real(KIND=rKind), dimension(a_row), intent(out) :: sing
real(KIND=rKind), dimension(a_row, a_col), intent(out) :: rm
real(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
integer, intent(out), optional :: errst
call eigsvd_lu_step1_real(um, sing, matin, a_row, a_col, errst)
!if(prop_error('eigsvd_lu_real: step 1 failed', &
! errst=errst)) return
call eigsvd_lu_step2_real(um, rm, matin, a_row, a_col, a_row, &
errst=errst)
!if(prop_error('eigsvd_lu_real: step 2 failed', &
! errst=errst)) return
end subroutine eigsvd_lu_real
"""
return
[docs]def eigsvd_lu_complex():
"""
fortran-subroutine - Eigenvalue decomposition as singular value decomposition.
?? (dj, updated)
**Arguments**
um : complex(a_row, a_row), out
Unitary matrix for decomposition.
sing : real(a_row), out
Singular values (squared?)
rm : complex(a_row, a_col), out
Right matrix in decomposition, i.e., singular values times
right unitary matrix.
matin : complex(a_row, a_col), in
Input matrix to be decomposed.
a_row : INTEGER, in
Number of row in the original matrix.
a_col : INTEGER, in
Number of columns in the original matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_lu_complex(um, sing, rm, matin, a_row, a_col, errst)
integer, intent(in) :: a_row, a_col
complex(KIND=rKIND), dimension(a_row, a_row), intent(out) :: um
real(KIND=rKind), dimension(a_row), intent(out) :: sing
complex(KIND=rKind), dimension(a_row, a_col), intent(out) :: rm
complex(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
integer, intent(out), optional :: errst
call eigsvd_lu_step1_complex(um, sing, matin, a_row, a_col, errst)
!if(prop_error('eigsvd_lu_complex: step 1 failed', &
! errst=errst)) return
call eigsvd_lu_step2_complex(um, rm, matin, a_row, a_col, a_row, &
errst=errst)
!if(prop_error('eigsvd_lu_complex: step 2 failed', &
! errst=errst)) return
end subroutine eigsvd_lu_complex
"""
return
[docs]def eigsvd_lu_step1_real():
"""
fortran-subroutine - ?? (dj, updated)
First step of eigenvalue decomposition as singular value decomposition.
This subroutine has the unitary matrix on the left.
**Arguments**
um : real(a_row, a_row), out
Unitary matrix on the left of the decomposition on exit.
sing : real(a_row), out
Singular values squared on exit.
matin : real(a_row, a_col), in
Input matrix to be decomposed.
a_row : INTEGER, in
Number of rows in matrix to be decomposed
a_col : INTEGER, in
Number of columns in matrix to be decomposed.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_lu_step1_real(um, sing, matin, a_row, a_col, errst)
integer, intent(in) :: a_row, a_col
real(KIND=rKind), dimension(a_row, a_row), intent(out) :: um
real(KIND=rKind), dimension(a_row), intent(out) :: sing
real(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
integer, intent(out), optional :: errst
call dsyrk('U', 'N', a_row, a_col, done, matin, a_row, dzero, &
um, a_row)
call eigd(um, sing, a_row, errst=errst)
!if(prop_error('eigsvd_lu_step1_real: eigd failed', &
! errst=errst)) return
end subroutine eigsvd_lu_step1_real
"""
return
[docs]def eigsvd_lu_step1_complex():
"""
fortran-subroutine - ?? (dj, updated)
First step of eigenvalue decomposition as singular value decomposition.
This subroutine has the unitary matrix on the left.
**Arguments**
um : complex(a_row, a_row), out
Unitary matrix on the left of the decomposition on exit.
sing : real(a_row), out
Singular values squared on exit.
matin : complex(a_row, a_col), in
Input matrix to be decomposed.
a_row : INTEGER, in
Number of rows in matrix to be decomposed
a_col : INTEGER, in
Number of columns in matrix to be decomposed.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_lu_step1_complex(um, sing, matin, a_row, a_col, errst)
integer, intent(in) :: a_row, a_col
complex(KIND=rKind), dimension(a_row, a_row), intent(out) :: um
real(KIND=rKind), dimension(a_row), intent(out) :: sing
complex(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
integer, intent(out), optional :: errst
call zherk('U', 'N', a_row, a_col, zone, matin, a_row, zzero, &
um, a_row)
call eigd(um, sing, a_row, errst=errst)
!if(prop_error('eigsvd_lu_step1_complex: eigd failed', &
! errst=errst)) return
end subroutine eigsvd_lu_step1_complex
"""
return
[docs]def eigsvd_lu_step2_real():
"""
fortran-subroutine - ?? (dj, updated)
Second step of eigenvalue decomposition as singular value decomposition.
This subroutine has the unitary matrix on the left.
**Arguments**
um : real(a_row, a_row), inout
Unitary matrix, modified in case of truncation.
rt : real(nkeep, a_col), out
On exit, singular values times right unitary matrix.
matin : real(a_row, a_col), in
Original matrix to be decomposed.
a_row : INTEGER, in
Number of rows in original matrix.
a_col : INTEGER, in
Number of columns in original matrix.
nkeep : INTEGER, in
Number of singular values to be kept.
renorm : REAL, OPTIONAL, in
Renormalization for decomposition (or scaling).
Default to 1.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_lu_step2_real(um, rt, matin, a_row, &
a_col, nkeep, renorm, errst)
integer, intent(in) :: a_row, a_col, nkeep
real(KIND=rKind), dimension(a_row, a_row), intent(inout) :: um
real(KIND=rKind), dimension(nkeep, a_col), intent(out) :: rt
real(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
real(KIND=rKind), intent(in), optional :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! scalar for renormalizing
real(KIND=rKind) :: sc
!if(present(errst)) errst = 0
if(nkeep /= a_row) then
! Singular values were discarded
um(:, :nkeep) = um(:, a_row - nkeep + 1:)
!else
! ! No singular values are discarded
end if
! Take care of renormalization if necessary
sc = done
if(present(renorm)) sc = sc * renorm
call dgemm('T', 'N', nkeep, a_col, a_row, sc, um, a_row, &
matin, a_row, dzero, rt, nkeep)
end subroutine eigsvd_lu_step2_real
"""
return
[docs]def eigsvd_lu_step2_complex():
"""
fortran-subroutine - ?? (dj, updated)
Second step of eigenvalue decomposition as singular value decomposition.
This subroutine has the unitary matrix on the left.
**Arguments**
um : complex(a_row, a_row), inout
Unitary matrix, modified in case of truncation.
rt : complex(nkeep, a_col), out
On exit, singular values times right unitary matrix.
matin : complex(a_row, a_col), in
Original matrix to be decomposed.
a_row : INTEGER, in
Number of rows in original matrix.
a_col : INTEGER, in
Number of columns in original matrix.
nkeep : INTEGER, in
Number of singular values to be kept.
renorm : REAL, OPTIONAL, in
Renormalization for decomposition (or scaling).
Default to 1.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_lu_step2_complex(um, rt, matin, a_row, &
a_col, nkeep, renorm, errst)
integer, intent(in) :: a_row, a_col, nkeep
complex(KIND=rKind), dimension(a_row, a_row), intent(inout) :: um
complex(KIND=rKind), dimension(nkeep, a_col), intent(out) :: rt
complex(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
real(KIND=rKind), intent(in), optional :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! scalar for renormalizing
complex(KIND=rKind) :: sc
!if(present(errst)) errst = 0
if(nkeep /= a_row) then
! Singular values were discarded
um(:, :nkeep) = um(:, a_row - nkeep + 1:)
!else
! ! No singular values are discarded
end if
! Take care of renormalization if necessary
sc = zone
if(present(renorm)) sc = sc * renorm
call zgemm('C', 'N', nkeep, a_col, a_row, sc, um, a_row, &
matin, a_row, zzero, rt, nkeep)
end subroutine eigsvd_lu_step2_complex
"""
return
[docs]def eigsvd_ru_step2_real():
"""
fortran-subroutine - Second step of eigenvalue decomposition as singular value decomposition.
The unitary tensor is on the right.
**Arguments**
um : real(\*, \*), in
Contains on exit the eigenvectors.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_ru_step2_real(um, lt, matin, a_row, a_col, &
nkeep, renorm, errst)
integer, intent(in) :: a_row, a_col, nkeep
real(KIND=rKind), dimension(a_col, a_col), intent(inout) :: um
real(KIND=rKind), dimension(a_row, nkeep), intent(out) :: lt
real(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
real(KIND=rKind), intent(in), optional :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
real(KIND=rKind) :: sc
!if(present(errst)) errst = 0
if(nkeep /= a_col) then
! Singular values were discarded
um(:nkeep, :) = transpose((um(:, a_col - nkeep + 1:)))
else
! All singular values, but need to transpose
um = transpose((um))
end if
! Take care of renormalization if necessary
sc = done
if(present(renorm)) sc = sc * renorm
call dgemm('N', 'T', a_row, nkeep, a_col, sc, matin, a_row, &
um, a_col, dzero, lt, a_row)
end subroutine eigsvd_ru_step2_real
"""
return
[docs]def eigsvd_ru_step2_complex():
"""
fortran-subroutine - Second step of eigenvalue decomposition as singular value decomposition.
The unitary tensor is on the right.
**Arguments**
um : complex(\*, \*), in
Contains on exit the eigenvectors.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine eigsvd_ru_step2_complex(um, lt, matin, a_row, a_col, &
nkeep, renorm, errst)
integer, intent(in) :: a_row, a_col, nkeep
complex(KIND=rKind), dimension(a_col, a_col), intent(inout) :: um
complex(KIND=rKind), dimension(a_row, nkeep), intent(out) :: lt
complex(KIND=rKind), dimension(a_row, a_col), intent(in) :: matin
real(KIND=rKind), intent(in), optional :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
complex(KIND=rKind) :: sc
!if(present(errst)) errst = 0
if(nkeep /= a_col) then
! Singular values were discarded
um(:nkeep, :) = transpose(conjg(um(:, a_col - nkeep + 1:)))
else
! All singular values, but need to transpose
um = transpose(conjg(um))
end if
! Take care of renormalization if necessary
sc = zone
if(present(renorm)) sc = sc * renorm
call zgemm('N', 'C', a_row, nkeep, a_col, sc, matin, a_row, &
um, a_col, zzero, lt, a_row)
end subroutine eigsvd_ru_step2_complex
"""
return
[docs]def expm_real_real():
"""
fortran-subroutine - April 2016 (dj)
Take the exponential of a general matrix returning a complex matrix.
**Arguments**
mexp : complex(dim, dim), in
Exponential of the matrix sc * mat
sc : real, in
Additional scalar inside exp-function.
mat : real(dim, dim), inout
On entry containing the matrix to be exponentiated. On exit
destroyed (and filled with right eigenvectors).
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expm_real_real(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
real(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! eigenvalues
complex(KIND=rKind), dimension(:), allocatable :: eival
! inverse matrix
real(KIND=rKind), dimension(:, :), allocatable :: matinv
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(matinv(dim, dim), eival(dim))
call eigd(mat, eival, dim, matinv, errst=errst)
!if(prop_error('exp_real_real: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
eival = exp(sc * eival)
!if(ieee_is_nan(abs(sum(eival)))) then
! errst = raise_error('expm_real_real : nan.', &
! 99, 'LinearAlgebra_include.f90:1253', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(eival)))) then
! errst = raise_error('expm_real_real : inf.', &
! 99, 'LinearAlgebra_include.f90:1259', errst=errst)
! return
!end if
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(ii, :) * eival(:) * matinv(:, jj))
end do
end do
deallocate(matinv, eival)
end subroutine expm_real_real
"""
return
[docs]def expm_complex_real():
"""
fortran-subroutine - April 2016 (dj)
Take the exponential of a general matrix returning a complex matrix.
**Arguments**
mexp : complex(dim, dim), in
Exponential of the matrix sc * mat
sc : complex, in
Additional scalar inside exp-function.
mat : real(dim, dim), inout
On entry containing the matrix to be exponentiated. On exit
destroyed (and filled with right eigenvectors).
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expm_complex_real(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
complex(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! eigenvalues
complex(KIND=rKind), dimension(:), allocatable :: eival
! inverse matrix
real(KIND=rKind), dimension(:, :), allocatable :: matinv
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(matinv(dim, dim), eival(dim))
call eigd(mat, eival, dim, matinv, errst=errst)
!if(prop_error('exp_complex_real: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
eival = exp(sc * eival)
!if(ieee_is_nan(abs(sum(eival)))) then
! errst = raise_error('expm_complex_real : nan.', &
! 99, 'LinearAlgebra_include.f90:1253', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(eival)))) then
! errst = raise_error('expm_complex_real : inf.', &
! 99, 'LinearAlgebra_include.f90:1259', errst=errst)
! return
!end if
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(ii, :) * eival(:) * matinv(:, jj))
end do
end do
deallocate(matinv, eival)
end subroutine expm_complex_real
"""
return
[docs]def expm_real_complex():
"""
fortran-subroutine - April 2016 (dj)
Take the exponential of a general matrix returning a complex matrix.
**Arguments**
mexp : complex(dim, dim), in
Exponential of the matrix sc * mat
sc : real, in
Additional scalar inside exp-function.
mat : complex(dim, dim), inout
On entry containing the matrix to be exponentiated. On exit
destroyed (and filled with right eigenvectors).
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expm_real_complex(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
real(KIND=rKind), intent(in) :: sc
complex(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! eigenvalues
complex(KIND=rKind), dimension(:), allocatable :: eival
! inverse matrix
complex(KIND=rKind), dimension(:, :), allocatable :: matinv
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(matinv(dim, dim), eival(dim))
call eigd(mat, eival, dim, matinv, errst=errst)
!if(prop_error('exp_real_complex: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
eival = exp(sc * eival)
!if(ieee_is_nan(abs(sum(eival)))) then
! errst = raise_error('expm_real_complex : nan.', &
! 99, 'LinearAlgebra_include.f90:1253', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(eival)))) then
! errst = raise_error('expm_real_complex : inf.', &
! 99, 'LinearAlgebra_include.f90:1259', errst=errst)
! return
!end if
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(ii, :) * eival(:) * matinv(:, jj))
end do
end do
deallocate(matinv, eival)
end subroutine expm_real_complex
"""
return
[docs]def expm_complex_complex():
"""
fortran-subroutine - April 2016 (dj)
Take the exponential of a general matrix returning a complex matrix.
**Arguments**
mexp : complex(dim, dim), in
Exponential of the matrix sc * mat
sc : complex, in
Additional scalar inside exp-function.
mat : complex(dim, dim), inout
On entry containing the matrix to be exponentiated. On exit
destroyed (and filled with right eigenvectors).
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expm_complex_complex(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
complex(KIND=rKind), intent(in) :: sc
complex(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! eigenvalues
complex(KIND=rKind), dimension(:), allocatable :: eival
! inverse matrix
complex(KIND=rKind), dimension(:, :), allocatable :: matinv
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(matinv(dim, dim), eival(dim))
call eigd(mat, eival, dim, matinv, errst=errst)
!if(prop_error('exp_complex_complex: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
eival = exp(sc * eival)
!if(ieee_is_nan(abs(sum(eival)))) then
! errst = raise_error('expm_complex_complex : nan.', &
! 99, 'LinearAlgebra_include.f90:1253', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(eival)))) then
! errst = raise_error('expm_complex_complex : inf.', &
! 99, 'LinearAlgebra_include.f90:1259', errst=errst)
! return
!end if
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(ii, :) * eival(:) * matinv(:, jj))
end do
end do
deallocate(matinv, eival)
end subroutine expm_complex_complex
"""
return
[docs]def expmh_real_real():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc * math) for hermitian math.
**Arguments**
mexp : real(dim, dim), out
Contains on exit exp(- eye * sc * math)
sc : real, in
Weight in the exponential.
mat : real(\*, \*), inout
Take exponential of this hermitian matrix. On exit, overwritten
with eigenvectors.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_real_real(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
real(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ===============
! for looping
integer :: ii, jj
! eigenvalues
real(KIND=rKind), dimension(:), allocatable :: ev
! temporary exp(eigenvalue)
real(KIND=rKind), dimension(:), allocatable :: ctmp
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(ev(dim), ctmp(dim))
call eigd(mat, ev, dim, errst=errst)
!if(prop_error('exph_real_real: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
ctmp = exp(sc * ev)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_real_real'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1361', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_real_real'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1367', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
mat = transpose(mat)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(:, ii) * ctmp(:) &
* (mat(:, jj)))
end do
end do
!if(ieee_is_nan(abs(sum(mexp)))) then
! errst = raise_error('exph_real_real : nan.', &
! 99, 'LinearAlgebra_include.f90:1383', errst=errst)
! return
!end if
deallocate(ev, ctmp)
end subroutine expmh_real_real
"""
return
[docs]def expmh_complex_real():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc * math) for hermitian math.
**Arguments**
mexp : complex(dim, dim), out
Contains on exit exp(- eye * sc * math)
sc : complex, in
Weight in the exponential.
mat : real(\*, \*), inout
Take exponential of this hermitian matrix. On exit, overwritten
with eigenvectors.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_complex_real(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
complex(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ===============
! for looping
integer :: ii, jj
! eigenvalues
real(KIND=rKind), dimension(:), allocatable :: ev
! temporary exp(eigenvalue)
complex(KIND=rKind), dimension(:), allocatable :: ctmp
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(ev(dim), ctmp(dim))
call eigd(mat, ev, dim, errst=errst)
!if(prop_error('exph_complex_real: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
ctmp = exp(sc * ev)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_complex_real'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1361', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_complex_real'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1367', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
mat = transpose(mat)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(:, ii) * ctmp(:) &
* (mat(:, jj)))
end do
end do
!if(ieee_is_nan(abs(sum(mexp)))) then
! errst = raise_error('exph_complex_real : nan.', &
! 99, 'LinearAlgebra_include.f90:1383', errst=errst)
! return
!end if
deallocate(ev, ctmp)
end subroutine expmh_complex_real
"""
return
[docs]def expmh_real_complex():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc * math) for hermitian math.
**Arguments**
mexp : complex(dim, dim), out
Contains on exit exp(- eye * sc * math)
sc : real, in
Weight in the exponential.
mat : complex(\*, \*), inout
Take exponential of this hermitian matrix. On exit, overwritten
with eigenvectors.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_real_complex(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
real(KIND=rKind), intent(in) :: sc
complex(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ===============
! for looping
integer :: ii, jj
! eigenvalues
real(KIND=rKind), dimension(:), allocatable :: ev
! temporary exp(eigenvalue)
real(KIND=rKind), dimension(:), allocatable :: ctmp
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(ev(dim), ctmp(dim))
call eigd(mat, ev, dim, errst=errst)
!if(prop_error('exph_real_complex: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
ctmp = exp(sc * ev)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_real_complex'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1361', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_real_complex'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1367', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
mat = transpose(mat)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(:, ii) * ctmp(:) &
* conjg(mat(:, jj)))
end do
end do
!if(ieee_is_nan(abs(sum(mexp)))) then
! errst = raise_error('exph_real_complex : nan.', &
! 99, 'LinearAlgebra_include.f90:1383', errst=errst)
! return
!end if
deallocate(ev, ctmp)
end subroutine expmh_real_complex
"""
return
[docs]def expmh_complex_complex():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc * math) for hermitian math.
**Arguments**
mexp : complex(dim, dim), out
Contains on exit exp(- eye * sc * math)
sc : complex, in
Weight in the exponential.
mat : complex(\*, \*), inout
Take exponential of this hermitian matrix. On exit, overwritten
with eigenvectors.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_complex_complex(mexp, sc, mat, dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
complex(KIND=rKind), intent(in) :: sc
complex(KIND=rKind), dimension(dim, dim), intent(inout) :: mat
integer, intent(out), optional :: errst
! Local variables
! ===============
! for looping
integer :: ii, jj
! eigenvalues
real(KIND=rKind), dimension(:), allocatable :: ev
! temporary exp(eigenvalue)
complex(KIND=rKind), dimension(:), allocatable :: ctmp
! Quick return for zero matrix
if(sum(abs(mat)) < 1e-15) then
mexp = 0.0_rKind
do jj = 1, dim
mexp(jj, jj) = 1.0_rKind
end do
return
end if
allocate(ev(dim), ctmp(dim))
call eigd(mat, ev, dim, errst=errst)
!if(prop_error('exph_complex_complex: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues
ctmp = exp(sc * ev)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_complex_complex'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1361', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_complex_complex'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1367', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
mat = transpose(mat)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(mat(:, ii) * ctmp(:) &
* conjg(mat(:, jj)))
end do
end do
!if(ieee_is_nan(abs(sum(mexp)))) then
! errst = raise_error('exph_complex_complex : nan.', &
! 99, 'LinearAlgebra_include.f90:1383', errst=errst)
! return
!end if
deallocate(ev, ctmp)
end subroutine expmh_complex_complex
"""
return
[docs]def expmh_tri_real_real():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.
**Arguments**
mexp : real(dim, dim), out
Contains on exit exp(- eye * sc * mat)
sc : real, in
Weight in the exponential.
diag : REAL(\*), in
Diagonal entries of the matrix to be exponentiated.
supdiag : real(\*), in
Super-diagonal entries of the matrix to be exponentiated.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_tri_real_real(mexp, sc, diag, supdiag, &
dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
real(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim), intent(in) ::diag
real(KIND=rKind), dimension(dim - 1), intent(in) :: supdiag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! temporary arrays to save intent(in)
real(KIND=rKind), dimension(:), allocatable :: cdiag
real(KIND=rKind), dimension(:), allocatable :: csupdiag
! temporary exp(eigenvalue)
real(KIND=rKind), dimension(:), allocatable :: ctmp
! eigenvectors
real(KIND=rKind), dimension(:, :), allocatable :: evec
allocate(cdiag(dim), csupdiag(dim), evec(dim, dim), ctmp(dim))
cdiag = diag
csupdiag = supdiag
call eigd(evec, cdiag, csupdiag, dim, errst=errst)
!if(prop_error('exp_tri_real_real: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues (eigenvalues must be real for hermitian
! matrix)
ctmp = exp(sc * cdiag)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_real_real'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1483', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_real_real'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1489', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
evec = transpose(evec)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(evec(:, ii) * ctmp(:) &
* (evec(:, jj)))
end do
end do
deallocate(cdiag, csupdiag, evec, ctmp)
end subroutine expmh_tri_real_real
"""
return
[docs]def expmh_tri_complex_real():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.
**Arguments**
mexp : complex(dim, dim), out
Contains on exit exp(- eye * sc * mat)
sc : complex, in
Weight in the exponential.
diag : REAL(\*), in
Diagonal entries of the matrix to be exponentiated.
supdiag : real(\*), in
Super-diagonal entries of the matrix to be exponentiated.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_tri_complex_real(mexp, sc, diag, supdiag, &
dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
complex(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim), intent(in) ::diag
real(KIND=rKind), dimension(dim - 1), intent(in) :: supdiag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! temporary arrays to save intent(in)
real(KIND=rKind), dimension(:), allocatable :: cdiag
real(KIND=rKind), dimension(:), allocatable :: csupdiag
! temporary exp(eigenvalue)
complex(KIND=rKind), dimension(:), allocatable :: ctmp
! eigenvectors
real(KIND=rKind), dimension(:, :), allocatable :: evec
allocate(cdiag(dim), csupdiag(dim), evec(dim, dim), ctmp(dim))
cdiag = diag
csupdiag = supdiag
call eigd(evec, cdiag, csupdiag, dim, errst=errst)
!if(prop_error('exp_tri_complex_real: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues (eigenvalues must be real for hermitian
! matrix)
ctmp = exp(sc * cdiag)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_complex_real'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1483', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_complex_real'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1489', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
evec = transpose(evec)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(evec(:, ii) * ctmp(:) &
* (evec(:, jj)))
end do
end do
deallocate(cdiag, csupdiag, evec, ctmp)
end subroutine expmh_tri_complex_real
"""
return
[docs]def expmh_tri_real_complex():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.
**Arguments**
mexp : complex(dim, dim), out
Contains on exit exp(- eye * sc * mat)
sc : real, in
Weight in the exponential.
diag : REAL(\*), in
Diagonal entries of the matrix to be exponentiated.
supdiag : complex(\*), in
Super-diagonal entries of the matrix to be exponentiated.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_tri_real_complex(mexp, sc, diag, supdiag, &
dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
real(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim), intent(in) ::diag
complex(KIND=rKind), dimension(dim - 1), intent(in) :: supdiag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! temporary arrays to save intent(in)
real(KIND=rKind), dimension(:), allocatable :: cdiag
complex(KIND=rKind), dimension(:), allocatable :: csupdiag
! temporary exp(eigenvalue)
real(KIND=rKind), dimension(:), allocatable :: ctmp
! eigenvectors
complex(KIND=rKind), dimension(:, :), allocatable :: evec
allocate(cdiag(dim), csupdiag(dim), evec(dim, dim), ctmp(dim))
cdiag = diag
csupdiag = supdiag
call eigd(evec, cdiag, csupdiag, dim, errst=errst)
!if(prop_error('exp_tri_real_complex: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues (eigenvalues must be real for hermitian
! matrix)
ctmp = exp(sc * cdiag)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_real_complex'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1483', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_real_complex'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1489', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
evec = transpose(evec)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(evec(:, ii) * ctmp(:) &
* conjg(evec(:, jj)))
end do
end do
deallocate(cdiag, csupdiag, evec, ctmp)
end subroutine expmh_tri_real_complex
"""
return
[docs]def expmh_tri_complex_complex():
"""
fortran-subroutine - April 2016 (dj, updated)
Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.
**Arguments**
mexp : complex(dim, dim), out
Contains on exit exp(- eye * sc * mat)
sc : complex, in
Weight in the exponential.
diag : REAL(\*), in
Diagonal entries of the matrix to be exponentiated.
supdiag : complex(\*), in
Super-diagonal entries of the matrix to be exponentiated.
dim : INTEGER, in
Dimension of the square matrix.
**Details**
(Template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine expmh_tri_complex_complex(mexp, sc, diag, supdiag, &
dim, errst)
!use IEEE_ARITHMETIC, only : IEEE_IS_FINITE, IEEE_IS_NAN
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: mexp
complex(KIND=rKind), intent(in) :: sc
real(KIND=rKind), dimension(dim), intent(in) ::diag
complex(KIND=rKind), dimension(dim - 1), intent(in) :: supdiag
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj
! temporary arrays to save intent(in)
real(KIND=rKind), dimension(:), allocatable :: cdiag
complex(KIND=rKind), dimension(:), allocatable :: csupdiag
! temporary exp(eigenvalue)
complex(KIND=rKind), dimension(:), allocatable :: ctmp
! eigenvectors
complex(KIND=rKind), dimension(:, :), allocatable :: evec
allocate(cdiag(dim), csupdiag(dim), evec(dim, dim), ctmp(dim))
cdiag = diag
csupdiag = supdiag
call eigd(evec, cdiag, csupdiag, dim, errst=errst)
!if(prop_error('exp_tri_complex_complex: eigd failed', &
! errst=errst)) return
! Exponentiate eigenvalues (eigenvalues must be real for hermitian
! matrix)
ctmp = exp(sc * cdiag)
!if(ieee_is_nan(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_complex_complex'//&
! ' : nan.', 99, 'LinearAlgebra_include.f90:1483', errst=errst)
! return
!end if
!if(.not. ieee_is_finite(abs(sum(ctmp)))) then
! errst = raise_error('expmh_tri_complex_complex'//&
! ' : inf.', 99, 'LinearAlgebra_include.f90:1489', errst=errst)
! return
!end if
! Transpose for hopefully faster matrix-matrix multiplication
evec = transpose(evec)
do jj = 1, dim
do ii = 1, dim
mexp(ii, jj) = sum(evec(:, ii) * ctmp(:) &
* conjg(evec(:, jj)))
end do
end do
deallocate(cdiag, csupdiag, evec, ctmp)
end subroutine expmh_tri_complex_complex
"""
return
[docs]def kron_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R`.
**Arguments**
Mat : real(\*, \*), inout
On exit the tensor product between the two matrices. Has to be
allocated on entry unless optional argument is set otherwise.
Matl : _real(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
Matr : real(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
transl : CHARACTER, in
Defines the transformation executed on the left matrix :math:`M_L`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
transr : CHARACTER, in
Defines the transformation executed on the left matrix :math:`M_R`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_real_real(Mat, Matl, Matr, r1, c1, r2, c2, &
transl, transr, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(r1*r2*c1*c2), intent(inout) :: Mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: Matl
real(KIND=rKind), dimension(r2, c2), intent(in) :: Matr
character, intent(in) :: transl, transr
character, intent(in), optional :: op
! Local variables
! ---------------
! Key for calling correct subroutine: 0 : M x M, 1 : M^t x M,
! 10 : M x M^t, 11 M^t x M^t, 2 : M^d x M, 20 : M x M^d
! 22 : M^d x M^d, 21 : M^d x M^t, 12 : M^t x M^d
integer :: key
! dublette for optional argument
character :: op_
key = 0
if(transl == 'T') key = key + 1
if(transl == 'D') key = key + 2
if(transr == 'T') key = key + 10
if(transr == 'D') key = key + 20
op_ = 'N'
if(present(op)) op_ = op
select case(key)
case(0)
call kron_mm_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(1)
call kron_mtm_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(2)
call kron_mdm_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(10)
call kron_mmt_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(11)
call kron_mtmt_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(12)
call kron_mdmt_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(20)
call kron_mmd_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(21)
call kron_mtmd_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(22)
call kron_mdmd_real_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
end select
end subroutine kron_real_real
"""
return
[docs]def kron_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R`.
**Arguments**
Mat : complex(\*, \*), inout
On exit the tensor product between the two matrices. Has to be
allocated on entry unless optional argument is set otherwise.
Matl : _complex(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
Matr : complex(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
transl : CHARACTER, in
Defines the transformation executed on the left matrix :math:`M_L`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
transr : CHARACTER, in
Defines the transformation executed on the left matrix :math:`M_R`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_complex_complex(Mat, Matl, Matr, r1, c1, r2, c2, &
transl, transr, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*r2*c1*c2), intent(inout) :: Mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: Matl
complex(KIND=rKind), dimension(r2, c2), intent(in) :: Matr
character, intent(in) :: transl, transr
character, intent(in), optional :: op
! Local variables
! ---------------
! Key for calling correct subroutine: 0 : M x M, 1 : M^t x M,
! 10 : M x M^t, 11 M^t x M^t, 2 : M^d x M, 20 : M x M^d
! 22 : M^d x M^d, 21 : M^d x M^t, 12 : M^t x M^d
integer :: key
! dublette for optional argument
character :: op_
key = 0
if(transl == 'T') key = key + 1
if(transl == 'D') key = key + 2
if(transr == 'T') key = key + 10
if(transr == 'D') key = key + 20
op_ = 'N'
if(present(op)) op_ = op
select case(key)
case(0)
call kron_mm_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(1)
call kron_mtm_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(2)
call kron_mdm_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(10)
call kron_mmt_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(11)
call kron_mtmt_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(12)
call kron_mdmt_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(20)
call kron_mmd_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(21)
call kron_mtmd_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(22)
call kron_mdmd_complex_complex(mat, matl, matr, &
r1, c1, r2, c2, op_)
end select
end subroutine kron_complex_complex
"""
return
[docs]def kron_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R`.
**Arguments**
Mat : complex(\*, \*), inout
On exit the tensor product between the two matrices. Has to be
allocated on entry unless optional argument is set otherwise.
Matl : _real(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
Matr : real(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
transl : CHARACTER, in
Defines the transformation executed on the left matrix :math:`M_L`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
transr : CHARACTER, in
Defines the transformation executed on the left matrix :math:`M_R`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_complex_real(Mat, Matl, Matr, r1, c1, r2, c2, &
transl, transr, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*r2*c1*c2), intent(inout) :: Mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: Matl
real(KIND=rKind), dimension(r2, c2), intent(in) :: Matr
character, intent(in) :: transl, transr
character, intent(in), optional :: op
! Local variables
! ---------------
! Key for calling correct subroutine: 0 : M x M, 1 : M^t x M,
! 10 : M x M^t, 11 M^t x M^t, 2 : M^d x M, 20 : M x M^d
! 22 : M^d x M^d, 21 : M^d x M^t, 12 : M^t x M^d
integer :: key
! dublette for optional argument
character :: op_
key = 0
if(transl == 'T') key = key + 1
if(transl == 'D') key = key + 2
if(transr == 'T') key = key + 10
if(transr == 'D') key = key + 20
op_ = 'N'
if(present(op)) op_ = op
select case(key)
case(0)
call kron_mm_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(1)
call kron_mtm_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(2)
call kron_mdm_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(10)
call kron_mmt_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(11)
call kron_mtmt_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(12)
call kron_mdmt_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(20)
call kron_mmd_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(21)
call kron_mtmd_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
case(22)
call kron_mdmd_complex_real(mat, matl, matr, &
r1, c1, r2, c2, op_)
end select
end subroutine kron_complex_real
"""
return
[docs]def kron_mm_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mm_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(r1*r2, c1*c2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mm_real_real
"""
return
[docs]def kron_mm_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mm_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*r2, c1*c2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mm_complex_complex
"""
return
[docs]def kron_mm_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mm_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*r2, c1*c2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mm_complex_real
"""
return
[docs]def kron_mmt_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R^t`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. Is transposed while
taking tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mmt_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(r1*c2, c1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mmt_real_real
"""
return
[docs]def kron_mmt_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R^t`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. Is transposed while
taking tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mmt_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*c2, c1*r2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mmt_complex_complex
"""
return
[docs]def kron_mmt_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L \otimes M_R^t`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. Is transposed while
taking tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mmt_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*c2, c1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mmt_complex_real
"""
return
[docs]def kron_mmd_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. Is daggered while
taking tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mmd_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(r1*c2, c1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = (mat2(j2, j1)) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ (mat2(j2, j1)) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mmd_real_real
"""
return
[docs]def kron_mmd_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. Is daggered while
taking tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mmd_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*c2, c1*r2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = conjg(mat2(j2, j1)) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ conjg(mat2(j2, j1)) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mmd_complex_complex
"""
return
[docs]def kron_mmd_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. Is daggered while
taking tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mmd_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(r1*c2, c1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = (mat2(j2, j1)) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, c1
k3 = k3 + 1
k1 = 1
k2 = r1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ (mat2(j2, j1)) * mat1(:, i2)
k1 = k1 + r1
k2 = k2 + r1
end do
end do
end do
end if
end subroutine kron_mmd_complex_real
"""
return
[docs]def kron_mtm_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L^t \otimes M_R`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. Is transposed while
taking the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtm_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(c1*r2, r1*c2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtm_real_real
"""
return
[docs]def kron_mtm_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L^t \otimes M_R`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. Is transposed while
taking the tensor product.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtm_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*r2, r1*c2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtm_complex_complex
"""
return
[docs]def kron_mtm_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices :math:`M_L^t \otimes M_R`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. Is transposed while
taking the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtm_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*r2, r1*c2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtm_complex_real
"""
return
[docs]def kron_mdm_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. Is daggered while
taking the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdm_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(c1*r2, r1*c2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdm_real_real
"""
return
[docs]def kron_mdm_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. Is daggered while
taking the tensor product.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdm_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*r2, r1*c2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * conjg(mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * conjg(mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdm_complex_complex
"""
return
[docs]def kron_mdm_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. Is daggered while
taking the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdm_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*r2, r1*c2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat2(j1, j2) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, c2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, r2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j1, j2) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdm_complex_real
"""
return
[docs]def kron_mtmt_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^t \otimes M_R^t`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtmt_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtmt_real_real
"""
return
[docs]def kron_mtmt_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^t \otimes M_R^t`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtmt_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtmt_complex_complex
"""
return
[docs]def kron_mtmt_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^t \otimes M_R^t`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtmt_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtmt_complex_real
"""
return
[docs]def kron_mdmd_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is daggered during the
subroutine.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is daggered during the
subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdmd_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = (mat2(j2, j1)) &
* (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ (mat2(j2, j1)) &
* (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdmd_real_real
"""
return
[docs]def kron_mdmd_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is daggered during the
subroutine.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is daggered during the
subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdmd_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = conjg(mat2(j2, j1)) &
* conjg(mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ conjg(mat2(j2, j1)) &
* conjg(mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdmd_complex_complex
"""
return
[docs]def kron_mdmd_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is daggered during the
subroutine.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is daggered during the
subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdmd_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = (mat2(j2, j1)) &
* (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ (mat2(j2, j1)) &
* (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdmd_complex_real
"""
return
[docs]def kron_mdmt_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R^{t}`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is daggered
during the subroutine.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is transposed
during the subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdmt_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(c1 * c2, r1 * r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdmt_real_real
"""
return
[docs]def kron_mdmt_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R^{t}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is daggered
during the subroutine.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is transposed
during the subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdmt_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1 * c2, r1 * r2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * conjg(mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * conjg(mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdmt_complex_complex
"""
return
[docs]def kron_mdmt_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{\dagger} \otimes M_R^{t}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is daggered
during the subroutine.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is transposed
during the subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdmt_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1 * c2, r1 * r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat2(j2, j1) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, k3) &
+ mat2(j2, j1) * (mat1(i2, :))
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mdmt_complex_real
"""
return
[docs]def kron_mtmd_real_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{t} \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(real)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is transposed
during the subroutine.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is daggered
during the subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtmd_real_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
real(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = (mat2(j2, j1)) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, 3) &
+ (mat2(j2, j1)) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtmd_real_real
"""
return
[docs]def kron_mtmd_complex_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{t} \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is transposed
during the subroutine.
mat2 : TYPE(complex)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is daggered
during the subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtmd_complex_complex(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = conjg(mat2(j2, j1)) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, 3) &
+ conjg(mat2(j2, j1)) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtmd_complex_complex
"""
return
[docs]def kron_mtmd_complex_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between two matrices
:math:`M_L^{t} \otimes M_R^{\dagger}`.
**Arguments**
mat : TYPE(complex)(\*, \*), inout
On exit the tensor product between the two matrices.
mat1 : TYPE(real)(\*, \*), in
Left matrix :math:`M_L` in the tensor product. It is transposed
during the subroutine.
mat2 : TYPE(real)(\*, \*), in
Left matrix :math:`M_R` in the tensor product. It is daggered
during the subroutine.
r1 : INTEGER, in
First dimension of :math:`M_L`
c1 : INTEGER, in
Second dimension of :math:`M_L`
r2 : INTEGER, in
First dimension of :math:`M_R`
c2 : INTEGER, in
Second dimension of :math:`M_R`
op : CHARACTER, in
Either 'N' for usual Kronecker product, or '+' for adding
to the result.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mtmd_complex_real(mat, mat1, mat2, r1, c1, &
r2, c2, op)
integer, intent(in) :: r1, c1, r2, c2
complex(KIND=rKind), dimension(c1*c2, r1*r2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
character, intent(in) :: op
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k2, k3
if(op == 'N') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = (mat2(j2, j1)) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
elseif(op == '+') then
k3 = 0
do j2 = 1, r2
do i2 = 1, r1
k3 = k3 + 1
k1 = 1
k2 = c1
do j1 = 1, c2
mat(k1:k2, k3) = mat(k1:k2, 3) &
+ (mat2(j2, j1)) * mat1(i2, :)
k1 = k1 + c1
k2 = k2 + c1
end do
end do
end do
end if
end subroutine kron_mtmd_complex_real
"""
return
[docs]def kron_id_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity, so
either :math:`M_{in} \otimes 1` or :math:`1 \otimes M_{in}`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
matin : real(\*, \*), in
Matrix in the tensor product with the identity
d1 : INTEGER, in
First dimension of matin.
d2 : INTEGER, in
Second dimension of matin.
did : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
id : CHARACTER, in
Defines the position of the identity matrix. Either 'L' for
:math:`1 \otimes M_{in}` or 'R' for :math:`M_{in} \otimes 1`.
trans : CHARACTER, in
Defines the transformation executed on matrix :math:`M_{in}`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_id_real(mat, matin, d1, d2, did, id, trans)
integer, intent(in) :: d1, d2, did
real(KIND=rKind), dimension(d1 * did, d2 * did), intent(inout) :: mat
real(KIND=rKind), dimension(d1, d2), intent(in) :: matin
character, intent(in) :: id, trans
! Local variables
! ---------------
! Key for calling correct subroutine: 0 : M x Id, 1 : M^t x Id,
! 10 : Id x M, 11 Id x M^t, 2 : M^d x Id, 12 Id x M^d
integer :: key
key = 0
if(trans == 'T') key = key + 1
if(trans == 'D') key = key + 2
if(id == 'L') key = key + 10
select case(key)
case(0)
call kron_mi_real(mat, matin, d1, d2, did)
case(1)
call kron_mti_real(mat, matin, d1, d2, did)
case(2)
call kron_mdi_real(mat, matin, d1, d2, did)
case(10)
call kron_im_real(mat, matin, did, d1, d2)
case(11)
call kron_imt_real(mat, matin, did, d1, d2)
case(12)
call kron_imd_real(mat, matin, did, d1, d2)
end select
end subroutine kron_id_real
"""
return
[docs]def kron_id_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity, so
either :math:`M_{in} \otimes 1` or :math:`1 \otimes M_{in}`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
matin : complex(\*, \*), in
Matrix in the tensor product with the identity
d1 : INTEGER, in
First dimension of matin.
d2 : INTEGER, in
Second dimension of matin.
did : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
id : CHARACTER, in
Defines the position of the identity matrix. Either 'L' for
:math:`1 \otimes M_{in}` or 'R' for :math:`M_{in} \otimes 1`.
trans : CHARACTER, in
Defines the transformation executed on matrix :math:`M_{in}`
when taking the tensor product. Either 'N' (none), 'T' (transposed),
or 'D' (dagger / conjugate transposed).
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_id_complex(mat, matin, d1, d2, did, id, trans)
integer, intent(in) :: d1, d2, did
complex(KIND=rKind), dimension(d1 * did, d2 * did), intent(inout) :: mat
complex(KIND=rKind), dimension(d1, d2), intent(in) :: matin
character, intent(in) :: id, trans
! Local variables
! ---------------
! Key for calling correct subroutine: 0 : M x Id, 1 : M^t x Id,
! 10 : Id x M, 11 Id x M^t, 2 : M^d x Id, 12 Id x M^d
integer :: key
key = 0
if(trans == 'T') key = key + 1
if(trans == 'D') key = key + 2
if(id == 'L') key = key + 10
select case(key)
case(0)
call kron_mi_complex(mat, matin, d1, d2, did)
case(1)
call kron_mti_complex(mat, matin, d1, d2, did)
case(2)
call kron_mdi_complex(mat, matin, d1, d2, did)
case(10)
call kron_im_complex(mat, matin, did, d1, d2)
case(11)
call kron_imt_complex(mat, matin, did, d1, d2)
case(12)
call kron_imd_complex(mat, matin, did, d1, d2)
end select
end subroutine kron_id_complex
"""
return
[docs]def kron_mi_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`M_{in} \otimes 1`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat1 : real(\*, \*), in
Matrix in the tensor product with the identity
r1 : INTEGER, in
First dimension of the matrix.
c1 : INTEGER, in
Second dimension of the matrix.
d2 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mi_real(mat, mat1, r1, c1, d2)
integer, intent(in) :: r1, c1, d2
real(KIND=rKind), dimension(r1 * d2, c1 * d2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
! Local variables
! ---------------
! for looping
integer :: i2, jj, k1, k2, k3
mat = 0.0_rKind
k3 = 0
k1 = 1
k2 = r1
do jj = 1, d2
do i2 = 1, c1
k3 = k3 + 1
mat(k1:k2, k3) = mat1(:, i2)
end do
k1 = k1 + r1
k2 = k2 + r1
end do
end subroutine kron_mi_real
"""
return
[docs]def kron_mi_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`M_{in} \otimes 1`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat1 : complex(\*, \*), in
Matrix in the tensor product with the identity
r1 : INTEGER, in
First dimension of the matrix.
c1 : INTEGER, in
Second dimension of the matrix.
d2 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mi_complex(mat, mat1, r1, c1, d2)
integer, intent(in) :: r1, c1, d2
complex(KIND=rKind), dimension(r1 * d2, c1 * d2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
! Local variables
! ---------------
! for looping
integer :: i2, jj, k1, k2, k3
mat = 0.0_rKind
k3 = 0
k1 = 1
k2 = r1
do jj = 1, d2
do i2 = 1, c1
k3 = k3 + 1
mat(k1:k2, k3) = mat1(:, i2)
end do
k1 = k1 + r1
k2 = k2 + r1
end do
end subroutine kron_mi_complex
"""
return
[docs]def kron_im_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`1 \otimes M_{in}`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat2 : real(\*, \*), in
Matrix in the tensor product with the identity
d1 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
r2 : INTEGER, in
First dimension of the matrix.
c2 : INTEGER, in
Second dimension of the matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_im_real(mat, mat2, d1, r2, c2)
integer, intent(in) :: d1, r2, c2
real(KIND=rKind), dimension(d1 * r2, d1 * c2), intent(out) :: mat
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k3
mat = 0.0_rKind
k3 = 0
do j2 = 1, c2
do i2 = 1, d1
k3 = k3 + 1
k1 = 0
do j1 = 1, r2
mat(k1 + i2, k3) = mat2(j1, j2)
k1 = k1 + d1
end do
end do
end do
end subroutine kron_im_real
"""
return
[docs]def kron_im_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`1 \otimes M_{in}`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat2 : complex(\*, \*), in
Matrix in the tensor product with the identity
d1 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
r2 : INTEGER, in
First dimension of the matrix.
c2 : INTEGER, in
Second dimension of the matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_im_complex(mat, mat2, d1, r2, c2)
integer, intent(in) :: d1, r2, c2
complex(KIND=rKind), dimension(d1 * r2, d1 * c2), intent(out) :: mat
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k3
mat = 0.0_rKind
k3 = 0
do j2 = 1, c2
do i2 = 1, d1
k3 = k3 + 1
k1 = 0
do j1 = 1, r2
mat(k1 + i2, k3) = mat2(j1, j2)
k1 = k1 + d1
end do
end do
end do
end subroutine kron_im_complex
"""
return
[docs]def kron_mti_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`M_{in}^t \otimes 1`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat1 : real(\*, \*), in
Matrix in the tensor product with the identity. Transposed while
taking the tensor product.
r1 : INTEGER, in
First dimension of the matrix.
c1 : INTEGER, in
Second dimension of the matrix.
d2 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mti_real(mat, mat1, r1, c1, d2)
integer, intent(in) :: r1, c1, d2
real(KIND=rKind), dimension(c1 * d2, r1 * d2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
! Local variables
! ---------------
! for looping
integer :: i2, jj, k1, k2, k3
mat = 0.0_rKind
k3 = 0
k1 = 1
k2 = c1
do jj = 1, d2
do i2 = 1, r1
k3 = k3 + 1
mat(k1:k2, k3) = mat1(i2, :)
end do
k1 = k1 + c1
k2 = k2 + c1
end do
end subroutine kron_mti_real
"""
return
[docs]def kron_mti_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`M_{in}^t \otimes 1`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat1 : complex(\*, \*), in
Matrix in the tensor product with the identity. Transposed while
taking the tensor product.
r1 : INTEGER, in
First dimension of the matrix.
c1 : INTEGER, in
Second dimension of the matrix.
d2 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mti_complex(mat, mat1, r1, c1, d2)
integer, intent(in) :: r1, c1, d2
complex(KIND=rKind), dimension(c1 * d2, r1 * d2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
! Local variables
! ---------------
! for looping
integer :: i2, jj, k1, k2, k3
mat = 0.0_rKind
k3 = 0
k1 = 1
k2 = c1
do jj = 1, d2
do i2 = 1, r1
k3 = k3 + 1
mat(k1:k2, k3) = mat1(i2, :)
end do
k1 = k1 + c1
k2 = k2 + c1
end do
end subroutine kron_mti_complex
"""
return
[docs]def kron_imt_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`1 \otimes M_{in}^t`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat2 : real(\*, \*), in
Matrix in the tensor product with the identity. Transposed while
taking tensor product.
d1 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
r2 : INTEGER, in
First dimension of the matrix.
c2 : INTEGER, in
Second dimension of the matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_imt_real(mat, mat2, d1, r2, c2)
integer, intent(in) :: d1, r2, c2
real(KIND=rKind), dimension(d1 * c2, d1 * r2), intent(out) :: mat
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k3
mat = 0.0_rKind
k3 = 0
do j2 = 1, r2
do i2 = 1, d1
k3 = k3 + 1
k1 = 0
do j1 = 1, c2
mat(k1 + i2, k3) = mat2(j2, j1)
k1 = k1 + d1
end do
end do
end do
end subroutine kron_imt_real
"""
return
[docs]def kron_imt_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`1 \otimes M_{in}^t`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat2 : complex(\*, \*), in
Matrix in the tensor product with the identity. Transposed while
taking tensor product.
d1 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
r2 : INTEGER, in
First dimension of the matrix.
c2 : INTEGER, in
Second dimension of the matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_imt_complex(mat, mat2, d1, r2, c2)
integer, intent(in) :: d1, r2, c2
complex(KIND=rKind), dimension(d1 * c2, d1 * r2), intent(out) :: mat
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k3
mat = 0.0_rKind
k3 = 0
do j2 = 1, r2
do i2 = 1, d1
k3 = k3 + 1
k1 = 0
do j1 = 1, c2
mat(k1 + i2, k3) = mat2(j2, j1)
k1 = k1 + d1
end do
end do
end do
end subroutine kron_imt_complex
"""
return
[docs]def kron_mdi_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`M_{in}^{\dagger} \otimes 1`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat1 : real(\*, \*), in
Matrix in the tensor product with the identity. Daggered while
taking the tensor product.
r1 : INTEGER, in
First dimension of the matrix.
c1 : INTEGER, in
Second dimension of the matrix.
d2 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdi_real(mat, mat1, r1, c1, d2)
integer, intent(in) :: r1, c1, d2
real(KIND=rKind), dimension(c1 * d2, r1 * d2), intent(out) :: mat
real(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
! Local variables
! ---------------
! for looping
integer :: i2, jj, k1, k2, k3
mat = 0.0_rKind
k3 = 0
k1 = 1
k2 = c1
do jj = 1, d2
do i2 = 1, r1
k3 = k3 + 1
mat(k1:k2, k3) = (mat1(i2, :))
end do
k1 = k1 + c1
k2 = k2 + c1
end do
end subroutine kron_mdi_real
"""
return
[docs]def kron_mdi_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`M_{in}^{\dagger} \otimes 1`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat1 : complex(\*, \*), in
Matrix in the tensor product with the identity. Daggered while
taking the tensor product.
r1 : INTEGER, in
First dimension of the matrix.
c1 : INTEGER, in
Second dimension of the matrix.
d2 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_mdi_complex(mat, mat1, r1, c1, d2)
integer, intent(in) :: r1, c1, d2
complex(KIND=rKind), dimension(c1 * d2, r1 * d2), intent(out) :: mat
complex(KIND=rKind), dimension(r1, c1), intent(in) :: mat1
! Local variables
! ---------------
! for looping
integer :: i2, jj, k1, k2, k3
mat = 0.0_rKind
k3 = 0
k1 = 1
k2 = c1
do jj = 1, d2
do i2 = 1, r1
k3 = k3 + 1
mat(k1:k2, k3) = conjg(mat1(i2, :))
end do
k1 = k1 + c1
k2 = k2 + c1
end do
end subroutine kron_mdi_complex
"""
return
[docs]def kron_imd_real():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`1 \otimes M_{in}^{\dagger}`.
**Arguments**
mat : real(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat2 : real(\*, \*), in
Matrix in the tensor product with the identity. Daggered while
taking tensor product.
d1 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
r2 : INTEGER, in
First dimension of the matrix.
c2 : INTEGER, in
Second dimension of the matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_imd_real(mat, mat2, d1, r2, c2)
integer, intent(in) :: d1, r2, c2
real(KIND=rKind), dimension(d1 * c2, d1 * r2), intent(out) :: mat
real(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k3
mat = 0.0_rKind
k3 = 0
do j2 = 1, r2
do i2 = 1, d1
k3 = k3 + 1
k1 = 0
do j1 = 1, c2
mat(k1 + i2, k3) = (mat2(j2, j1))
k1 = k1 + d1
end do
end do
end do
end subroutine kron_imd_real
"""
return
[docs]def kron_imd_complex():
"""
fortran-subroutine - May 2016 (dj)
Kronecker or tensor product between a matrix and the identity
:math:`1 \otimes M_{in}^{\dagger}`.
**Arguments**
mat : complex(\*, \*), in
On exit the tensor product between the two matrices. Has to be
allocated on entry.
mat2 : complex(\*, \*), in
Matrix in the tensor product with the identity. Daggered while
taking tensor product.
d1 : INTEGER, in
Specifies the dimension of the identity matrix. Alwyas square matrix.
r2 : INTEGER, in
First dimension of the matrix.
c2 : INTEGER, in
Second dimension of the matrix.
**Details**
(template defined in LinearAlgebra_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine kron_imd_complex(mat, mat2, d1, r2, c2)
integer, intent(in) :: d1, r2, c2
complex(KIND=rKind), dimension(d1 * c2, d1 * r2), intent(out) :: mat
complex(KIND=rKind), dimension(r2, c2), intent(in) :: mat2
! Local variables
! ---------------
! for looping
integer :: i2, j1, j2, k1, k3
mat = 0.0_rKind
k3 = 0
do j2 = 1, r2
do i2 = 1, d1
k3 = k3 + 1
k1 = 0
do j1 = 1, c2
mat(k1 + i2, k3) = conjg(mat2(j2, j1))
k1 = k1 + d1
end do
end do
end do
end subroutine kron_imd_complex
"""
return
[docs]def qr_inplace_real():
"""
fortran-subroutine - September 2014 (dj)
Compute thin QR-decomposition inplace for nxm-matrix with n >=m
**Arguments**
mat : REAL(ld_mat, \\*), inout
Matrix to be deomposed in QR with nr >= nc. On exit
matrix Q.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : REAL(nc, \\*), out
On exit triangular matrix R of QR-decomposition.
**Details**
(template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine qr_inplace_real(mat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
real(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
real(KIND=rKind), dimension(nc, nc), intent(out) :: rmat
integer, intent(out), optional :: errst
! for LAPACK
real(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
real(KIND=rKind), dimension(:), allocatable :: work
real(KIND=rKind), dimension(1) :: query
! for looping / workspace query
integer ii
!if(present(errst)) errst = 0
! ZGEQRF (query and call)
! =======================
lwork = -1
call dgeqrf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_real: query dgeqrf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call dgeqrf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_real: dgeqrf failed', &
! 6, errst=errst)
! return
!end if
! Get the triangular
! ==================
do ii = 1, nc
rmat(:ii, ii) = mat(:ii, ii)
rmat(ii+1:, ii) = 0.0_rKind
end do
! ZUNGQR (query and call)
! =======================
ii = -1
call dorgqr(nr, nc, nc, mat, ld_mat, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_real: query dorgqr failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call dorgqr(nr, nc, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_real: dorgqr failed', &
! 6, errst=errst)
! return
!end if
deallocate(work)
end subroutine qr_inplace_real
"""
return
[docs]def qr_trapez_real():
"""
fortran-subroutine - April 2016 (dj)
Compute thin QR-decomposition inplace for nxm-matrix with n < m
**Arguments**
mat : REAL(ld_mat, \\*), inout
Matrix to be deomposed in QR with nr < nc. On exit
entries are destroyed.
qmat : REAL(nr, nr), out
On exit, unitary matrix Q of QR-decomposition.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : REAL(nc, \\*), out
On exit triangular matrix R of QR-decomposition.
**Details**
(template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine qr_trapez_real(mat, qmat, nr, nc, ld_mat, rmat, errst)
real(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
real(KIND=rKind), dimension(nr, nr), intent(out) :: qmat
integer, intent(in) :: nr, nc, ld_mat
real(KIND=rKind), dimension(nr, nc), intent(out) :: rmat
integer, intent(out), optional :: errst
! for LAPACK
real(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
real(KIND=rKind), dimension(:), allocatable :: work
real(KIND=rKind), dimension(1) :: query
! for looping / workspace query
integer ii
! ZGEQRF (query and call)
! =======================
lwork = -1
call dgeqrf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_real: query dgeqrf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call dgeqrf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_real: dgeqrf failed', &
! 6, errst=errst)
! return
!end if
! Get the triangular
! ==================
do ii = 1, nr
rmat(:ii, ii) = mat(:ii, ii)
rmat(ii+1:, ii) = dzero
end do
do ii = nr + 1, nc
rmat(:, ii) = mat(:, ii)
end do
! ZUNGQR (query and call)
! =======================
ii = -1
call dorgqr(nr, nr, nr, mat, ld_mat, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_real: query dorgqr failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call dorgqr(nr, nr, nr, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_real: dorgqr failed', &
! 6, errst=errst)
! return
!end if
qmat = mat(:, :nr)
deallocate(work)
end subroutine qr_trapez_real
"""
return
[docs]def qr_inplace_complex():
"""
fortran-subroutine - September 2014 (dj)
Compute thin QR-decomposition inplace for nxm-matrix with n >=m
**Arguments**
mat : COMPLEX(ld_mat, \\*), inout
Matrix to be deomposed in QR with nr >= nc. On exit
matrix Q.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : COMPLEX(nc, \\*), out
On exit triangular matrix R of QR-decomposition.
**Details**
(template defined in LinearAlgebar_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine qr_inplace_complex(mat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
complex(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
complex(KIND=rKind), dimension(nc, nc), intent(out) :: rmat
integer, intent(out), optional :: errst
! for LAPACK
complex(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
complex(KIND=rKind), dimension(:), allocatable :: work
complex(KIND=rKind), dimension(1) :: query
! for looping / workspace query
integer ii
!if(present(errst)) errst = 0
! ZGEQRF (query and call)
! =======================
lwork = -1
call zgeqrf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_complex: query zgeqrf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call zgeqrf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_complex: zgeqrf failed', &
! 6, errst=errst)
! return
!end if
! Get the triangular
! ==================
do ii = 1, nc
rmat(:ii, ii) = mat(:ii, ii)
rmat(ii+1:, ii) = 0.0_rKind
end do
! ZUNGQR (query and call)
! =======================
ii = -1
call zungqr(nr, nc, nc, mat, ld_mat, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_complex: query zungqr failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call zungqr(nr, nc, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_inplace_complex: zungqr failed', &
! 6, errst=errst)
! return
!end if
deallocate(work)
end subroutine qr_inplace_complex
"""
return
[docs]def qr_trapez_complex():
"""
fortran-subroutine - April 2016 (dj)
Compute thin QR-decomposition inplace for nxm-matrix with n >=m
**Arguments**
mat : COMPLEX(ld_mat, \\*), inout
Matrix to be deomposed in QR with nr < nc. On exit
entries are destroyed.
qmat : COMPLEX(nr, nr), out
On exit, unitary matrix Q of QR decomposition.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : COMPLEX(nc, \\*), out
On exit triangular matrix R of QR-decomposition.
**Details**
(template defined in LinearAlgebar_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine qr_trapez_complex(mat, qmat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
complex(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
complex(KIND=rKind), dimension(nr, nr), intent(out) :: qmat
complex(KIND=rKind), dimension(nr, nc), intent(out) :: rmat
integer, intent(out), optional :: errst
! for LAPACK
complex(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
complex(KIND=rKind), dimension(:), allocatable :: work
complex(KIND=rKind), dimension(1) :: query
! for looping / workspace query
integer ii
!if(present(errst)) errst = 0
! ZGEQRF (query and call)
! =======================
lwork = -1
call zgeqrf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_complex: query zgeqrf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call zgeqrf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_complex: zgeqrf failed', &
! 6, errst=errst)
! return
!end if
! Get the triangular
! ==================
do ii = 1, nr
rmat(:ii, ii) = mat(:ii, ii)
rmat(ii+1:, ii) = zzero
end do
do ii = nr + 1, nc
rmat(:, ii) = mat(:, ii)
end do
! ZUNGQR (query and call)
! =======================
ii = -1
call zungqr(nr, nr, nr, mat, ld_mat, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_complex: query zungqr failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call zungqr(nr, nr, nr, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('qr_trapez_complex: zungqr failed', &
! 6, errst=errst)
! return
!end if
qmat = mat(:, :nr)
deallocate(work)
end subroutine qr_trapez_complex
"""
return
[docs]def rq_inplace_real():
"""
fortran-subroutine - September 2014 (dj)
Compute thin RQ-decomposition inplace for nxm-matrix with n <= m
**Arguments**
mat : REAL(ld_mat, \\*), inout
Matrix to be deomposed in RQ with nr <= nc. On exit
matrix Q.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : REAL(nc, \\*), out
On exit triangular matrix R of RQ-decomposition.
**Details**
(template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine rq_inplace_real(mat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
real(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
real(KIND=rKind), dimension(nr, nr), intent(out) :: rmat
integer, intent(out), optional :: errst
! for looping/query
integer :: ii
! for LAPACK
real(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
real(KIND=rKind), dimension(:), allocatable :: work
real(KIND=rKind), dimension(1) :: query
!if(present(errst)) errst = errst
! ZGERQF (query and run)
! ----------------------
lwork = -1
call dgerqf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_real: query dgerqf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call dgerqf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_real: dgerqf failed', &
! 6, errst=errst)
! return
!end if
! Copy triangular
! ---------------
do ii = 1, nr
rmat(1:ii, ii) = mat(1:ii, nc - nr + ii)
rmat(ii+1:nr, ii) = 0.0_rKind
end do
! ZUNGRQ (query and run)
! ----------------------
ii = -1
call dorgrq(nr, nc, nr, mat, ld_mat, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_real: query dorgrq failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call dorgrq(nr, nc, nr, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_real: dorgrq failed', &
! 6, errst=errst)
! return
!end if
deallocate(work)
end subroutine rq_inplace_real
"""
return
[docs]def rq_trapez_real():
"""
fortran-subroutine - April 2016 (dj)
Compute thin RQ-decomposition inplace for nxm-matrix with n > m
**Arguments**
mat : REAL(ld_mat, \\*), inout
Matrix to be deomposed in RQ with nr > nc. On exit
entries are destroyed.
qmat : REAL(nc, nc), out
On exit, qmat is unitary matrix Q or RQ-decomposition.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : REAL(nc, \\*), out
On exit triangular matrix R of RQ-decomposition.
**Details**
(template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine rq_trapez_real(mat, qmat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
real(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
real(KIND=rKind), dimension(nc, nc), intent(out) :: qmat
real(KIND=rKind), dimension(nr, nc), intent(out) :: rmat
integer, intent(out), optional :: errst
! for looping/query
integer :: ii
! for LAPACK
real(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
real(KIND=rKind), dimension(:), allocatable :: work
real(KIND=rKind), dimension(1) :: query
!if(present(errst)) errst = 0
! ZGERQF (query and run)
! ----------------------
lwork = -1
call dgerqf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_real: query dgerqf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call dgerqf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_real: dgerqf failed', &
! 6, errst=errst)
! return
!end if
! Copy triangular
! ---------------
do ii = 1, nc
rmat(:nr - ii + 1, nc - ii + 1) = &
mat(:nr - ii + 1, nc - ii + 1)
rmat(nr - ii + 2:, nc - ii + 1) = dzero
end do
! ZUNGRQ (query and run)
! ----------------------
qmat = mat(nr - nc + 1:nr, :nc)
ii = -1
call dorgrq(nc, nc, nc, qmat, &
nc, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_real: query dorgrq failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call dorgrq(nc, nc, nc, qmat, &
nc, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_real: dorgrq failed', &
! 6, errst=errst)
! return
!end if
deallocate(work)
end subroutine rq_trapez_real
"""
return
[docs]def rq_inplace_complex():
"""
fortran-subroutine - September 2014 (dj)
Compute thin RQ-decomposition inplace for nxm-matrix with n <= m
**Arguments**
mat : COMPLEX(ld_mat, \\*), inout
Matrix to be deomposed in RQ with nr <= nc. On exit
matrix Q.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : COMPLEX(nc, \\*), out
On exit triangular matrix R of RQ-decomposition.
**Details**
(template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine rq_inplace_complex(mat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
complex(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
complex(KIND=rKind), dimension(nr, nr), intent(out) :: rmat
integer, intent(out), optional :: errst
! for looping/query
integer :: ii
! for LAPACK
complex(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
complex(KIND=rKind), dimension(:), allocatable :: work
complex(KIND=rKind), dimension(1) :: query
!if(present(errst)) errst = 0
! ZGERQF (query and run)
! ----------------------
lwork = -1
call zgerqf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_complex: query zgerqf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call zgerqf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_complex: zgerqf failed', &
! 6, errst=errst)
! return
!end if
! Copy triangular
! ---------------
do ii = 1, nr
rmat(1:ii, ii) = mat(1:ii, nc - nr + ii)
rmat(ii+1:nr, ii) = 0.0_rKind
end do
! ZUNGRQ (query and run)
! ----------------------
ii = -1
call zungrq(nr, nc, nr, mat, ld_mat, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_complex: query zungrq failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call zungrq(nr, nc, nr, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_inplace_complex: zungrq failed', &
! 6, errst=errst)
! return
!end if
deallocate(work)
end subroutine rq_inplace_complex
"""
return
[docs]def rq_trapez_complex():
"""
fortran-subroutine - April 2016 (dj)
RQ-decomposition for the case that complex matrix has more rows than
colums.
**Arguments**
mat : COMPLEX(ld_mat, \\*), inout
Matrix to be deomposed in RQ with nr > nc. On exit
entries destroyed.
qmat : COMPLEX(nc, nc), inout
On exit, unitary matrix Q of RQ decomposition.
nr : INTEGER, in
number of rows in `mat`
nc : INTEGER, in
number of columns in `mat`
ld_mat : INTEGER, in
leading dimension of the array `mat`
rmat : COMPLEX(nc, \\*), out
On exit triangular matrix R of RQ-decomposition.
**Details**
(template defined in LinearAlgebra_template.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine rq_trapez_complex(mat, qmat, nr, nc, ld_mat, rmat, errst)
integer, intent(in) :: nr, nc, ld_mat
complex(KIND=rKind), dimension(ld_mat, nc), intent(inout) :: mat
complex(KIND=rKind), dimension(nc, nc), intent(out) :: qmat
complex(KIND=rKind), dimension(nr, nc), intent(out) :: rmat
integer, intent(out), optional :: errst
! for looping/query
integer :: ii
! for LAPACK
complex(KIND=rKind), dimension(min(nr, nc)) :: tau
integer :: lwork, info
complex(KIND=rKind), dimension(:), allocatable :: work
complex(KIND=rKind), dimension(1) :: query
!if(present(errst)) errst = 0
! ZGERQF (query and run)
! ----------------------
lwork = -1
call zgerqf(nr, nc, mat, ld_mat, tau, query, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_complex: query zgerqf failed', &
! 6, errst=errst)
! return
!end if
lwork = int(query(1))
allocate(work(lwork))
call zgerqf(nr, nc, mat, ld_mat, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_complex: zgerqf failed', &
! 6, errst=errst)
! return
!end if
! Copy triangular
! ---------------
do ii = 1, nc
rmat(:nr - ii + 1, nc - ii + 1) = &
mat(:nr - ii + 1, nc - ii + 1)
rmat(nr - ii + 2:, nc - ii + 1) = zzero
end do
! ZUNGRQ (query and run)
! ----------------------
qmat = mat(nr - nc + 1:nr, :nc)
ii = -1
call zungrq(nc, nc, nc, qmat, &
nc, tau, query, ii, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_complex: query zungrq failed', &
! 6, errst=errst)
! return
!end if
ii = int(query(1))
if(ii > lwork) then
lwork = ii
deallocate(work)
allocate(work(lwork))
end if
call zungrq(nc, nc, nc, qmat, &
nc, tau, work, lwork, info)
!if(info /= 0) then
! errst = raise_error('rq_trapez_complex: query zungrq failed', &
! 6, errst=errst)
! return
!end if
deallocate(work)
end subroutine rq_trapez_complex
"""
return
[docs]def svd_elem_real():
"""
fortran-subroutine - September 2014 (dj)
Perform a SVD without truncation on real arrays with DGESVD
**Arguments**
uleft : REAL(\\*, \\*), out
on exit the left unitary matrix
sing : REAL(\\*), out
contains on exit the singular values
vright : REAL(\\*,\\*), out
on exit the right unitary matrix
matin : REAL(\\*,\\*), inout
contains the matrix to be decomposed in a SVD. On exit the
entries are destroyed.
a_row : INTEGER, in
number of rows in matrix `matin`
a_col : INTEGER, in
number of columns in `matin`
info : INTEGER, inout
exit status of ZGESDD (0: sucessful, <0 : wrong argument,
>0 : not converged)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine svd_elem_real(uleft, sing, vright, matin, a_row, a_col, info)
integer, intent(in) :: a_row, a_col
real(KIND=rKind), dimension(a_row, a_col), intent(inout) :: matin
real(KIND=rKind), dimension(a_row, min(a_row, a_col)), &
intent(out) :: uleft
real(KIND=rKind), dimension(min(a_row, a_col), a_col), &
intent(out) :: vright
real(KIND=rKind), dimension(min(a_row, a_col)), intent(out) :: sing
integer, intent(inout) :: info
! number of singular values
integer :: s_dim
! for LAPACK workspace query
integer :: lwork
real(KIND=rKind), dimension(:), allocatable :: work
real(KIND=rKind), dimension(1) :: query
! Get the missing dimension not contained in the dummy arguments
s_dim = min(a_row, a_col)
call dgesvd('S', 'S', a_row, a_col, matin, a_row, sing, uleft, a_row, &
vright, s_dim, query, -1, info)
lwork = int(query(1))
allocate(work(lwork))
call dgesvd('S', 'S', a_row, a_col, matin, a_row, sing, uleft, a_row, &
vright, s_dim, work, lwork, info)
deallocate(work)
end subroutine svd_elem_real
"""
return
[docs]def svdd_elem_real():
"""
fortran-subroutine - September 2014 (dj)
Perform an SVD without truncation on real arrays with DGESDD.
**Arguments**
u : REAL(\\*, \\*), out
on exit the left unitary matrix
s : REAL(\\*), out
contains on exit the singular values
v : REAL(\\*,\\*), out
on exit the right unitary matrix
a : REAL(\\*,\\*), inout
contains the matrix to be decomposed in a SVD. On exit the
entries are destroyed.
a_row : INTEGER, in
number of rows in matrix a
a_col : INTEGER, in
number of columns in a
jobz : CHARACTER, in
'O' : either U or VT are written in input matrix a
'S' : economic SVD without overwritting
'A' : all rows/columns in U and VT are calculated
'N' : only singular values
further see LAPACK description
info : INTEGER, inout
exit status of ZGESDD (0: sucessful, <0: wrong argument,
>0: not converged)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine svdd_elem_real(u, s, v, a, a_row, a_col, jobz, info)
integer, intent(in) :: a_row, a_col
real(KIND=rKind), dimension(a_row, a_col), intent(inout) :: a
real(KIND=rKind), dimension(a_row, min(a_row, a_col)), intent(out) :: u
real(KIND=rKind), dimension(min(a_row, a_col), a_col), intent(out) :: v
real(KIND=rKind), dimension(min(a_row, a_col)), intent(out) :: s
character, intent(in) :: jobz
integer, intent(inout) :: info
! number of singular values
integer :: s_dim
! Variables/Workspaces for ZGESDD
integer :: lwork
real(KIND=rKind), dimension(:), allocatable :: work
integer, dimension(8 * min(a_row, a_col)) :: iwork
real(KIND=rKind), dimension(1) :: query
! Get the missing dimension not contained in the dummy arguments
s_dim = min(a_row, a_col)
! Workspace query and allocation of work
call dgesdd(jobz, a_row, a_col, a, a_row, s, u, a_row, v, s_dim, &
query, -1, iwork, info)
lwork = int(query(1))
allocate(work(lwork))
! Run SVD
call dgesdd(jobz, a_row, a_col, a, a_row, s, u, a_row, v, s_dim, &
work, lwork, iwork, info)
deallocate(work)
end subroutine svdd_elem_real
"""
return
[docs]def svd_elem_complex():
"""
fortran-subroutine - September 2014 (dj)
Perform a SVD without truncation on complex arrays with ZGESVD
**Arguments**
uleft : COMPLEX(\\*, \\*), out
on exit the left unitary matrix
sing : REAL(\\*), out
contains on exit the singular values
vright : COMPLEX(\\*,\\*), out
on exit the right unitary matrix
matin : COMPLEX(\\*,\\*), inout
contains the matrix to be decomposed in a SVD. On exit the
entries are destroyed.
a_row : INTEGER, in
number of rows in matrix `matin`
a_col : INTEGER, in
number of columns in `matin`
info : INTEGER, inout
exit status of ZGESDD (0: sucessful, <0: wrong argument,
>0: not converged)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine svd_elem_complex(uleft, sing, vright, matin, a_row, a_col, info)
integer, intent(in) :: a_row, a_col
complex(KIND=rKind), dimension(a_row, a_col), intent(inout) :: matin
complex(KIND=rKind), dimension(a_row, min(a_row, a_col)), &
intent(out) :: uleft
complex(KIND=rKind), dimension(min(a_row, a_col), a_col), &
intent(out) :: vright
real(KIND=rKind), dimension(min(a_row, a_col)), intent(out) :: sing
integer, intent(inout) :: info
! number of singular values
integer :: s_dim
! for LAPACK workspace query
integer :: lwork
complex(KIND=rKind), dimension(:), allocatable :: work
complex(KIND=rKind), dimension(1) :: query
! workspace
real(KIND=rKIND), dimension(:), allocatable :: rwork
! Get the missing dimension not contained in the dummy arguments
s_dim = min(a_row, a_col)
allocate(rwork(5 * s_dim))
call zgesvd('S', 'S', a_row, a_col, matin, a_row, sing, uleft, a_row, &
vright, s_dim, query, -1, rwork, info)
lwork = int(query(1))
allocate(work(lwork))
call zgesvd('S', 'S', a_row, a_col, matin, a_row, sing, uleft, a_row, &
vright, s_dim, work, lwork, rwork, info)
deallocate(rwork, work)
end subroutine svd_elem_complex
"""
return
[docs]def svdd_elem_complex():
"""
fortran-subroutine - September 2014 (dj)
Perform a SVD without truncation on complex arrays with ZGESDD
**Arguments**
u : COMPLEX(\\*, \\*), out
on exit the left unitary matrix
s : REAL(\\*), out
contains on exit the singular values
v : COMPLEX(\\*,\\*), out
on exit the right unitary matrix
a : COMPLEX(\\*,\\*), inout
contains the matrix to be decomposed in a SVD. On exit the
entries are destroyed.
a_row : INTEGER, in
number of rows in matrix a
a_col : INTEGER, in
number of columns in a
jobz : CHARACTER, in
'O' : either U or VT are written in input matrix a
'S' : economic SVD without overwritting
'A' : all rows/columns in U and VT are calculated
'N' : only singular values
further see LAPACK description
info : INTEGER, inout
exit status of ZGESDD (0: sucessful, <0: wrong argument,
>0: not converged)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine svdd_elem_complex(u, s, v, a, a_row, a_col, jobz, info)
integer, intent(in) :: a_row, a_col
complex(KIND=rKind), dimension(a_row,*), intent(inout) :: a
complex(KIND=rKind), dimension(a_row, min(a_row, a_col)), intent(out) :: u
complex(KIND=rKind), dimension(min(a_row, a_col), a_col), intent(out) :: v
real(KIND=rKind), dimension(min(a_row, a_col)), intent(out) :: s
character, intent(in) :: jobz
integer, intent(inout) :: info
! number of singular values
integer :: s_dim
! Variables/Workspaces for ZGESDD
integer :: lwork
integer, dimension(8 * min(a_row, a_col)) :: iwork
real(KIND=rKind), dimension(5 * min(a_col, a_row)**2 + &
5 * min(a_col, a_row)) :: rwork
complex(KIND=rKind), dimension(1) :: query
complex(KIND=rKind), dimension(:), allocatable :: work
! Get the missing dimension not contained in the dummy arguments
s_dim = min(a_row, a_col)
! Workspace query and allocation of work
call zgesdd(jobz, a_row, a_col, a, a_row, s, u, a_row, v, s_dim, &
query, -1, rwork, iwork, info)
lwork = int(query(1))
allocate(work(lwork))
! Run SVD
call zgesdd(jobz, a_row, a_col, a, a_row, s, u, a_row, v, s_dim, &
work, lwork, rwork, iwork, info)
deallocate(work)
end subroutine svdd_elem_complex
"""
return
[docs]def trace_rho_x_mat_real_real():
"""
fortran-function - November 2016 (dj)
Trace of the matrix-matrix multiplication of a density matrix and
a matrix.
**Arguments**
rho : RHO_TYOE(\*, \*), in
Density matrix. Must be hermitian since this property is used.
mat : real(\*, \*), in
Operator matrix.
dim : INTEGER, in
Dimension of the two square matrices.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
function trace_rho_x_mat_real_real(rho, mat, dim) result(tr)
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(in) :: rho
real(KIND=rKind), dimension(dim, dim), intent(in) :: mat
real(KIND=rKind) :: tr
! Local variables
! ---------------
! for looping
integer :: ii
tr = dzero
do ii = 1, dim
tr = tr + sum((rho(:, ii)) * mat(:, ii))
!tr = tr + sum(rho(ii, :) * mat(:, ii))
end do
end function trace_rho_x_mat_real_real
"""
return
[docs]def trace_rho_x_mat_complex_real():
"""
fortran-function - November 2016 (dj)
Trace of the matrix-matrix multiplication of a density matrix and
a matrix.
**Arguments**
rho : RHO_TYOE(\*, \*), in
Density matrix. Must be hermitian since this property is used.
mat : real(\*, \*), in
Operator matrix.
dim : INTEGER, in
Dimension of the two square matrices.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
function trace_rho_x_mat_complex_real(rho, mat, dim) result(tr)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(in) :: rho
real(KIND=rKind), dimension(dim, dim), intent(in) :: mat
complex(KIND=rKind) :: tr
! Local variables
! ---------------
! for looping
integer :: ii
tr = zzero
do ii = 1, dim
tr = tr + sum(conjg(rho(:, ii)) * mat(:, ii))
!tr = tr + sum(rho(ii, :) * mat(:, ii))
end do
end function trace_rho_x_mat_complex_real
"""
return
[docs]def trace_rho_x_mat_complex_complex():
"""
fortran-function - November 2016 (dj)
Trace of the matrix-matrix multiplication of a density matrix and
a matrix.
**Arguments**
rho : RHO_TYOE(\*, \*), in
Density matrix. Must be hermitian since this property is used.
mat : complex(\*, \*), in
Operator matrix.
dim : INTEGER, in
Dimension of the two square matrices.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
function trace_rho_x_mat_complex_complex(rho, mat, dim) result(tr)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(in) :: rho
complex(KIND=rKind), dimension(dim, dim), intent(in) :: mat
complex(KIND=rKind) :: tr
! Local variables
! ---------------
! for looping
integer :: ii
tr = zzero
do ii = 1, dim
tr = tr + sum(conjg(rho(:, ii)) * mat(:, ii))
!tr = tr + sum(rho(ii, :) * mat(:, ii))
end do
end function trace_rho_x_mat_complex_complex
"""
return
[docs]def tracematmul_real():
"""
fortran-function - November 2017 (dj)
Calculate the trace of a matrix multiplication without contructing
the intermediate matrix.
**Arguments**
matl : real(\*, \*), in
Left matrix in the multiplication.
matr : real(\*, \*), in
Right matrix in the multuplication.
dim : INTEGER, in
Dimension of the square matrices (right now we restrict the
function to square matrices).
transl : CHARACTER, in
Transformation of the left matrix. 'N' for no transformation,
'T' for transposed, and 'C' for the hermitian conjugated.
transr : CHARACTER, in
Transformation of the right matrix. 'N' for no transformation,
'T' for transposed, and 'C' for the hermitian conjugated.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
function tracematmul_real(matl, matr, dim, transl, transr, &
errst) result(tr)
integer, intent(in) :: dim
real(KIND=rKind), dimension(dim, dim), intent(in) :: matl, matr
character, intent(in), optional :: transl, transr
integer, intent(out), optional :: errst
real(KIND=rKind) :: tr
! Local variables
! ---------------
! for looping
integer :: ii, jj
! dupletes for optional arguments
character :: transl_, transr_
!if(present(errst)) errst = 0
transl_ = 'N'
if(present(transl)) transl_ = transl
transr_ = 'N'
if(present(transr)) transr_ = transr
tr = 0.0_rKind
if((transl_ == 'N') .and. (transr_ == 'N')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + matl(ii, jj) * matr(jj, ii)
end do
end do
elseif((transl_ == 'N') .and. (transr_ == 'T')) then
do jj = 1, dim
tr = tr + sum(matl(:, jj) * matr(:, jj))
end do
elseif((transl_ == 'N') .and. (transr_ == 'C')) then
do jj = 1, dim
tr = tr + sum(matl(:, jj) * (matr(:, jj)))
end do
elseif((transl_ == 'T') .and. (transr_ == 'N')) then
do ii = 1, dim
tr = tr + sum(matl(:, ii) * matr(:, ii))
end do
elseif((transl_ == 'C') .and. (transr_ == 'N')) then
do ii = 1, dim
tr = tr + sum((matl(:, ii)) * matr(:, ii))
end do
elseif((transl_ == 'T') .and. (transr_ == 'T')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + matl(jj, ii) * matr(ii, jj)
end do
end do
elseif((transl_ == 'T') .and. (transr_ == 'C')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + matl(jj, ii) * (matr(ii, jj))
end do
end do
elseif((transl_ == 'C') .and. (transr_ == 'T')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + (matl(jj, ii)) * matr(ii, jj)
end do
end do
elseif((transl_ == 'C') .and. (transr_ == 'C')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + (matl(jj, ii)) * (matr(ii, jj))
end do
end do
else
end if
end function tracematmul_real
"""
return
[docs]def tracematmul_complex():
"""
fortran-function - November 2017 (dj)
Calculate the trace of a matrix multiplication without contructing
the intermediate matrix.
**Arguments**
matl : complex(\*, \*), in
Left matrix in the multiplication.
matr : complex(\*, \*), in
Right matrix in the multuplication.
dim : INTEGER, in
Dimension of the square matrices (right now we restrict the
function to square matrices).
transl : CHARACTER, in
Transformation of the left matrix. 'N' for no transformation,
'T' for transposed, and 'C' for the hermitian conjugated.
transr : CHARACTER, in
Transformation of the right matrix. 'N' for no transformation,
'T' for transposed, and 'C' for the hermitian conjugated.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
function tracematmul_complex(matl, matr, dim, transl, transr, &
errst) result(tr)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(in) :: matl, matr
character, intent(in), optional :: transl, transr
integer, intent(out), optional :: errst
complex(KIND=rKind) :: tr
! Local variables
! ---------------
! for looping
integer :: ii, jj
! dupletes for optional arguments
character :: transl_, transr_
!if(present(errst)) errst = 0
transl_ = 'N'
if(present(transl)) transl_ = transl
transr_ = 'N'
if(present(transr)) transr_ = transr
tr = 0.0_rKind
if((transl_ == 'N') .and. (transr_ == 'N')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + matl(ii, jj) * matr(jj, ii)
end do
end do
elseif((transl_ == 'N') .and. (transr_ == 'T')) then
do jj = 1, dim
tr = tr + sum(matl(:, jj) * matr(:, jj))
end do
elseif((transl_ == 'N') .and. (transr_ == 'C')) then
do jj = 1, dim
tr = tr + sum(matl(:, jj) * conjg(matr(:, jj)))
end do
elseif((transl_ == 'T') .and. (transr_ == 'N')) then
do ii = 1, dim
tr = tr + sum(matl(:, ii) * matr(:, ii))
end do
elseif((transl_ == 'C') .and. (transr_ == 'N')) then
do ii = 1, dim
tr = tr + sum(conjg(matl(:, ii)) * matr(:, ii))
end do
elseif((transl_ == 'T') .and. (transr_ == 'T')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + matl(jj, ii) * matr(ii, jj)
end do
end do
elseif((transl_ == 'T') .and. (transr_ == 'C')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + matl(jj, ii) * conjg(matr(ii, jj))
end do
end do
elseif((transl_ == 'C') .and. (transr_ == 'T')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + conjg(matl(jj, ii)) * matr(ii, jj)
end do
end do
elseif((transl_ == 'C') .and. (transr_ == 'C')) then
do ii = 1, dim
do jj = 1, dim
tr = tr + conjg(matl(jj, ii)) * conjg(matr(ii, jj))
end do
end do
else
end if
end function tracematmul_complex
"""
return