LinearAlgebra

Fortran module LinearAlgebra: May 2017 (dj)

Containing basic linear algebra operators defined for arrays.

Authors

• 4. Jaschke

• 13. 12. 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.

LinearAlgebra_f90.diag_times_mat_real_real()

fortran-subroutine - November 2016 (dj) Multiplication of diagonal matrix with matrix.

Arguments

matreal(nrow, ncol), inout

Matrix to be multiplied with diagonal matrix in vector form.

diagreal(nrow), in

Diagonal matrix in vector form.

nrowINTEGER, in

Number of rows in matrix and length of vector for diagonal matrix.

ncolINTEGER, in

Number of columns in matrix.

Source Code

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

LinearAlgebra_f90.diag_times_mat_complex_real()

fortran-subroutine - November 2016 (dj) Multiplication of diagonal matrix with matrix.

Arguments

matcomplex(nrow, ncol), inout

Matrix to be multiplied with diagonal matrix in vector form.

diagreal(nrow), in

Diagonal matrix in vector form.

nrowINTEGER, in

Number of rows in matrix and length of vector for diagonal matrix.

ncolINTEGER, in

Number of columns in matrix.

Source Code

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

LinearAlgebra_f90.diag_times_mat_complex_complex()

fortran-subroutine - November 2016 (dj) Multiplication of diagonal matrix with matrix.

Arguments

matcomplex(nrow, ncol), inout

Matrix to be multiplied with diagonal matrix in vector form.

diagcomplex(nrow), in

Diagonal matrix in vector form.

nrowINTEGER, in

Number of rows in matrix and length of vector for diagonal matrix.

ncolINTEGER, in

Number of columns in matrix.

Source Code

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

LinearAlgebra_f90.mat_times_diag_real_real()

fortran-subroutine - November 2016 (dj) Multiplication of matrix with diagonal matrix.

Arguments

matreal(nrow, ncol), inout

Matrix to be multiplied with diagonal matrix in vector form.

diagreal(nrow), in

Diagonal matrix in vector form.

nrowINTEGER, in

Number of rows in matrix.

ncolINTEGER, in

Number of columns in matrix and length of vector for diagonal matrix.

Source Code

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

LinearAlgebra_f90.mat_times_diag_complex_real()

fortran-subroutine - November 2016 (dj) Multiplication of matrix with diagonal matrix.

Arguments

matcomplex(nrow, ncol), inout

Matrix to be multiplied with diagonal matrix in vector form.

diagreal(nrow), in

Diagonal matrix in vector form.

nrowINTEGER, in

Number of rows in matrix.

ncolINTEGER, in

Number of columns in matrix and length of vector for diagonal matrix.

Source Code

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

LinearAlgebra_f90.mat_times_diag_complex_complex()

fortran-subroutine - November 2016 (dj) Multiplication of matrix with diagonal matrix.

Arguments

matcomplex(nrow, ncol), inout

Matrix to be multiplied with diagonal matrix in vector form.

diagcomplex(nrow), in

Diagonal matrix in vector form.

nrowINTEGER, in

Number of rows in matrix.

ncolINTEGER, in

Number of columns in matrix and length of vector for diagonal matrix.

Source Code

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

LinearAlgebra_f90.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

tensreal(n1, n2, n3), inout

Rank three tensor to be contracted with diagonal matrix in vector form over second index.

diagreal(n2), in

Diagonal matrix in vector form.

n1INTEGER, in

Dimension of the first index of the tensor.

n2INTEGER, in

Dimension of the second index of the tensor and the diagonal matrix.

n3INTEGER, in

Dimension of the third index of the tensor.

Source Code

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

LinearAlgebra_f90.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

tenscomplex(n1, n2, n3), inout

Rank three tensor to be contracted with diagonal matrix in vector form over second index.

diagreal(n2), in

Diagonal matrix in vector form.

n1INTEGER, in

Dimension of the first index of the tensor.

n2INTEGER, in

Dimension of the second index of the tensor and the diagonal matrix.

n3INTEGER, in

Dimension of the third index of the tensor.

Source Code

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

LinearAlgebra_f90.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

tenscomplex(n1, n2, n3), inout

Rank three tensor to be contracted with diagonal matrix in vector form over second index.

diagcomplex(n2), in

Diagonal matrix in vector form.

n1INTEGER, in

Dimension of the first index of the tensor.

n2INTEGER, in

Dimension of the second index of the tensor and the diagonal matrix.

n3INTEGER, in

Dimension of the third index of the tensor.

Source Code

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

LinearAlgebra_f90.eigd_symm_real()

fortran-subroutine - April 2016 (dj) Diagonalize symmetric mata using LAPACK’s DSYEVD. Return the eigenvectors in mata, eigenvalues in ev.

Arguments

mataREAL(,), inout

Diagonalize this matrix. On exit overwritten with eigenvectors.

evREAL(*), out

Eigenvalues of diagonlization are stored in this vector.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.eigd_symm_complex()

fortran-subroutine - April 2016 (dj) Diagonalize hermitian mata using LAPACK’s ZHEEV. Return the eigenvectors in mata, eigenvalues in ev.

Arguments

mataCOMPLEX(,), inout

Diagonalize this matrix. On exit overwritten with eigenvectors.

evREAL(*), out

Eigenvalues of diagonlization are stored in this vector.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.eigd_real()

fortran-subroutine - November 2014 (dj) Diagonalization of a general real matrix

Arguments

matREAL(,), inout

on input containing the matrix to be diagonalized on exit right eigenvectors of the matrix

eivalCOMPLEX(*), out

containing on exit the (possibly complex) eigenvalues

dimINTEGER, in

dimension of the eigenvalue problem.

matinvREAL(,), 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

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

LinearAlgebra_f90.eigd_complex()

fortran-subroutine - November 2014 (dj) Diagonalization of a general complex matrix

Arguments

matCOMPLEX(,), inout

on input containing the matrix to be diagonalized on exit right eigenvectors of the matrix

eivalCOMPLEX(*), out

containing on exit the (possibly complex) eigenvalues

dimINTEGER, in

dimension of the eigenvalue problem.

matinvCOMPLEX(,), 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

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

LinearAlgebra_f90.eigd_tri_symm_real()

fortran-subroutine - April 2016 (dj) Diagonalize symmetric tridiagonal matrix with LAPACKs DSTEQR.

Arguments

mataREAL(*, *), out

Contains on exit the eigenvectors of the symmetric tridiagonal matrix.

diagREAL(*), inout

On entry the diagonal elements of the symmetric matrix, on exit the eigenvalues in ascending order.

supdiagREAL(*), inout

Super-diagonal of the symmetric matrix. Destroyed on exit.

dimINTEGER, in

Dimension of the square matrix

Details

(Template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.eigd_tri_symm_complex()

fortran-subroutine - April 2016 (dj) Diagonalize hermitian tridiagonal matrix using LAPACKS dense method ZHEEV.

Arguments

mataCOMPLEX(*, *), out

Contains on exit the eigenvectors of the symmetric tridiagonal matrix.

diagREAL(*), inout

On entry the diagonal elements of the hermitian matrix, on exit the eigenvalues in ascending order.

supdiagCOMPLEX(*), in

Super-diagonal of the hermitian matrix.

dimINTEGER, in

Dimension of the square matrix

Details

(Template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.eigd_symm_tri_first_real()

fortran-subroutine - July 2017 (dj, updated) Get the first eigenvalue and vector of a tridiagonal and real matrix.

Arguments

diagsREAL(*), inout

The diagonal entries in the tridiagonal matrix.

odiagsREAL(*), inout

The first off-diagonal entries in the triadiagonal matrix.

eigvalREAL, out

The first eigenvalue.

eigvecREAL(*), inout

The first eigenvector.

Source Code

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

LinearAlgebra_f90.eigsvd_lu_real()

fortran-subroutine - Eigenvalue decomposition as singular value decomposition. ?? (dj, updated)

Arguments

umreal(a_row, a_row), out

Unitary matrix for decomposition.

singreal(a_row), out

Singular values (squared?)

rmreal(a_row, a_col), out

Right matrix in decomposition, i.e., singular values times right unitary matrix.

matinreal(a_row, a_col), in

Input matrix to be decomposed.

a_rowINTEGER, in

Number of row in the original matrix.

a_colINTEGER, in

Number of columns in the original matrix.

Source Code

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

LinearAlgebra_f90.eigsvd_lu_complex()

fortran-subroutine - Eigenvalue decomposition as singular value decomposition. ?? (dj, updated)

Arguments

umcomplex(a_row, a_row), out

Unitary matrix for decomposition.

singreal(a_row), out

Singular values (squared?)

rmcomplex(a_row, a_col), out

Right matrix in decomposition, i.e., singular values times right unitary matrix.

matincomplex(a_row, a_col), in

Input matrix to be decomposed.

a_rowINTEGER, in

Number of row in the original matrix.

a_colINTEGER, in

Number of columns in the original matrix.

Source Code

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

LinearAlgebra_f90.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

umreal(a_row, a_row), out

Unitary matrix on the left of the decomposition on exit.

singreal(a_row), out

Singular values squared on exit.

matinreal(a_row, a_col), in

Input matrix to be decomposed.

a_rowINTEGER, in

Number of rows in matrix to be decomposed

a_colINTEGER, in

Number of columns in matrix to be decomposed.

Source Code

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

LinearAlgebra_f90.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

umcomplex(a_row, a_row), out

Unitary matrix on the left of the decomposition on exit.

singreal(a_row), out

Singular values squared on exit.

matincomplex(a_row, a_col), in

Input matrix to be decomposed.

a_rowINTEGER, in

Number of rows in matrix to be decomposed

a_colINTEGER, in

Number of columns in matrix to be decomposed.

Source Code

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

LinearAlgebra_f90.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

umreal(a_row, a_row), inout

Unitary matrix, modified in case of truncation.

rtreal(nkeep, a_col), out

On exit, singular values times right unitary matrix.

matinreal(a_row, a_col), in

Original matrix to be decomposed.

a_rowINTEGER, in

Number of rows in original matrix.

a_colINTEGER, in

Number of columns in original matrix.

nkeepINTEGER, in

Number of singular values to be kept.

renormREAL, OPTIONAL, in

Renormalization for decomposition (or scaling). Default to 1.

Source Code

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

LinearAlgebra_f90.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

umcomplex(a_row, a_row), inout

Unitary matrix, modified in case of truncation.

rtcomplex(nkeep, a_col), out

On exit, singular values times right unitary matrix.

matincomplex(a_row, a_col), in

Original matrix to be decomposed.

a_rowINTEGER, in

Number of rows in original matrix.

a_colINTEGER, in

Number of columns in original matrix.

nkeepINTEGER, in

Number of singular values to be kept.

renormREAL, OPTIONAL, in

Renormalization for decomposition (or scaling). Default to 1.

Source Code

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

LinearAlgebra_f90.eigsvd_ru_step2_real()

fortran-subroutine - Second step of eigenvalue decomposition as singular value decomposition. The unitary tensor is on the right.

Arguments

umreal(*, *), in

Contains on exit the eigenvectors.

Source Code

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

LinearAlgebra_f90.eigsvd_ru_step2_complex()

fortran-subroutine - Second step of eigenvalue decomposition as singular value decomposition. The unitary tensor is on the right.

Arguments

umcomplex(*, *), in

Contains on exit the eigenvectors.

Source Code

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

LinearAlgebra_f90.expm_real_real()

fortran-subroutine - April 2016 (dj) Take the exponential of a general matrix returning a complex matrix.

Arguments

mexpcomplex(dim, dim), in

Exponential of the matrix sc * mat

screal, in

Additional scalar inside exp-function.

matreal(dim, dim), inout

On entry containing the matrix to be exponentiated. On exit destroyed (and filled with right eigenvectors).

dimINTEGER, in

Dimension of the square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expm_complex_real()

fortran-subroutine - April 2016 (dj) Take the exponential of a general matrix returning a complex matrix.

Arguments

mexpcomplex(dim, dim), in

Exponential of the matrix sc * mat

sccomplex, in

Additional scalar inside exp-function.

matreal(dim, dim), inout

On entry containing the matrix to be exponentiated. On exit destroyed (and filled with right eigenvectors).

dimINTEGER, in

Dimension of the square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expm_real_complex()

fortran-subroutine - April 2016 (dj) Take the exponential of a general matrix returning a complex matrix.

Arguments

mexpcomplex(dim, dim), in

Exponential of the matrix sc * mat

screal, in

Additional scalar inside exp-function.

matcomplex(dim, dim), inout

On entry containing the matrix to be exponentiated. On exit destroyed (and filled with right eigenvectors).

dimINTEGER, in

Dimension of the square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expm_complex_complex()

fortran-subroutine - April 2016 (dj) Take the exponential of a general matrix returning a complex matrix.

Arguments

mexpcomplex(dim, dim), in

Exponential of the matrix sc * mat

sccomplex, in

Additional scalar inside exp-function.

matcomplex(dim, dim), inout

On entry containing the matrix to be exponentiated. On exit destroyed (and filled with right eigenvectors).

dimINTEGER, in

Dimension of the square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_real_real()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc * math) for hermitian math.

Arguments

mexpreal(dim, dim), out

Contains on exit exp(- eye * sc * math)

screal, in

Weight in the exponential.

matreal(*, *), inout

Take exponential of this hermitian matrix. On exit, overwritten with eigenvectors.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_complex_real()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc * math) for hermitian math.

Arguments

mexpcomplex(dim, dim), out

Contains on exit exp(- eye * sc * math)

sccomplex, in

Weight in the exponential.

matreal(*, *), inout

Take exponential of this hermitian matrix. On exit, overwritten with eigenvectors.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_real_complex()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc * math) for hermitian math.

Arguments

mexpcomplex(dim, dim), out

Contains on exit exp(- eye * sc * math)

screal, in

Weight in the exponential.

matcomplex(*, *), inout

Take exponential of this hermitian matrix. On exit, overwritten with eigenvectors.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_complex_complex()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc * math) for hermitian math.

Arguments

mexpcomplex(dim, dim), out

Contains on exit exp(- eye * sc * math)

sccomplex, in

Weight in the exponential.

matcomplex(*, *), inout

Take exponential of this hermitian matrix. On exit, overwritten with eigenvectors.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_tri_real_real()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.

Arguments

mexpreal(dim, dim), out

Contains on exit exp(- eye * sc * mat)

screal, in

Weight in the exponential.

diagREAL(*), in

Diagonal entries of the matrix to be exponentiated.

supdiagreal(*), in

Super-diagonal entries of the matrix to be exponentiated.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_tri_complex_real()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.

Arguments

mexpcomplex(dim, dim), out

Contains on exit exp(- eye * sc * mat)

sccomplex, in

Weight in the exponential.

diagREAL(*), in

Diagonal entries of the matrix to be exponentiated.

supdiagreal(*), in

Super-diagonal entries of the matrix to be exponentiated.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_tri_real_complex()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.

Arguments

mexpcomplex(dim, dim), out

Contains on exit exp(- eye * sc * mat)

screal, in

Weight in the exponential.

diagREAL(*), in

Diagonal entries of the matrix to be exponentiated.

supdiagcomplex(*), in

Super-diagonal entries of the matrix to be exponentiated.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.expmh_tri_complex_complex()

fortran-subroutine - April 2016 (dj, updated) Exponentiate U=exp(sc* math) for hermitian tridiagonal matrix math.

Arguments

mexpcomplex(dim, dim), out

Contains on exit exp(- eye * sc * mat)

sccomplex, in

Weight in the exponential.

diagREAL(*), in

Diagonal entries of the matrix to be exponentiated.

supdiagcomplex(*), in

Super-diagonal entries of the matrix to be exponentiated.

dimINTEGER, in

Dimension of the square matrix.

Details

(Template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

Matreal(*, *), 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 in the tensor product.

Matrreal(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

translCHARACTER, in

Defines the transformation executed on the left matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

transrCHARACTER, in

Defines the transformation executed on the left matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

Matcomplex(*, *), 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 in the tensor product.

Matrcomplex(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

translCHARACTER, in

Defines the transformation executed on the left matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

transrCHARACTER, in

Defines the transformation executed on the left matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

Matcomplex(*, *), 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 in the tensor product.

Matrreal(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

translCHARACTER, in

Defines the transformation executed on the left matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

transrCHARACTER, in

Defines the transformation executed on the left matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mm_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mm_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mm_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mmt_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. Is transposed while taking tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mmt_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product. Is transposed while taking tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mmt_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. Is transposed while taking tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mmd_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. Is daggered while taking tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mmd_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product. Is daggered while taking tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mmd_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. Is daggered while taking tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtm_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. Is transposed while taking the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtm_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product. Is transposed while taking the tensor product.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtm_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. Is transposed while taking the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdm_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. Is daggered while taking the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdm_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product. Is daggered while taking the tensor product.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdm_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. Is daggered while taking the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtmt_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtmt_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtmt_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdmd_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdmd_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdmd_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdmt_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. It is transposed during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdmt_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product. It is transposed during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdmt_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. It is transposed during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtmd_real_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(real)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. It is transposed during the subroutine.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtmd_complex_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(complex)(*, *), in

Left matrix in the tensor product. It is transposed during the subroutine.

mat2TYPE(complex)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mtmd_complex_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between two matrices .

Arguments

matTYPE(complex)(*, *), inout

On exit the tensor product between the two matrices.

mat1TYPE(real)(*, *), in

Left matrix in the tensor product. It is transposed during the subroutine.

mat2TYPE(real)(*, *), in

Left matrix in the tensor product. It is daggered during the subroutine.

r1INTEGER, in

First dimension of

c1INTEGER, in

Second dimension of

r2INTEGER, in

First dimension of

c2INTEGER, in

Second dimension of

opCHARACTER, in

Either ‘N’ for usual Kronecker product, or ‘+’ for adding to the result.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_id_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity, so either or .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

matinreal(*, *), in

Matrix in the tensor product with the identity

d1INTEGER, in

First dimension of matin.

d2INTEGER, in

Second dimension of matin.

didINTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

idCHARACTER, in

Defines the position of the identity matrix. Either ‘L’ for or ‘R’ for .

transCHARACTER, in

Defines the transformation executed on matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_id_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity, so either or .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

matincomplex(*, *), in

Matrix in the tensor product with the identity

d1INTEGER, in

First dimension of matin.

d2INTEGER, in

Second dimension of matin.

didINTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

idCHARACTER, in

Defines the position of the identity matrix. Either ‘L’ for or ‘R’ for .

transCHARACTER, in

Defines the transformation executed on matrix when taking the tensor product. Either ‘N’ (none), ‘T’ (transposed), or ‘D’ (dagger / conjugate transposed).

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mi_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat1real(*, *), in

Matrix in the tensor product with the identity

r1INTEGER, in

First dimension of the matrix.

c1INTEGER, in

Second dimension of the matrix.

d2INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mi_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat1complex(*, *), in

Matrix in the tensor product with the identity

r1INTEGER, in

First dimension of the matrix.

c1INTEGER, in

Second dimension of the matrix.

d2INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_im_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat2real(*, *), in

Matrix in the tensor product with the identity

d1INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

r2INTEGER, in

First dimension of the matrix.

c2INTEGER, in

Second dimension of the matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_im_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat2complex(*, *), in

Matrix in the tensor product with the identity

d1INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

r2INTEGER, in

First dimension of the matrix.

c2INTEGER, in

Second dimension of the matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mti_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat1real(*, *), in

Matrix in the tensor product with the identity. Transposed while taking the tensor product.

r1INTEGER, in

First dimension of the matrix.

c1INTEGER, in

Second dimension of the matrix.

d2INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mti_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat1complex(*, *), in

Matrix in the tensor product with the identity. Transposed while taking the tensor product.

r1INTEGER, in

First dimension of the matrix.

c1INTEGER, in

Second dimension of the matrix.

d2INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_imt_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat2real(*, *), in

Matrix in the tensor product with the identity. Transposed while taking tensor product.

d1INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

r2INTEGER, in

First dimension of the matrix.

c2INTEGER, in

Second dimension of the matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_imt_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat2complex(*, *), in

Matrix in the tensor product with the identity. Transposed while taking tensor product.

d1INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

r2INTEGER, in

First dimension of the matrix.

c2INTEGER, in

Second dimension of the matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdi_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat1real(*, *), in

Matrix in the tensor product with the identity. Daggered while taking the tensor product.

r1INTEGER, in

First dimension of the matrix.

c1INTEGER, in

Second dimension of the matrix.

d2INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_mdi_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat1complex(*, *), in

Matrix in the tensor product with the identity. Daggered while taking the tensor product.

r1INTEGER, in

First dimension of the matrix.

c1INTEGER, in

Second dimension of the matrix.

d2INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_imd_real()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matreal(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat2real(*, *), in

Matrix in the tensor product with the identity. Daggered while taking tensor product.

d1INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

r2INTEGER, in

First dimension of the matrix.

c2INTEGER, in

Second dimension of the matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.kron_imd_complex()

fortran-subroutine - May 2016 (dj) Kronecker or tensor product between a matrix and the identity .

Arguments

matcomplex(*, *), in

On exit the tensor product between the two matrices. Has to be allocated on entry.

mat2complex(*, *), in

Matrix in the tensor product with the identity. Daggered while taking tensor product.

d1INTEGER, in

Specifies the dimension of the identity matrix. Alwyas square matrix.

r2INTEGER, in

First dimension of the matrix.

c2INTEGER, in

Second dimension of the matrix.

Details

(template defined in LinearAlgebra_include.f90)

Source Code

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

LinearAlgebra_f90.qr_inplace_real()

fortran-subroutine - September 2014 (dj) Compute thin QR-decomposition inplace for nxm-matrix with n >=m

Arguments

matREAL(ld_mat, *), inout

Matrix to be deomposed in QR with nr >= nc. On exit matrix Q.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatREAL(nc, *), out

On exit triangular matrix R of QR-decomposition.

Details

(template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.qr_trapez_real()

fortran-subroutine - April 2016 (dj) Compute thin QR-decomposition inplace for nxm-matrix with n < m

Arguments

matREAL(ld_mat, *), inout

Matrix to be deomposed in QR with nr < nc. On exit entries are destroyed.

qmatREAL(nr, nr), out

On exit, unitary matrix Q of QR-decomposition.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatREAL(nc, *), out

On exit triangular matrix R of QR-decomposition.

Details

(template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.qr_inplace_complex()

fortran-subroutine - September 2014 (dj) Compute thin QR-decomposition inplace for nxm-matrix with n >=m

Arguments

matCOMPLEX(ld_mat, *), inout

Matrix to be deomposed in QR with nr >= nc. On exit matrix Q.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatCOMPLEX(nc, *), out

On exit triangular matrix R of QR-decomposition.

Details

(template defined in LinearAlgebar_template.f90)

Source Code

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

LinearAlgebra_f90.qr_trapez_complex()

fortran-subroutine - April 2016 (dj) Compute thin QR-decomposition inplace for nxm-matrix with n >=m

Arguments

matCOMPLEX(ld_mat, *), inout

Matrix to be deomposed in QR with nr < nc. On exit entries are destroyed.

qmatCOMPLEX(nr, nr), out

On exit, unitary matrix Q of QR decomposition.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatCOMPLEX(nc, *), out

On exit triangular matrix R of QR-decomposition.

Details

(template defined in LinearAlgebar_template.f90)

Source Code

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

LinearAlgebra_f90.rq_inplace_real()

fortran-subroutine - September 2014 (dj) Compute thin RQ-decomposition inplace for nxm-matrix with n <= m

Arguments

matREAL(ld_mat, *), inout

Matrix to be deomposed in RQ with nr <= nc. On exit matrix Q.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatREAL(nc, *), out

On exit triangular matrix R of RQ-decomposition.

Details

(template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.rq_trapez_real()

fortran-subroutine - April 2016 (dj) Compute thin RQ-decomposition inplace for nxm-matrix with n > m

Arguments

matREAL(ld_mat, *), inout

Matrix to be deomposed in RQ with nr > nc. On exit entries are destroyed.

qmatREAL(nc, nc), out

On exit, qmat is unitary matrix Q or RQ-decomposition.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatREAL(nc, *), out

On exit triangular matrix R of RQ-decomposition.

Details

(template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.rq_inplace_complex()

fortran-subroutine - September 2014 (dj) Compute thin RQ-decomposition inplace for nxm-matrix with n <= m

Arguments

matCOMPLEX(ld_mat, *), inout

Matrix to be deomposed in RQ with nr <= nc. On exit matrix Q.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatCOMPLEX(nc, *), out

On exit triangular matrix R of RQ-decomposition.

Details

(template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.rq_trapez_complex()

fortran-subroutine - April 2016 (dj) RQ-decomposition for the case that complex matrix has more rows than colums.

Arguments

matCOMPLEX(ld_mat, *), inout

Matrix to be deomposed in RQ with nr > nc. On exit entries destroyed.

qmatCOMPLEX(nc, nc), inout

On exit, unitary matrix Q of RQ decomposition.

nrINTEGER, in

number of rows in mat

ncINTEGER, in

number of columns in mat

ld_matINTEGER, in

leading dimension of the array mat

rmatCOMPLEX(nc, *), out

On exit triangular matrix R of RQ-decomposition.

Details

(template defined in LinearAlgebra_template.f90)

Source Code

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

LinearAlgebra_f90.svd_elem_real()

fortran-subroutine - September 2014 (dj) Perform a SVD without truncation on real arrays with DGESVD

Arguments

uleftREAL(*, *), out

on exit the left unitary matrix

singREAL(*), out

contains on exit the singular values

vrightREAL(*,*), out

on exit the right unitary matrix

matinREAL(*,*), inout

contains the matrix to be decomposed in a SVD. On exit the entries are destroyed.

a_rowINTEGER, in

number of rows in matrix matin

a_colINTEGER, in

number of columns in matin

infoINTEGER, inout

exit status of ZGESDD (0: sucessful, <0 : wrong argument, >0 : not converged)

Source Code

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

LinearAlgebra_f90.svdd_elem_real()

fortran-subroutine - September 2014 (dj) Perform an SVD without truncation on real arrays with DGESDD.

Arguments

uREAL(*, *), out

on exit the left unitary matrix

sREAL(*), out

contains on exit the singular values

vREAL(*,*), out

on exit the right unitary matrix

aREAL(*,*), inout

contains the matrix to be decomposed in a SVD. On exit the entries are destroyed.

a_rowINTEGER, in

number of rows in matrix a

a_colINTEGER, in

number of columns in a

jobzCHARACTER, 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

infoINTEGER, inout

exit status of ZGESDD (0: sucessful, <0: wrong argument, >0: not converged)

Source Code

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

LinearAlgebra_f90.svd_elem_complex()

fortran-subroutine - September 2014 (dj) Perform a SVD without truncation on complex arrays with ZGESVD

Arguments

uleftCOMPLEX(*, *), out

on exit the left unitary matrix

singREAL(*), out

contains on exit the singular values

vrightCOMPLEX(*,*), out

on exit the right unitary matrix

matinCOMPLEX(*,*), inout

contains the matrix to be decomposed in a SVD. On exit the entries are destroyed.

a_rowINTEGER, in

number of rows in matrix matin

a_colINTEGER, in

number of columns in matin

infoINTEGER, inout

exit status of ZGESDD (0: sucessful, <0: wrong argument, >0: not converged)

Source Code

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

LinearAlgebra_f90.svdd_elem_complex()

fortran-subroutine - September 2014 (dj) Perform a SVD without truncation on complex arrays with ZGESDD

Arguments

uCOMPLEX(*, *), out

on exit the left unitary matrix

sREAL(*), out

contains on exit the singular values

vCOMPLEX(*,*), out

on exit the right unitary matrix

aCOMPLEX(*,*), inout

contains the matrix to be decomposed in a SVD. On exit the entries are destroyed.

a_rowINTEGER, in

number of rows in matrix a

a_colINTEGER, in

number of columns in a

jobzCHARACTER, 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

infoINTEGER, inout

exit status of ZGESDD (0: sucessful, <0: wrong argument, >0: not converged)

Source Code

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

LinearAlgebra_f90.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

rhoRHO_TYOE(*, *), in

Density matrix. Must be hermitian since this property is used.

matreal(*, *), in

Operator matrix.

dimINTEGER, in

Dimension of the two square matrices.

Source Code

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

LinearAlgebra_f90.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

rhoRHO_TYOE(*, *), in

Density matrix. Must be hermitian since this property is used.

matreal(*, *), in

Operator matrix.

dimINTEGER, in

Dimension of the two square matrices.

Source Code

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

LinearAlgebra_f90.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

rhoRHO_TYOE(*, *), in

Density matrix. Must be hermitian since this property is used.

matcomplex(*, *), in

Operator matrix.

dimINTEGER, in

Dimension of the two square matrices.

Source Code

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

LinearAlgebra_f90.tracematmul_real()

fortran-function - November 2017 (dj) Calculate the trace of a matrix multiplication without contructing the intermediate matrix.

Arguments

matlreal(*, *), in

Left matrix in the multiplication.

matrreal(*, *), in

Right matrix in the multuplication.

dimINTEGER, in

Dimension of the square matrices (right now we restrict the function to square matrices).

translCHARACTER, in

Transformation of the left matrix. ‘N’ for no transformation, ‘T’ for transposed, and ‘C’ for the hermitian conjugated.

transrCHARACTER, in

Transformation of the right matrix. ‘N’ for no transformation, ‘T’ for transposed, and ‘C’ for the hermitian conjugated.

Source Code

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

LinearAlgebra_f90.tracematmul_complex()

fortran-function - November 2017 (dj) Calculate the trace of a matrix multiplication without contructing the intermediate matrix.

Arguments

matlcomplex(*, *), in

Left matrix in the multiplication.

matrcomplex(*, *), in

Right matrix in the multuplication.

dimINTEGER, in

Dimension of the square matrices (right now we restrict the function to square matrices).

translCHARACTER, in

Transformation of the left matrix. ‘N’ for no transformation, ‘T’ for transposed, and ‘C’ for the hermitian conjugated.

transrCHARACTER, in

Transformation of the right matrix. ‘N’ for no transformation, ‘T’ for transposed, and ‘C’ for the hermitian conjugated.

Source Code

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