Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/matrix.f90
93a1e9712d04feature/refactoring

Last Change
Raw File

matrix.f90
View Options

!>
!> @file matrix.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE matrix
!
! MATRIX: Simple interface to the direct solver LAPACK.
!
! T.M. Tran, CRPP-EPFL
! February 2007
!
IMPLICIT NONE
!
TYPE gbmat ! Lapack General Band matrix storage
INTEGER :: kl, ku, rank
INTEGER :: mrows, ncols
INTEGER :: nterms, kmat
INTEGER, DIMENSION(:), POINTER :: piv => null()
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: val => null()
END TYPE gbmat
!
TYPE gemat ! Lapack General DENSE matrix storage
INTEGER :: rank
INTEGER :: mrows, ncols
INTEGER :: nterms, kmat
INTEGER, DIMENSION(:), POINTER :: piv => null()
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: val => null()
END TYPE gemat
!
TYPE pbmat ! Lapack Pack Band matrix storage (super-diagonals)
INTEGER :: ku, rank
INTEGER :: nterms, kmat
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: val => null()
END TYPE pbmat
!
TYPE zgbmat ! Lapack General Band matrix storage
INTEGER :: kl, ku, rank
INTEGER :: mrows, ncols
INTEGER :: nterms, kmat
INTEGER, DIMENSION(:), POINTER :: piv => null()
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: val => NULL()
END TYPE zgbmat
!
TYPE zgemat ! Lapack General DENSE matrix storage
INTEGER :: rank
INTEGER :: mrows, ncols
INTEGER :: nterms, kmat
INTEGER, DIMENSION(:), POINTER :: piv => null()
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: val => null()
END TYPE zgemat
!
TYPE zpbmat ! Lapack Pack Band matrix storage (super-diagonals)
INTEGER :: ku, rank
INTEGER :: nterms, kmat
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: val => null()
END TYPE zpbmat
!
TYPE periodic_mat
TYPE(gbmat) :: mat
INTEGER :: nterms
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: &
& matu => null(), &
& matvt => null()
END TYPE periodic_mat
!
TYPE zperiodic_mat
TYPE(zgbmat) :: mat
INTEGER :: nterms
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: &
& matu => null(), &
& matvt => null()
END TYPE zperiodic_mat
!
!--------------------------------------------------------------------------------
INTERFACE init
MODULE PROCEDURE init_gb, init_ge, init_pb, &
& init_zgb, init_zge, init_zpb, &
& init_periodic, init_zperiodic
END INTERFACE
INTERFACE getvalp
MODULE PROCEDURE getvalp_gb, getvalp_ge, getvalp_pb, &
& getvalp_zgb, getvalp_zge, getvalp_zpb
END INTERFACE
INTERFACE mcopy
MODULE PROCEDURE mcopy_gb, mcopy_ge, mcopy_pb, &
& mcopy_zgb, mcopy_zge, mcopy_zpb, &
& mcopy_periodic, mcopy_zperiodic
END INTERFACE
INTERFACE maddto
MODULE PROCEDURE maddto_gb, maddto_ge, maddto_pb, &
& maddto_zgb, maddto_zge, maddto_zpb, &
& maddto_periodic, maddto_zperiodic
END INTERFACE
INTERFACE destroy
MODULE PROCEDURE destroy_gb, destroy_ge, destroy_pb, &
& destroy_zgb, destroy_zge, destroy_zpb, &
& destroy_periodic, destroy_zperiodic
END INTERFACE
INTERFACE updtmat
MODULE PROCEDURE updt_gb, updt_ge, updt_pb, &
& updt_zgb, updt_zpb, &
& updt_periodic, updt_zperiodic
END INTERFACE
INTERFACE getele
MODULE PROCEDURE getele_gb, getele_pb, &
& getele_zgb, getele_zpb, &
& getele_periodic, getele_zperiodic
END INTERFACE
INTERFACE putele
MODULE PROCEDURE putele_gb, putele_pb, &
& putele_zgb, putele_zpb, &
& putele_periodic, putele_zperiodic
END INTERFACE
INTERFACE getcol
MODULE PROCEDURE getcol_gb, getcol_pb, &
& getcol_zgb, getcol_zpb, &
& getcol_periodic, getcol_zperiodic
END INTERFACE
INTERFACE getrow
MODULE PROCEDURE getrow_gb, getrow_ge, getrow_pb, &
& getrow_zgb, getrow_zpb, &
& getrow_periodic, getrow_zperiodic
END INTERFACE
INTERFACE putcol
MODULE PROCEDURE putcol_gb, putcol_ge, putcol_pb, &
& putcol_zgb, putcol_zpb, &
& putcol_periodic, putcol_zperiodic
END INTERFACE
INTERFACE putrow
MODULE PROCEDURE putrow_gb, putrow_ge, putrow_pb, &
& putrow_zgb, putrow_zpb, &
& putrow_periodic, putrow_zperiodic
END INTERFACE
INTERFACE factor
MODULE PROCEDURE factor_gb, factor_ge, factor_pb, &
& factor_zgb, factor_zge, factor_zpb, &
& factor_periodic, factor_zperiodic
END INTERFACE
INTERFACE bsolve
MODULE PROCEDURE bsolve_gb1, bsolve_gbn, bsolve_ge1, bsolve_gen, &
& bsolve_pb1, bsolve_pbn, &
& bsolve_periodic1, bsolve_periodicn, &
& bsolve_zperiodic1, bsolve_zperiodicn, &
& bsolve_zgb1, bsolve_zgbn, bsolve_zge1, bsolve_zgen, &
& bsolve_zpb1, bsolve_zpbn
END INTERFACE
INTERFACE vmx
MODULE PROCEDURE vmx_gb, vmx_gbn, vmx_pb, vmx_pbn, &
& vmx_zgb, vmx_zgbn, vmx_zpb, vmx_zpbn, &
& vmx_ge, vmx_gen, vmx_zge, vmx_zgen, &
& vmx_periodic, vmx_zperiodic
END INTERFACE
INTERFACE determinant
MODULE PROCEDURE determinant_ge, determinant_gb, determinant_pb, &
& determinant_zge, determinant_zgb, determinant_zpb
END INTERFACE
INTERFACE putmat
MODULE PROCEDURE putmat_gb
END INTERFACE
INTERFACE getmat
MODULE PROCEDURE getmat_gb
END INTERFACE
INTERFACE kron
MODULE PROCEDURE kron_ge
END INTERFACE kron
!
CONTAINS
!===========================================================================
SUBROUTINE init_ge(n, nterms, mat, kmat, mrows)
!
! Initialize Lapack General Dense matrice
!
INTEGER, INTENT(in) :: n, nterms
INTEGER, OPTIONAL :: kmat
INTEGER, OPTIONAL :: mrows
TYPE(gemat) :: mat
!
mat%ncols = n
mat%mrows = n
IF(PRESENT(mrows)) THEN
mat%mrows = mrows
END IF
mat%rank = n ! Warning: ok if square matrix
mat%nterms = nterms
IF(PRESENT(kmat)) mat%kmat = kmat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
ALLOCATE(mat%val(mat%mrows,mat%ncols))
ALLOCATE(mat%piv(MIN(mat%mrows,mat%ncols)))
mat%val = 0.0d0
mat%piv = 0
END SUBROUTINE init_ge
!===========================================================================
SUBROUTINE init_gb(kl, ku, n, nterms, mat, kmat, mrows)
!
! Initialize Lapack General Banded matrice
!
INTEGER, INTENT(in) :: kl, ku, n, nterms
INTEGER, OPTIONAL :: kmat
INTEGER, OPTIONAL :: mrows
TYPE(gbmat) :: mat
INTEGER :: lda
!
mat%kl = kl
mat%ku = ku
mat%ncols = n
mat%mrows = n
IF(PRESENT(mrows)) THEN
mat%mrows = mrows
END IF
mat%rank = n ! Warning: ok if square matrix
mat%nterms = nterms
IF(PRESENT(kmat)) mat%kmat = kmat
lda = 2*kl + ku + 1
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
ALLOCATE(mat%val(lda,n))
ALLOCATE(mat%piv(n))
mat%val = 0.0d0
mat%piv = 0
END SUBROUTINE init_gb
!===========================================================================
SUBROUTINE init_periodic(kl, ku, n, nterms, mat, kmat)
!
! Initialize Lapack Periodic General Banded matrice
!
INTEGER, INTENT(in) :: kl, ku, n, nterms
INTEGER, OPTIONAL :: kmat
TYPE(periodic_mat) :: mat
INTEGER :: i,j
DOUBLE PRECISION, PARAMETER :: zero=0.0d0, one=1.0d0
!
! In band matrix matp%mat is a GB matrix
IF( PRESENT(kmat)) THEN
CALL init(kl, ku, n, nterms, mat%mat, kmat)
ELSE
CALL init(kl, ku, n, nterms, mat%mat)
END IF
mat%nterms = nterms
!
! Off band matrices
IF( ASSOCIATED(mat%matu)) DEALLOCATE(mat%matu)
IF( ASSOCIATED(mat%matvt)) DEALLOCATE(mat%matvt)
ALLOCATE(mat%matu(n, kl+ku))
ALLOCATE(mat%matvt(kl+ku,n))
!
mat%matu = zero
mat%matvt = zero ! kl=3, ku=2
DO j=1,kl ! [ 1 0 0 . . ]
mat%matu(j,j) = one ! [ 0 1 0 . . ]
END DO ! [ 0 0 1 . . ]
! ! [ 0 . . . . ]
DO j=1,ku ! [ . . . . . ]
i=n-ku+j ! [ . . . 1 0 ]
mat%matu(i,ku+j) = one ! [ . . . 0 1 ]
END DO
END SUBROUTINE init_periodic
!===========================================================================
SUBROUTINE init_zperiodic(kl, ku, n, nterms, mat, kmat)
!
! Initialize Lapack Periodic General Banded matrice
!
INTEGER, INTENT(in) :: kl, ku, n, nterms
INTEGER, OPTIONAL :: kmat
TYPE(zperiodic_mat) :: mat
INTEGER :: i,j
DOUBLE PRECISION, PARAMETER :: zero=0.0d0, one=1.0d0
!
! In band matrix matp%mat is a GB matrix
IF( PRESENT(kmat)) THEN
CALL init(kl, ku, n, nterms, mat%mat, kmat)
ELSE
CALL init(kl, ku, n, nterms, mat%mat)
END IF
mat%nterms = nterms
!
! Off band matrices
IF( ASSOCIATED(mat%matu)) DEALLOCATE(mat%matu)
IF( ASSOCIATED(mat%matvt)) DEALLOCATE(mat%matvt)
ALLOCATE(mat%matu(n, kl+ku))
ALLOCATE(mat%matvt(kl+ku,n))
!
mat%matu = zero
mat%matvt = zero ! kl=3, ku=2
DO j=1,kl ! [ 1 0 0 . . ]
mat%matu(j,j) = one ! [ 0 1 0 . . ]
END DO ! [ 0 0 1 . . ]
! ! [ 0 . . . . ]
DO j=1,ku ! [ . . . . . ]
i=n-ku+j ! [ . . . 1 0 ]
mat%matu(i,ku+j) = one ! [ . . . 0 1 ]
END DO
END SUBROUTINE init_zperiodic
!===========================================================================
SUBROUTINE init_pb(ku, n, nterms, mat, kmat)
!
! Initialize Lapack Packed Banded matrice
!
INTEGER, INTENT(in) :: ku, n, nterms
INTEGER, OPTIONAL :: kmat
TYPE(pbmat) :: mat
INTEGER :: lda
!
mat%ku = ku
mat%rank = n
mat%nterms = nterms
IF(PRESENT(kmat)) mat%kmat = kmat
lda = ku + 1
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(lda,n))
mat%val = 0.0d0
END SUBROUTINE init_pb
!===========================================================================
SUBROUTINE init_zge(n, nterms, mat, kmat, mrows)
!
! Initialize Lapack General Dense matrice
!
INTEGER, INTENT(in) :: n, nterms
INTEGER, OPTIONAL :: kmat
INTEGER, OPTIONAL :: mrows
TYPE(zgemat) :: mat
!
mat%ncols = n
mat%mrows = n
IF(PRESENT(mrows)) THEN
mat%mrows = mrows
END IF
mat%rank = n
mat%nterms = nterms
IF(PRESENT(kmat)) mat%kmat = kmat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
ALLOCATE(mat%val(mat%mrows,mat%ncols))
ALLOCATE(mat%piv(MIN(mat%mrows,mat%ncols)))
mat%val = 0.0d0
mat%piv = 0
END SUBROUTINE init_zge
!===========================================================================
SUBROUTINE init_zgb(kl, ku, n, nterms, mat, kmat, mrows)
!
! Initialize Lapack General Banded matrice
!
INTEGER, INTENT(in) :: kl, ku, n, nterms
INTEGER, OPTIONAL :: kmat
INTEGER, OPTIONAL :: mrows
TYPE(zgbmat) :: mat
INTEGER :: lda
!
mat%kl = kl
mat%ku = ku
mat%ncols = n
mat%mrows = n
IF(PRESENT(mrows)) THEN
mat%mrows = mrows
END IF
mat%rank = n ! Warning: ok if square matrix
mat%nterms = nterms
IF(PRESENT(kmat)) mat%kmat = kmat
lda = 2*kl + ku + 1
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
ALLOCATE(mat%val(lda,n))
ALLOCATE(mat%piv(n))
mat%val = 0.0d0
mat%piv = 0
END SUBROUTINE init_zgb
!===========================================================================
SUBROUTINE init_zpb(ku, n, nterms, mat, kmat)
!
! Initialize Lapack Packed Banded matrice
!
INTEGER, INTENT(in) :: ku, n, nterms
INTEGER, OPTIONAL :: kmat
TYPE(zpbmat) :: mat
INTEGER :: lda
!
mat%ku = ku
mat%rank = n
mat%nterms = nterms
IF(PRESENT(kmat)) mat%kmat = kmat
lda = ku + 1
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(lda,n))
mat%val = 0.0d0
END SUBROUTINE init_zpb
!===========================================================================
SUBROUTINE mcopy_ge(mata, matb)
!
! Matrix copy: B = A
!
TYPE(gemat) :: mata, matb
!
matb%rank = mata%rank
matb%mrows = mata%mrows
matb%ncols = mata%ncols
matb%nterms = mata%nterms
matb%kmat = mata%kmat
IF( ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
IF( ASSOCIATED(matb%piv)) DEALLOCATE(matb%piv)
ALLOCATE(matb%val(matb%mrows,matb%ncols))
ALLOCATE(matb%piv(MIN(matb%mrows,matb%ncols)))
matb%val = mata%val
matb%piv = mata%piv
END SUBROUTINE mcopy_ge
!===========================================================================
SUBROUTINE mcopy_gb(mata, matb)
!
! Matrix copy: B = A
!
TYPE(gbmat) :: mata, matb
INTEGER :: n, lda
!
n = mata%rank
matb%kl = mata%kl
matb%ku = mata%ku
matb%rank = mata%rank
matb%mrows = mata%mrows
matb%ncols = mata%ncols
matb%nterms = mata%nterms
matb%kmat = mata%kmat
lda = 2*mata%kl + mata%ku + 1
IF( ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
IF( ASSOCIATED(matb%piv)) DEALLOCATE(matb%piv)
ALLOCATE(matb%val(lda,n))
ALLOCATE(matb%piv(n))
matb%val = mata%val
matb%piv = mata%piv
END SUBROUTINE mcopy_gb
!===========================================================================
SUBROUTINE mcopy_periodic(mata, matb)
!
! Matrix copy: B = A
!
TYPE(periodic_mat) :: mata, matb
INTEGER :: n, kl, ku
!
kl = mata%mat%kl
ku = mata%mat%ku
n = mata%mat%rank
!
CALL mcopy(mata%mat, matb%mat)
IF( ASSOCIATED(matb%matu)) DEALLOCATE(matb%matu)
IF( ASSOCIATED(matb%matvt)) DEALLOCATE(matb%matvt)
ALLOCATE(matb%matu(n,kl+ku))
ALLOCATE(matb%matvt(kl+ku,n))
matb%matu = mata%matu
matb%matvt = mata%matvt
END SUBROUTINE mcopy_periodic
!===========================================================================
SUBROUTINE mcopy_zperiodic(mata, matb)
!
! Matrix copy: B = A
!
TYPE(zperiodic_mat) :: mata, matb
INTEGER :: n, kl, ku
!
kl = mata%mat%kl
ku = mata%mat%ku
n = mata%mat%rank
!
CALL mcopy(mata%mat, matb%mat)
IF( ASSOCIATED(matb%matu)) DEALLOCATE(matb%matu)
IF( ASSOCIATED(matb%matvt)) DEALLOCATE(matb%matvt)
ALLOCATE(matb%matu(n,kl+ku))
ALLOCATE(matb%matvt(kl+ku,n))
matb%matu = mata%matu
matb%matvt = mata%matvt
END SUBROUTINE mcopy_zperiodic
!===========================================================================
SUBROUTINE mcopy_pb(mata, matb)
!
! Matrix copy: B = A
!
TYPE(pbmat) :: mata, matb
INTEGER :: n, lda
!
n = mata%rank
matb%ku = mata%ku
matb%rank = mata%rank
matb%nterms = mata%nterms
lda = mata%ku + 1
IF( ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
ALLOCATE(matb%val(lda,n))
matb%val = mata%val
END SUBROUTINE mcopy_pb
!===========================================================================
SUBROUTINE mcopy_zge(mata, matb)
!
! Matrix copy: B = A
!
TYPE(zgemat) :: mata, matb
INTEGER :: n
!
n = mata%rank
matb%rank = n
IF( ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
IF( ASSOCIATED(matb%piv)) DEALLOCATE(matb%piv)
ALLOCATE(matb%val(n,n))
ALLOCATE(matb%piv(n))
matb%val = mata%val
matb%piv = mata%piv
END SUBROUTINE mcopy_zge
!===========================================================================
SUBROUTINE mcopy_zgb(mata, matb)
!
! Matrix copy: B = A
!
TYPE(zgbmat) :: mata, matb
INTEGER :: n, lda
!
n = mata%rank
matb%kl = mata%kl
matb%ku = mata%ku
matb%rank = mata%rank
matb%nterms = mata%nterms
lda = 2*mata%kl + mata%ku + 1
IF( ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
IF( ASSOCIATED(matb%piv)) DEALLOCATE(matb%piv)
ALLOCATE(matb%val(lda,n))
ALLOCATE(matb%piv(n))
matb%val = mata%val
matb%piv = mata%piv
END SUBROUTINE mcopy_zgb
!===========================================================================
SUBROUTINE mcopy_zpb(mata, matb)
!
! Matrix copy: B = A
!
TYPE(zpbmat) :: mata, matb
INTEGER :: n, lda
!
n = mata%rank
matb%ku = mata%ku
matb%rank = mata%rank
matb%nterms = mata%nterms
lda = mata%ku + 1
IF( ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
ALLOCATE(matb%val(lda,n))
matb%val = mata%val
END SUBROUTINE mcopy_zpb
!===========================================================================
SUBROUTINE maddto_ge(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(gemat) :: mata, matb
DOUBLE PRECISION :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_ge
!===========================================================================
SUBROUTINE maddto_gb(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(gbmat) :: mata, matb
DOUBLE PRECISION :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_gb
!===========================================================================
SUBROUTINE maddto_periodic(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(periodic_mat) :: mata, matb
DOUBLE PRECISION :: alpha
!
mata%mat%val = mata%mat%val + alpha*matb%mat%val
mata%matvt = mata%matvt + alpha*matb%matvt
END SUBROUTINE maddto_periodic
!===========================================================================
SUBROUTINE maddto_zperiodic(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(zperiodic_mat) :: mata, matb
DOUBLE COMPLEX :: alpha
!
mata%mat%val = mata%mat%val + alpha*matb%mat%val
mata%matvt = mata%matvt + alpha*matb%matvt
END SUBROUTINE maddto_zperiodic
!===========================================================================
SUBROUTINE maddto_pb(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(pbmat) :: mata, matb
DOUBLE PRECISION :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_pb
!===========================================================================
SUBROUTINE maddto_zge(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(zgemat) :: mata, matb
DOUBLE COMPLEX :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_zge
!===========================================================================
SUBROUTINE maddto_zgb(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(zgbmat) :: mata, matb
DOUBLE COMPLEX :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_zgb
!===========================================================================
SUBROUTINE maddto_zpb(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(zpbmat) :: mata, matb
DOUBLE COMPLEX :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_zpb
!===========================================================================
SUBROUTINE getvalp_ge(mat, p)
!
! Get pointer to matrix coefficients
!
TYPE(gemat) :: mat
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: p
!
p => mat%val
END SUBROUTINE getvalp_ge
!===========================================================================
SUBROUTINE getvalp_gb(mat, p)
!
! Get pointer to matrix coefficients
!
TYPE(gbmat) :: mat
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: p
!
p => mat%val
END SUBROUTINE getvalp_gb
!===========================================================================
SUBROUTINE getvalp_pb(mat, p)
!
! Get pointer to matrix coefficients
!
TYPE(pbmat) :: mat
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: p
!
p => mat%val
END SUBROUTINE getvalp_pb
!===========================================================================
SUBROUTINE getvalp_zge(mat, p)
!
! Get pointer to matrix coefficients
!
TYPE(zgemat) :: mat
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: p
!
p => mat%val
END SUBROUTINE getvalp_zge
!===========================================================================
SUBROUTINE getvalp_zgb(mat, p)
!
! Get pointer to matrix coefficients
!
TYPE(zgbmat) :: mat
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: p
!
p => mat%val
END SUBROUTINE getvalp_zgb
!===========================================================================
SUBROUTINE getvalp_zpb(mat, p)
!
! Get pointer to matrix coefficients
!
TYPE(zpbmat) :: mat
DOUBLE COMPLEX, DIMENSION(:,:), POINTER :: p
!
p => mat%val
END SUBROUTINE getvalp_zpb
!===========================================================================
SUBROUTINE destroy_gb(mat)
!
! Deallocate pointers in mat
!
TYPE(gbmat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
END SUBROUTINE destroy_gb
!===========================================================================
SUBROUTINE destroy_ge(mat)
!
! Deallocate pointers in mat
!
TYPE(gemat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
END SUBROUTINE destroy_ge
!===========================================================================
SUBROUTINE destroy_pb(mat)
!
! Deallocate pointers in mat
!
TYPE(pbmat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
END SUBROUTINE destroy_pb
!===========================================================================
SUBROUTINE destroy_zgb(mat)
!
! Deallocate pointers in mat
!
TYPE(zgbmat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
END SUBROUTINE destroy_zgb
!===========================================================================
SUBROUTINE destroy_zge(mat)
!
! Deallocate pointers in mat
!
TYPE(zgemat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%piv)) DEALLOCATE(mat%piv)
END SUBROUTINE destroy_zge
!===========================================================================
SUBROUTINE destroy_zpb(mat)
!
! Deallocate pointers in mat
!
TYPE(zpbmat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
END SUBROUTINE destroy_zpb
!===========================================================================
SUBROUTINE destroy_periodic(mat)
!
! Deallocate pointers in mat
!
TYPE(periodic_mat) :: mat
IF( ASSOCIATED(mat%matu)) DEALLOCATE(mat%matu)
IF( ASSOCIATED(mat%matvt)) DEALLOCATE(mat%matvt)
CALL destroy(mat%mat)
END SUBROUTINE destroy_periodic
!===========================================================================
SUBROUTINE destroy_zperiodic(mat)
!
! Deallocate pointers in mat
!
TYPE(zperiodic_mat) :: mat
IF( ASSOCIATED(mat%matu)) DEALLOCATE(mat%matu)
IF( ASSOCIATED(mat%matvt)) DEALLOCATE(mat%matvt)
CALL destroy(mat%mat)
END SUBROUTINE destroy_zperiodic
!===========================================================================
SUBROUTINE updt_gb(mat, i, j, val)
!
! Update element Aij into banded matrix
!
TYPE(gbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
ib = mat%kl + mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n) .OR. (j .LT. 1)) THEN
WRITE(*,*) 'UPDT: i, j out of range ', i, j, lda, n
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(ib,j) = mat%val(ib,j) + val
END SUBROUTINE updt_gb
!===========================================================================
SUBROUTINE updt_ge(mat, i, j, val)
!
! Update element Aij into banded matrix
!
TYPE(gemat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
!
IF( (i .GT. mat%mrows) .OR. (j .GT. mat%ncols) .OR. (j .LT. 1) .OR. (i.LT.1)) THEN
WRITE(*,*) 'UPDT: i, j out of range ', i, j, mat%mrows, mat%ncols
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(i,j) = mat%val(i,j) + val
END SUBROUTINE updt_ge
!===========================================================================
SUBROUTINE updt_periodic(mat, i, j, val)
!
! Update element Aij into periodic banded matrix
!
TYPE(periodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: n, kl, ku
!
n = mat%mat%rank
kl = mat%mat%kl
ku = mat%mat%ku
!
IF( i.LE.kl .AND. j.GE.n-kl+1 ) THEN
!
! Put into C(1:kl, 1:kl) = V^T(1:kl, n-kl+1:n)
mat%matvt(i,j) = mat%matvt(i,j) + val
!
ELSE IF( i.GE.n-ku+1 .AND. j.LE.ku ) THEN
!
! Put into D(1:ku,1:ku) = V^T(kl+1:kl+ku, 1:ku)
mat%matvt(i-n+kl+ku,j) = mat%matvt(i-n+kl+ku,j) + val
!
ELSE
!
! Put into the banded matrix
CALL updtmat(mat%mat, i, j, val)
!
END IF
END SUBROUTINE updt_periodic
!===========================================================================
SUBROUTINE updt_zperiodic(mat, i, j, val)
!
! Update element Aij into periodic banded matrix
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
INTEGER :: n, kl, ku
!
n = mat%mat%rank
kl = mat%mat%kl
ku = mat%mat%ku
!
IF( i.LE.kl .AND. j.GE.n-kl+1 ) THEN
!
! Put into C(1:kl, 1:kl) = V^T(1:kl, n-kl+1:n)
mat%matvt(i,j) = mat%matvt(i,j) + val
!
ELSE IF( i.GE.n-ku+1 .AND. j.LE.ku ) THEN
!
! Put into D(1:ku,1:ku) = V^T(kl+1:kl+ku, 1:ku)
mat%matvt(i-n+kl+ku,j) = mat%matvt(i-n+kl+ku,j) + val
!
ELSE
!
! Put into the banded matrix
CALL updtmat(mat%mat, i, j, val)
!
END IF
END SUBROUTINE updt_zperiodic
!===========================================================================
SUBROUTINE updt_pb(mat, i, j, val)
!
! Update element Aij into banded matrix
!
TYPE(pbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
IF( i .LE. j ) THEN
ib = mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n) .OR. (j .LT. 1)) THEN
WRITE(*,*) 'UPDT: i, j out of range ', i, j, lda, n
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(ib,j) = mat%val(ib,j) + val
END IF
END SUBROUTINE updt_pb
!===========================================================================
SUBROUTINE updt_zgb(mat, i, j, val)
!
! Update element Aij into banded matrix
!
TYPE(zgbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
ib = mat%kl + mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n) .OR. (j .LT. 1)) THEN
WRITE(*,*) 'UPDT: i, j out of range ', i, j, lda, n
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(ib,j) = mat%val(ib,j) + val
END SUBROUTINE updt_zgb
!===========================================================================
SUBROUTINE updt_zpb(mat, i, j, val)
!
! Update element Aij into banded matrix
!
TYPE(zpbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
IF( i .LE. j ) THEN
ib = mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n) .OR. (j .LT. 1)) THEN
WRITE(*,*) 'UPDT: i, j out of range ', i, j, lda, n
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(ib,j) = mat%val(ib,j) + val
END IF
END SUBROUTINE updt_zpb
!===========================================================================
SUBROUTINE getele_gb(mat, i, j, val)
!
! Get element (i,j) of matrix
!
TYPE(gbmat), INTENT(in) :: mat
DOUBLE PRECISION, INTENT (OUT) :: val
INTEGER, INTENT (IN) :: i, j
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
ib = mat%kl + mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
val = mat%val(ib,j)
END SUBROUTINE getele_gb
!===========================================================================
SUBROUTINE getele_periodic(mat, i, j, val)
!
! Get element Aij of periodic banded matrix
!
TYPE(periodic_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(out) :: val
INTEGER :: n, kl, ku
!
n = mat%mat%rank
kl = mat%mat%kl
ku = mat%mat%ku
!
IF( i.LE.kl .AND. j.GE.n-kl+1 ) THEN
!
! Put into C(1:kl, 1:kl) = V^T(1:kl, n-kl+1:n)
val =mat%matvt(i,j)
!
ELSE IF( i.GE.n-ku+1 .AND. j.LE.ku ) THEN
!
! Put into D(1:ku,1:ku) = V^T(kl+1:kl+ku, 1:ku)
val = mat%matvt(i-n+kl+ku,j)
!
ELSE
!
! Put into the banded matrix
CALL getele(mat%mat, i, j, val)
!
END IF
END SUBROUTINE getele_periodic
!===========================================================================
SUBROUTINE getele_zperiodic(mat, i, j, val)
!
! Get element Aij of periodic banded matrix
!
TYPE(zperiodic_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(out) :: val
INTEGER :: n, kl, ku
!
n = mat%mat%rank
kl = mat%mat%kl
ku = mat%mat%ku
!
IF( i.LE.kl .AND. j.GE.n-kl+1 ) THEN
!
! Put into C(1:kl, 1:kl) = V^T(1:kl, n-kl+1:n)
val =mat%matvt(i,j)
!
ELSE IF( i.GE.n-ku+1 .AND. j.LE.ku ) THEN
!
! Put into D(1:ku,1:ku) = V^T(kl+1:kl+ku, 1:ku)
val = mat%matvt(i-n+kl+ku,j)
!
ELSE
!
! Put into the banded matrix
CALL getele(mat%mat, i, j, val)
!
END IF
END SUBROUTINE getele_zperiodic
!===========================================================================
SUBROUTINE getele_pb(mat, i, j, val)
!
! Get element (i,j) of matrix
!
TYPE(pbmat), INTENT(in) :: mat
DOUBLE PRECISION, INTENT (OUT) :: val
INTEGER, INTENT (IN) :: i, j
INTEGER :: lda, n, ib, irow, jcol
!
lda = SIZE(mat%val, 1)
n = mat%rank
IF( i .LE. j ) THEN ! Upper triangular matrix
irow = i; jcol = j
ELSE ! Lower triangular matrix
irow = j; jcol = i
END IF
ib = mat%ku + irow - jcol + 1
IF( (ib .GT. lda) .OR. (jcol .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
val = mat%val(ib,jcol)
END SUBROUTINE getele_pb
!===========================================================================
SUBROUTINE getele_zgb(mat, i, j, val)
!
! Get element (i,j) of matrix
!
TYPE(zgbmat), INTENT(in) :: mat
DOUBLE COMPLEX, INTENT (OUT) :: val
INTEGER, INTENT (IN) :: i, j
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
ib = mat%kl + mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
val = mat%val(ib,j)
END SUBROUTINE getele_zgb
!===========================================================================
SUBROUTINE getele_zpb(mat, i, j, val)
!
! Get element (i,j) of matrix
!
TYPE(zpbmat), INTENT(in) :: mat
DOUBLE COMPLEX, INTENT (OUT) :: val
INTEGER, INTENT (IN) :: i, j
INTEGER :: lda, n, ib, irow, jcol
!
lda = SIZE(mat%val, 1)
n = mat%rank
!
IF( i .LE. j ) THEN ! Upper triangular matrix
irow = i; jcol = j
ib = mat%ku + irow - jcol + 1
IF( (ib .GT. lda) .OR. (jcol .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
val = mat%val(ib,jcol)
RETURN
ELSE ! Lower triangular matrix
irow = j; jcol = i
ib = mat%ku + irow - jcol + 1
IF( (ib .GT. lda) .OR. (jcol .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
val = CONJG(mat%val(ib,jcol))
END IF
END SUBROUTINE getele_zpb
!===========================================================================
SUBROUTINE putele_gb(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(gbmat), INTENT(inout) :: mat
DOUBLE PRECISION, INTENT (in) :: val
INTEGER, INTENT (in) :: i, j
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
ib = mat%kl + mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(ib,j) = val
END SUBROUTINE putele_gb
!===========================================================================
SUBROUTINE putele_periodic(mat, i, j, val)
!
! Put element Aij into periodic banded matrix
!
TYPE(periodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: n, kl, ku
!
n = mat%mat%rank
kl = mat%mat%kl
ku = mat%mat%ku
!
IF( i.LE.kl .AND. j.GE.n-kl+1 ) THEN
!
! Put into C(1:kl, 1:kl) = V^T(1:kl, n-kl+1:n)
mat%matvt(i,j) = val
!
ELSE IF( i.GE.n-ku+1 .AND. j.LE.ku ) THEN
!
! Put into D(1:ku,1:ku) = V^T(kl+1:kl+ku, 1:ku)
mat%matvt(i-n+kl+ku,j) = val
!
ELSE
!
! Put into the banded matrix
CALL putele(mat%mat, i, j, val)
!
END IF
END SUBROUTINE putele_periodic
!===========================================================================
SUBROUTINE putele_zperiodic(mat, i, j, val)
!
! Put element Aij into periodic banded matrix
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
INTEGER :: n, kl, ku
!
n = mat%mat%rank
kl = mat%mat%kl
ku = mat%mat%ku
!
IF( i.LE.kl .AND. j.GE.n-kl+1 ) THEN
!
! Put into C(1:kl, 1:kl) = V^T(1:kl, n-kl+1:n)
mat%matvt(i,j) = val
!
ELSE IF( i.GE.n-ku+1 .AND. j.LE.ku ) THEN
!
! Put into D(1:ku,1:ku) = V^T(kl+1:kl+ku, 1:ku)
mat%matvt(i-n+kl+ku,j) = val
!
ELSE
!
! Put into the banded matrix
CALL putele(mat%mat, i, j, val)
!
END IF
END SUBROUTINE putele_zperiodic
!===========================================================================
SUBROUTINE putele_pb(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(pbmat), INTENT(inout) :: mat
DOUBLE PRECISION, INTENT (in) :: val
INTEGER, INTENT (IN) :: i, j
INTEGER :: lda, n, ib, irow, jcol
!
lda = SIZE(mat%val, 1)
n = mat%rank
IF( i .LE. j ) THEN ! Upper triangular matrix
irow = i; jcol = j
ELSE ! Lower triangular matrix
irow = j; jcol = i
END IF
ib = mat%ku + irow - jcol + 1
IF( (ib .GT. lda) .OR. (jcol .GT. n)) THEN
WRITE(*,*) 'PUTELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
mat%val(ib,jcol) = val
END SUBROUTINE putele_pb
!===========================================================================
SUBROUTINE putele_zgb(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(zgbmat), INTENT(inout) :: mat
DOUBLE COMPLEX, INTENT (in) :: val
INTEGER, INTENT (in) :: i, j
INTEGER :: lda, n, ib
!
lda = SIZE(mat%val, 1)
n = mat%rank
ib = mat%kl + mat%ku + i - j + 1
IF( (ib .GT. lda) .OR. (j .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(ib,j) = val
END SUBROUTINE putele_zgb
!===========================================================================
SUBROUTINE putele_zpb(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(zpbmat), INTENT(inout) :: mat
DOUBLE COMPLEX, INTENT (in) :: val
INTEGER, INTENT (IN) :: i, j
INTEGER :: lda, n, ib, irow, jcol
!
lda = SIZE(mat%val, 1)
n = mat%rank
IF( i .LE. j ) THEN ! Upper triangular matrix
irow = i; jcol = j
ib = mat%ku + irow - jcol + 1
IF( (ib .GT. lda) .OR. (jcol .GT. n)) THEN
WRITE(*,*) 'GETELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
mat%val(ib,jcol) = val
ELSE ! Lower triangular matrix
irow = j; jcol = i
ib = mat%ku + irow - jcol + 1
IF( (ib .GT. lda) .OR. (jcol .GT. n)) THEN
WRITE(*,*) 'PUTELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE matrix***'
END IF
mat%val(ib,jcol) = CONJG(val)
END IF
END SUBROUTINE putele_zpb
!===========================================================================
SUBROUTINE getcol_gb(mat, j, arr)
!
! Get a column from matrix
!
TYPE(gbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: m, kl, ku
INTEGER :: ibmin, ibmax, imin, imax
!
kl = mat%kl
ku = mat%ku
m = mat%mrows
IF( SIZE(arr) .LT. m ) THEN
WRITE(*,*) 'GETCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:m) = 0.0d0
imin = MAX(1,j-ku)
imax = MIN(m, j+kl)
ibmin = kl+ku+imin-j+1
ibmax = kl+ku+imax-j+1
arr(imin:imax) = mat%val(ibmin:ibmax,j)
END SUBROUTINE getcol_gb
!===========================================================================
SUBROUTINE getcol_periodic(mat, j, arr)
!
! Get a column from matrix
!
TYPE(periodic_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL getcol(mat%mat, j, arr)
!
IF( j.GE.n-kl+1 ) THEN
arr(1:kl) = mat%matvt(1:kl,j)
ELSE IF( j.LE.ku ) THEN
arr(n-ku+1:n) = mat%matvt(kl+1:kl+ku,j)
END IF
END SUBROUTINE getcol_periodic
!===========================================================================
SUBROUTINE getcol_zperiodic(mat, j, arr)
!
! Get a column from matrix
!
TYPE(zperiodic_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL getcol(mat%mat, j, arr)
!
IF( j.GE.n-kl+1 ) THEN
arr(1:kl) = mat%matvt(1:kl,j)
ELSE IF( j.LE.ku ) THEN
arr(n-ku+1:n) = mat%matvt(kl+1:kl+ku,j)
END IF
END SUBROUTINE getcol_zperiodic
!===========================================================================
SUBROUTINE getcol_pb(mat, j, arr)
!
! Get a column from matrix
!
TYPE(pbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, ku
INTEGER :: i, ib, ibmin, ibmax, imin, imax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
!
imin=MAX(1,j-ku); imax=j ! The column in the upper diagonal part
ibmin=ku+1+imin-j ; ibmax=ku+1+imax-j
arr(imin:imax) = mat%val(ibmin:ibmax,j)
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
arr(i) = mat%val(ib,i)
END DO
END SUBROUTINE getcol_pb
!===========================================================================
SUBROUTINE getcol_zgb(mat, j, arr)
!
! Get a column from matrix
!
TYPE(zgbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: arr
INTEGER :: m, kl, ku
INTEGER :: ibmin, ibmax, imin, imax
!
kl = mat%kl
ku = mat%ku
m = mat%mrows
IF( SIZE(arr) .LT. m ) THEN
WRITE(*,*) 'GETCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:m) = 0.0d0
imin = MAX(1,j-ku)
imax = MIN(m, j+kl)
ibmin = kl+ku+imin-j+1
ibmax = kl+ku+imax-j+1
arr(imin:imax) = mat%val(ibmin:ibmax,j)
END SUBROUTINE getcol_zgb
!===========================================================================
SUBROUTINE getcol_zpb(mat, j, arr)
!
! Get a column from matrix
!
TYPE(zpbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, ku
INTEGER :: i, ib, ibmin, ibmax, imin, imax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
!
imin=MAX(1,j-ku); imax=j ! The column in the upper diagonal part
ibmin=ku+1+imin-j ; ibmax=ku+1+imax-j
arr(imin:imax) = mat%val(ibmin:ibmax,j)
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
arr(i) = CONJG(mat%val(ib,i))
END DO
END SUBROUTINE getcol_zpb
!===========================================================================
SUBROUTINE getrow_gb(mat, i, arr)
!
! Get a row from matrix
!
TYPE(gbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, kl, ku
INTEGER :: j, ib, jmin, jmax
!
kl = mat%kl
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
jmin = MAX(1,i-kl)
jmax = MIN(n, i+ku)
DO j=jmin,jmax
ib = kl+ku+i-j+1
arr(j) = mat%val(ib,j)
END DO
END SUBROUTINE getrow_gb
!===========================================================================
SUBROUTINE getrow_ge(mat, i, arr)
!
! Get a row from matrix
!
TYPE(gemat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n
!
n = mat%ncols
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = mat%val(i,1:n)
END SUBROUTINE getrow_ge
!===========================================================================
SUBROUTINE getrow_periodic(mat, i, arr)
!
! Get a row from matrix
!
TYPE(periodic_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL getrow(mat%mat, i, arr)
!
IF( i.LE.kl ) THEN
arr(n-kl+1:n) = mat%matvt(i,n-kl+1:n)
ELSE IF( i.GE.n-ku+1 ) THEN
arr(1:ku) = mat%matvt(i-n+kl+ku,1:ku)
END IF
END SUBROUTINE getrow_periodic
!===========================================================================
SUBROUTINE getrow_zperiodic(mat, i, arr)
!
! Get a row from matrix
!
TYPE(zperiodic_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL getrow(mat%mat, i, arr)
!
IF( i.LE.kl ) THEN
arr(n-kl+1:n) = mat%matvt(i,n-kl+1:n)
ELSE IF( i.GE.n-ku+1 ) THEN
arr(1:ku) = mat%matvt(i-n+kl+ku,1:ku)
END IF
END SUBROUTINE getrow_zperiodic
!===========================================================================
SUBROUTINE getrow_pb(mat, i, arr)
!
! Get a row from matrix
!
TYPE(pbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, ku
INTEGER :: j, ib, ibmin, ibmax, jmin, jmax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
!
jmin=i; jmax=MIN(n,i+ku)
DO j=jmin,jmax
ib=ku+1+i-j
arr(j) = mat%val(ib,j)
END DO
!
jmin=MAX(1,i-ku); jmax=i-1
DO j=jmin,jmax
ib=ku+1+j-i
arr(j) = mat%val(ib,i)
END DO
END SUBROUTINE getrow_pb
!===========================================================================
SUBROUTINE getrow_zgb(mat, i, arr)
!
! Get a row from matrix
!
TYPE(zgbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, kl, ku
INTEGER :: j, ib, jmin, jmax
!
kl = mat%kl
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
jmin = MAX(1,i-kl)
jmax = MIN(n, i+ku)
DO j=jmin,jmax
ib = kl+ku+i-j+1
arr(j) = mat%val(ib,j)
END DO
END SUBROUTINE getrow_zgb
!===========================================================================
SUBROUTINE getrow_zpb(mat, i, arr)
!
! Get a row from matrix
!
TYPE(zpbmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: arr
INTEGER :: n, ku
INTEGER :: j, ib, ibmin, ibmax, jmin, jmax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
!
jmin=i; jmax=MIN(n,i+ku)
DO j=jmin,jmax
ib=ku+1+i-j
arr(j) = mat%val(ib,j)
END DO
!
jmin=MAX(1,i-ku); jmax=i-1
DO j=jmin,jmax
ib=ku+1+j-i
arr(j) = CONJG(mat%val(ib,i))
END DO
END SUBROUTINE getrow_zpb
!===========================================================================
SUBROUTINE putcol_gb(mat, j, arr)
!
! Put a column from matrix
!
TYPE(gbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: m, kl, ku
INTEGER :: ibmin, ibmax, imin, imax
!
kl = mat%kl
ku = mat%ku
m = mat%mrows
IF( SIZE(arr) .LT. m ) THEN
WRITE(*,*) 'PUTCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
imin = MAX(1,j-ku)
imax = MIN(m, j+kl)
ibmin = kl+ku+imin-j+1
ibmax = kl+ku+imax-j+1
mat%val(ibmin:ibmax,j) = arr(imin:imax)
END SUBROUTINE putcol_gb
!===========================================================================
SUBROUTINE putcol_ge(mat, j, arr)
!
! Put a column from matrix
!
TYPE(gemat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: m
!
m = mat%mrows
IF( SIZE(arr) .LT. m ) THEN
WRITE(*,*) 'PUTCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(:,j) = arr(:)
END SUBROUTINE putcol_ge
!===========================================================================
SUBROUTINE putrow_periodic(mat, i, arr)
!
! Put a row to matrix
!
TYPE(periodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL putrow(mat%mat, i, arr)
!
IF( i.LE.kl ) THEN
mat%matvt(i,n-kl+1:n) = arr(n-kl+1:n)
ELSE IF( i.GE.n-ku+1 ) THEN
mat%matvt(i-n+kl+ku,1:ku) = arr(1:ku)
END IF
END SUBROUTINE putrow_periodic
!===========================================================================
SUBROUTINE putrow_zperiodic(mat, i, arr)
!
! Put a row to matrix
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL putrow(mat%mat, i, arr)
!
IF( i.LE.kl ) THEN
mat%matvt(i,n-kl+1:n) = arr(n-kl+1:n)
ELSE IF( i.GE.n-ku+1 ) THEN
mat%matvt(i-n+kl+ku,1:ku) = arr(1:ku)
END IF
END SUBROUTINE putrow_zperiodic
!===========================================================================
SUBROUTINE putcol_periodic(mat, j, arr)
!
! Put a column into matrix
!
TYPE(periodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL putcol(mat%mat, j, arr)
!
IF( j.GE.n-kl+1 ) THEN
mat%matvt(1:kl,j) = arr(1:kl)
ELSE IF( j.LE.ku ) THEN
mat%matvt(kl+1:kl+ku,j) = arr(n-ku+1:n)
END IF
END SUBROUTINE putcol_periodic
!===========================================================================
SUBROUTINE putcol_zperiodic(mat, j, arr)
!
! Put a column into matrix
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
CALL putcol(mat%mat, j, arr)
!
IF( j.GE.n-kl+1 ) THEN
mat%matvt(1:kl,j) = arr(1:kl)
ELSE IF( j.LE.ku ) THEN
mat%matvt(kl+1:kl+ku,j) = arr(n-ku+1:n)
END IF
END SUBROUTINE putcol_zperiodic
!===========================================================================
SUBROUTINE putcol_pb(mat, j, arr)
!
! Put a column from matrix
!
TYPE(pbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, ku
INTEGER :: i, ib, ibmin, ibmax, imin, imax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
imin = MAX(1,j-ku); imax = j ! The column in the upper diagonal part
ibmin = ku+imin-j+1; ibmax = ku+imax-j+1
mat%val(ibmin:ibmax,j) = arr(imin:imax)
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
mat%val(ib,i) = arr(i)
END DO
END SUBROUTINE putcol_pb
!===========================================================================
SUBROUTINE putcol_zgb(mat, j, arr)
!
! Put a column from matrix
!
TYPE(zgbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: arr
INTEGER :: m, kl, ku
INTEGER :: ibmin, ibmax, imin, imax
!
kl = mat%kl
ku = mat%ku
m = mat%mrows
IF( SIZE(arr) .LT. m ) THEN
WRITE(*,*) 'PUTCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
imin = MAX(1,j-ku)
imax = MIN(m, j+kl)
ibmin = kl+ku+imin-j+1
ibmax = kl+ku+imax-j+1
mat%val(ibmin:ibmax,j) = arr(imin:imax)
END SUBROUTINE putcol_zgb
!===========================================================================
SUBROUTINE putcol_zpb(mat, j, arr)
!
! Put a column from matrix
!
TYPE(zpbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, ku, i, ib
INTEGER :: ibmin, ibmax, imin, imax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
!
imin=MAX(1,j-ku); imax=j ! The column in the upper diagonal part
ibmin=ku+1+imin-j ; ibmax=ku+1+imax-j
mat%val(ibmin:ibmax,j) = arr(imin:imax)
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
mat%val(ib,i) = CONJG(arr(i))
END DO
END SUBROUTINE putcol_zpb
!===========================================================================
SUBROUTINE putrow_gb(mat, i, arr)
!
! Put a row from matrix
!
TYPE(gbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
INTEGER :: j, ib, jmin, jmax
!
kl = mat%kl
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
jmin = MAX(1,i-kl)
jmax = MIN(n, i+ku)
DO j=jmin,jmax
ib = kl+ku+i-j+1
mat%val(ib,j) = arr(j)
END DO
END SUBROUTINE putrow_gb
!===========================================================================
SUBROUTINE putrow_ge(mat, i, arr)
!
! Put a row from matrix
!
TYPE(gemat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
INTEGER :: j, ib, jmin, jmax
!
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
mat%val(i,:) = arr(:)
END SUBROUTINE putrow_ge
!===========================================================================
SUBROUTINE putrow_pb(mat, i, arr)
!
! Put a row from matrix
!
TYPE(pbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, ku
INTEGER :: j, ib, jmin, jmax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
jmin = i
jmax = MIN(n, i+ku)
DO j=jmin,jmax
ib = ku+i-j+1
mat%val(ib,j) = arr(j)
END DO
!
jmin=MAX(1,i-ku); jmax=i-1
DO j=jmin,jmax
ib=ku+1+j-i
mat%val(ib,i) = arr(j)
END DO
END SUBROUTINE putrow_pb
!===========================================================================
SUBROUTINE putrow_zgb(mat, i, arr)
!
! Put a row from matrix
!
TYPE(zgbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, kl, ku
INTEGER :: j, ib, jmin, jmax
!
kl = mat%kl
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
jmin = MAX(1,i-kl)
jmax = MIN(n, i+ku)
DO j=jmin,jmax
ib = kl+ku+i-j+1
mat%val(ib,j) = arr(j)
END DO
END SUBROUTINE putrow_zgb
!===========================================================================
SUBROUTINE putrow_zpb(mat, i, arr)
!
! Put a row from matrix
!
TYPE(zpbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: arr
INTEGER :: n, ku
INTEGER :: j, ib, jmin, jmax
!
ku = mat%ku
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
jmin = i
jmax = MIN(n, i+ku)
DO j=jmin,jmax
ib = ku+i-j+1
mat%val(ib,j) = arr(j)
END DO
!
jmin=MAX(1,i-ku); jmax=i-1
DO j=jmin,jmax
ib=ku+1+j-i
mat%val(ib,i) = CONJG(arr(j))
END DO
END SUBROUTINE putrow_zpb
!===========================================================================
SUBROUTINE factor_gb(mat,flops)
!
! Factor the matrix, using Lapack
!
TYPE(gbmat), INTENT(inout) :: mat
DOUBLE PRECISION, OPTIONAL, INTENT (OUT) :: flops
INTEGER :: lda, n, m, kl, ku
INTEGER :: info
DOUBLE PRECISION :: dopgb
EXTERNAL dopgb
!
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
m = mat%mrows
n = mat%ncols
CALL dgbtrf(m, n, kl, ku, mat%val, lda, mat%piv, info)
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRF ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
IF( PRESENT(flops) ) THEN
flops = dopgb('DGBTRF',m, n, kl, ku, mat%piv)
END IF
END SUBROUTINE factor_gb
!===========================================================================
SUBROUTINE factor_periodic(mat)
!
! Factor the periodic GB matrix, using the
! Sherman-Morrisson-Woodburry formula
!
TYPE(periodic_mat), INTENT(inout) :: mat
TYPE(gemat) :: hmat
DOUBLE PRECISION :: one=1.0d0
INTEGER :: bandw, mr, nc, i
!
bandw = SIZE(mat%matvt,1)
!
! Factor A
CALL factor(mat%mat)
!
IF(bandw .EQ. 0 ) RETURN ! No off band terms
!
! U <-- A^(-1) * U
CALL bsolve(mat%mat, mat%matu)
!
! H <-- 1 + V^T * U
mr = SIZE(mat%matvt, 1)
nc = SIZE(mat%matvt, 2)
CALL init(mr, 0, hmat) ! hmat is initialized to 0!
DO i=1,mr
hmat%val(i,i) = one
END DO
CALL dgemm('N', 'N', mr, mr, nc, one, mat%matvt, mr, &
& mat%matu, nc, one, hmat%val, mr)
!
! V^T <-- H^(-1) V^T
CALL factor(hmat)
CALL bsolve(hmat, mat%matvt)
CALL destroy(hmat)
!
END SUBROUTINE factor_periodic
!===========================================================================
SUBROUTINE factor_zperiodic(mat)
!
! Factor the periodic GB matrix, using the
! Sherman-Morrisson-Woodburry formula
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
TYPE(zgemat) :: hmat
DOUBLE COMPLEX :: one=1.0d0
INTEGER :: bandw, mr, nc, i
!
bandw = SIZE(mat%matvt,1)
!
! Factor A
CALL factor(mat%mat)
!
IF(bandw .EQ. 0 ) RETURN ! No off band terms
!
! U <-- A^(-1) * U
CALL bsolve(mat%mat, mat%matu)
!
! H <-- 1 + V^T * U
mr = SIZE(mat%matvt, 1)
nc = SIZE(mat%matvt, 2)
CALL init(mr, 0, hmat) ! hmat is initialized to 0!
DO i=1,mr
hmat%val(i,i) = one
END DO
CALL zgemm('N', 'N', mr, mr, nc, one, mat%matvt, mr, &
& mat%matu, nc, one, hmat%val, mr)
!
! V^T <-- H^(-1) V^T
CALL factor(hmat)
CALL bsolve(hmat, mat%matvt)
CALL destroy(hmat)
!
END SUBROUTINE factor_zperiodic
!===========================================================================
SUBROUTINE factor_pb(mat,flops)
!
! Factor the matrix, using Lapack
!
TYPE(pbmat), INTENT(inout) :: mat
DOUBLE PRECISION, OPTIONAL, INTENT (OUT) :: flops
INTEGER :: lda, n, ku
INTEGER :: info
DOUBLE PRECISION :: dopla
EXTERNAL dopla
!
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
CALL dpbtrf('U', n, ku, mat%val, lda, info)
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from PBTRF ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
IF( PRESENT(flops) ) THEN
flops = dopla('DPBTRF', n, n, ku, ku, 1)
END IF
END SUBROUTINE factor_pb
!===========================================================================
SUBROUTINE factor_ge(mat,flops)
!
! Factor the matrix, using Lapack
!
TYPE(gemat), INTENT(inout) :: mat
DOUBLE PRECISION, OPTIONAL, INTENT (OUT) :: flops
INTEGER :: n, m
INTEGER :: info
DOUBLE PRECISION :: dopla
EXTERNAL dopla
!
m = mat%mrows
n = mat%ncols
CALL dgetrf(m, n, mat%val, m, mat%piv, info)
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GETRF ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
IF( PRESENT(flops) ) THEN
flops = dopla('DGETRF',m, n, 0, 0, 0)
END IF
END SUBROUTINE factor_ge
!===========================================================================
SUBROUTINE factor_zgb(mat,flops)
!
! Factor the matrix, using Lapack
!
TYPE(zgbmat), INTENT(inout) :: mat
DOUBLE PRECISION, OPTIONAL, INTENT (OUT) :: flops
INTEGER :: lda, n, m, kl, ku
INTEGER :: info
DOUBLE PRECISION :: dopgb
EXTERNAL dopgb
!
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
m = mat%mrows
n = mat%ncols
CALL zgbtrf(m, n, kl, ku, mat%val, lda, mat%piv, info)
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRF ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
IF( PRESENT(flops) ) THEN
flops = dopgb('ZGBTRF',m, n, kl, ku, mat%piv)
END IF
END SUBROUTINE factor_zgb
!===========================================================================
SUBROUTINE factor_zpb(mat,flops)
!
! Factor the matrix, using Lapack
!
TYPE(zpbmat), INTENT(inout) :: mat
DOUBLE PRECISION, OPTIONAL, INTENT (OUT) :: flops
INTEGER :: lda, n, ku
INTEGER :: info
DOUBLE PRECISION :: dopla
EXTERNAL dopla
!
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
CALL zpbtrf('U', n, ku, mat%val, lda, info)
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from PBTRF ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
IF( PRESENT(flops) ) THEN
flops = dopla('ZPBTRF', n, n, ku, ku, 1)
END IF
END SUBROUTINE factor_zpb
!===========================================================================
SUBROUTINE factor_zge(mat,flops)
!
! Factor the matrix, using Lapack
!
TYPE(zgemat), INTENT(inout) :: mat
DOUBLE PRECISION, OPTIONAL, INTENT (OUT) :: flops
INTEGER :: n, m
INTEGER :: info
DOUBLE PRECISION :: dopla
EXTERNAL dopla
!
m = mat%mrows
n = mat%ncols
CALL zgetrf(m, n, mat%val, m, mat%piv, info)
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GETRF ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
IF( PRESENT(flops) ) THEN
flops = dopla('ZGETRF',m, n, 0, 0, 0)
END IF
END SUBROUTINE factor_zge
!===========================================================================
SUBROUTINE bsolve_gb1(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(gbmat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:) :: rhs
DOUBLE PRECISION, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, kl, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
n = mat%rank
!
IF( PRESENT(sol) ) THEN
sol = rhs
CALL dgbtrs('N', n, kl, ku, 1, mat%val, lda, mat%piv, sol, n, info)
ELSE
CALL dgbtrs('N', n, kl, ku, 1, mat%val, lda, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_gb1
!===========================================================================
SUBROUTINE bsolve_periodic1(mat, rhs, sol)
!
! Backsolve, using the Sherman-Morrison-Woodburry formula
!
TYPE(periodic_mat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:) :: rhs
DOUBLE PRECISION, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
DOUBLE PRECISION :: one=1.0d0, zero=0.0d0, minus1=-1.0d0
DOUBLE PRECISION, ALLOCATABLE :: tarr(:,:)
INTEGER :: rank, bandw, nrhs
INTEGER :: info
!----------------------------------------------------------------------
rank = mat%mat%rank
bandw = SIZE(mat%matvt,1)
nrhs = 1
!
! Solve Ay = f
IF( PRESENT(sol) ) THEN
CALL bsolve(mat%mat, rhs, sol)
ELSE
CALL bsolve(mat%mat, rhs)
END IF
!
IF(bandw .EQ. 0 ) RETURN ! No off band terms
!
! t = V^T*y ( = W^T*y )
ALLOCATE(tarr(bandw,nrhs))
IF( PRESENT(sol) ) THEN
CALL dgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, sol, &
& rank, zero, tarr, bandw)
ELSE
CALL dgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, rhs, &
& rank, zero, tarr, bandw)
END IF
!
! y = y - U*t ( = y-Z*t)
IF( PRESENT(sol) ) THEN
CALL dgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, sol, rank)
ELSE
CALL dgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, rhs, rank)
END IF
!
DEALLOCATE(tarr)
END SUBROUTINE bsolve_periodic1
!===========================================================================
SUBROUTINE bsolve_zperiodic1(mat, rhs, sol)
!
! Backsolve, using the Sherman-Morrison-Woodburry formula
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:) :: rhs
DOUBLE COMPLEX, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
DOUBLE COMPLEX :: one=1.0d0, zero=0.0d0, minus1=-1.0d0
DOUBLE COMPLEX, ALLOCATABLE :: tarr(:,:)
INTEGER :: rank, bandw, nrhs
INTEGER :: info
!----------------------------------------------------------------------
rank = mat%mat%rank
bandw = SIZE(mat%matvt,1)
nrhs = 1
!
! Solve Ay = f
IF( PRESENT(sol) ) THEN
CALL bsolve(mat%mat, rhs, sol)
ELSE
CALL bsolve(mat%mat, rhs)
END IF
!
IF(bandw .EQ. 0 ) RETURN ! No off band terms
!
! t = V^T*y ( = W^T*y )
ALLOCATE(tarr(bandw,nrhs))
IF( PRESENT(sol) ) THEN
CALL zgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, sol, &
& rank, zero, tarr, bandw)
ELSE
CALL zgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, rhs, &
& rank, zero, tarr, bandw)
END IF
!
! y = y - U*t ( = y-Z*t)
IF( PRESENT(sol) ) THEN
CALL zgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, sol, rank)
ELSE
CALL zgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, rhs, rank)
END IF
!
DEALLOCATE(tarr)
END SUBROUTINE bsolve_zperiodic1
!===========================================================================
SUBROUTINE bsolve_pb1(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(pbmat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:) :: rhs
DOUBLE PRECISION, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
!
IF( PRESENT(sol) ) THEN
sol = rhs
CALL dpbtrs('U', n, ku, 1, mat%val, lda, sol, n, info)
ELSE
CALL dpbtrs('U', n, ku, 1, mat%val, lda, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from PBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_pb1
!===========================================================================
SUBROUTINE bsolve_ge1(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(gemat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:) :: rhs
DOUBLE PRECISION, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
INTEGER :: n
INTEGER :: info
!----------------------------------------------------------------------
n = mat%rank
!
IF( PRESENT(sol) ) THEN
sol = rhs
CALL dgetrs('N', n, 1, mat%val, n, mat%piv, sol, n, info)
ELSE
CALL dgetrs('N', n, 1, mat%val, n, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GETRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_ge1
!===========================================================================
SUBROUTINE bsolve_gbn(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(gbmat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:,:) :: rhs
DOUBLE PRECISION, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, nrhs, kl, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
n = mat%rank
nrhs = SIZE(rhs,2)
!
IF( PRESENT(sol) ) THEN
sol(:,1:nrhs) = rhs(:,1:nrhs)
CALL dgbtrs('N', n, kl, ku, nrhs, mat%val, lda, mat%piv, sol, n, info)
ELSE
CALL dgbtrs('N', n, kl, ku, nrhs, mat%val, lda, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_gbn
!===========================================================================
SUBROUTINE bsolve_periodicn(mat, rhs, sol)
!
! Backsolve, using the Sherman-Morrison-Woodburry formula
!
TYPE(periodic_mat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:,:) :: rhs
DOUBLE PRECISION, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
DOUBLE PRECISION :: one=1.0d0, zero=0.0d0, minus1=-1.0d0
DOUBLE PRECISION, ALLOCATABLE :: tarr(:,:)
INTEGER :: rank, bandw, nrhs
INTEGER :: info
!----------------------------------------------------------------------
rank = mat%mat%rank
bandw = SIZE(mat%matvt,1)
nrhs = SIZE(rhs,2)
!
! Solve Ay = f
IF( PRESENT(sol) ) THEN
CALL bsolve(mat%mat, rhs, sol)
ELSE
CALL bsolve(mat%mat, rhs)
END IF
!
IF(bandw .EQ. 0 ) RETURN ! No off band terms
!
! t = V^T*y ( = W^T*y )
ALLOCATE(tarr(bandw,nrhs))
IF( PRESENT(sol) ) THEN
CALL dgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, sol, &
& rank, zero, tarr, bandw)
ELSE
CALL dgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, rhs, &
& rank, zero, tarr, bandw)
END IF
!
! y = y - U*t ( = y-Z*t)
IF( PRESENT(sol) ) THEN
CALL dgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, sol, rank)
ELSE
CALL dgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, rhs, rank)
END IF
!
DEALLOCATE(tarr)
END SUBROUTINE bsolve_periodicn
!===========================================================================
SUBROUTINE bsolve_zperiodicn(mat, rhs, sol)
!
! Backsolve, using the Sherman-Morrison-Woodburry formula
!
TYPE(zperiodic_mat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:,:) :: rhs
DOUBLE COMPLEX, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
DOUBLE COMPLEX :: one=1.0d0, zero=0.0d0, minus1=-1.0d0
DOUBLE COMPLEX, ALLOCATABLE :: tarr(:,:)
INTEGER :: rank, bandw, nrhs
INTEGER :: info
!----------------------------------------------------------------------
rank = mat%mat%rank
bandw = SIZE(mat%matvt,1)
nrhs = SIZE(rhs,2)
!
! Solve Ay = f
IF( PRESENT(sol) ) THEN
CALL bsolve(mat%mat, rhs, sol)
ELSE
CALL bsolve(mat%mat, rhs)
END IF
!
IF(bandw .EQ. 0 ) RETURN ! No off band terms
!
! t = V^T*y ( = W^T*y )
ALLOCATE(tarr(bandw,nrhs))
IF( PRESENT(sol) ) THEN
CALL zgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, sol, &
& rank, zero, tarr, bandw)
ELSE
CALL zgemm('N', 'N', bandw, nrhs, rank, one, mat%matvt, bandw, rhs, &
& rank, zero, tarr, bandw)
END IF
!
! y = y - U*t ( = y-Z*t)
IF( PRESENT(sol) ) THEN
CALL zgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, sol, rank)
ELSE
CALL zgemm('N', 'N', rank, nrhs, bandw, minus1, mat%matu, rank, tarr, &
& bandw, one, rhs, rank)
END IF
!
DEALLOCATE(tarr)
END SUBROUTINE bsolve_zperiodicn
!===========================================================================
SUBROUTINE bsolve_pbn(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(pbmat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:,:) :: rhs
DOUBLE PRECISION, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, nrhs, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
nrhs = SIZE(rhs,2)
!
IF( PRESENT(sol) ) THEN
sol(:,1:nrhs) = rhs(:,1:nrhs)
CALL dpbtrs('U', n, ku, nrhs, mat%val, lda, sol, n, info)
ELSE
CALL dpbtrs('U', n, ku, nrhs, mat%val, lda, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_pbn
!===========================================================================
SUBROUTINE bsolve_gen(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(gemat), INTENT(inout) :: mat
DOUBLE PRECISION, DIMENSION (:,:) :: rhs
DOUBLE PRECISION, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
INTEGER :: n, nrhs
INTEGER :: info
!----------------------------------------------------------------------
n = mat%rank
nrhs = SIZE(rhs,2)
!
IF( PRESENT(sol) ) THEN
sol(:,1:nrhs) = rhs
CALL dgetrs('N', n, nrhs, mat%val, n, mat%piv, sol, n, info)
ELSE
CALL dgetrs('N', n, nrhs, mat%val, n, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GETRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_gen
!===========================================================================
SUBROUTINE bsolve_zgb1(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(zgbmat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:) :: rhs
DOUBLE COMPLEX, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, kl, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
n = mat%rank
!
IF( PRESENT(sol) ) THEN
sol = rhs
CALL zgbtrs('N', n, kl, ku, 1, mat%val, lda, mat%piv, sol, n, info)
ELSE
CALL zgbtrs('N', n, kl, ku, 1, mat%val, lda, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_zgb1
!===========================================================================
SUBROUTINE bsolve_zpb1(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(zpbmat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:) :: rhs
DOUBLE COMPLEX, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
!
IF( PRESENT(sol) ) THEN
sol = rhs
CALL zpbtrs('U', n, ku, 1, mat%val, lda, sol, n, info)
ELSE
CALL zpbtrs('U', n, ku, 1, mat%val, lda, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from PBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_zpb1
!===========================================================================
SUBROUTINE bsolve_zge1(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(zgemat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:) :: rhs
DOUBLE COMPLEX, DIMENSION (:), OPTIONAL, INTENT (out) :: sol
INTEGER :: n
INTEGER :: info
!----------------------------------------------------------------------
n = mat%rank
!
IF( PRESENT(sol) ) THEN
sol = rhs
CALL zgetrs('N', n, 1, mat%val, n, mat%piv, sol, n, info)
ELSE
CALL zgetrs('N', n, 1, mat%val, n, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GETRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_zge1
!===========================================================================
SUBROUTINE bsolve_zgbn(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(zgbmat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:,:) :: rhs
DOUBLE COMPLEX, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, nrhs, kl, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
n = mat%rank
nrhs = SIZE(rhs,2)
!
IF( PRESENT(sol) ) THEN
sol(:,1:nrhs) = rhs(:,1:nrhs)
CALL zgbtrs('N', n, kl, ku, nrhs, mat%val, lda, mat%piv, sol, n, info)
ELSE
CALL zgbtrs('N', n, kl, ku, nrhs, mat%val, lda, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_zgbn
!===========================================================================
SUBROUTINE bsolve_zpbn(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(zpbmat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:,:) :: rhs
DOUBLE COMPLEX, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
INTEGER :: lda, n, nrhs, ku
INTEGER :: info
!----------------------------------------------------------------------
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
nrhs = SIZE(rhs,2)
!
IF( PRESENT(sol) ) THEN
sol(:,1:nrhs) = rhs(:,1:nrhs)
CALL zpbtrs('U', n, ku, nrhs, mat%val, lda, sol, n, info)
ELSE
CALL zpbtrs('U', n, ku, nrhs, mat%val, lda, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GBTRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_zpbn
!===========================================================================
SUBROUTINE bsolve_zgen(mat, rhs, sol)
!
! Backsolve, using Lapack
!
TYPE(zgemat), INTENT(inout) :: mat
DOUBLE COMPLEX, DIMENSION (:,:) :: rhs
DOUBLE COMPLEX, DIMENSION (:,:), OPTIONAL, INTENT (out) :: sol
INTEGER :: n, nrhs
INTEGER :: info
!----------------------------------------------------------------------
n = mat%rank
nrhs = SIZE(rhs,2)
!
IF( PRESENT(sol) ) THEN
sol(:,1:nrhs) = rhs
CALL zgetrs('N', n, nrhs, mat%val, n, mat%piv, sol, n, info)
ELSE
CALL zgetrs('N', n, nrhs, mat%val, n, mat%piv, rhs, n, info)
END IF
IF( info .NE. 0) THEN
WRITE(*,*) 'FACTOR: info from GETRS ', info
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
END SUBROUTINE bsolve_zgen
!===========================================================================
FUNCTION vmx_gb(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(gbmat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE PRECISION, ALLOCATABLE :: vmx(:)
DOUBLE PRECISION :: one=1.0d0, zero=0.0d0
INTEGER :: lda, kl, ku, m, n, j, imin, imax, ibmin, ibmax
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
m = mat%mrows
n = mat%ncols
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
ALLOCATE(vmx(m))
ELSE
ALLOCATE(vmx(n))
END IF
!
CALL dgbmv(trans_loc, m, n, kl, ku, one, mat%val(kl+1,1), lda, x, 1, zero,&
& vmx, 1)
END FUNCTION vmx_gb
!===========================================================================
FUNCTION vmx_ge(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(gemat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE PRECISION, ALLOCATABLE :: vmx(:)
DOUBLE PRECISION :: one=1.0d0, zero=0.0d0
INTEGER :: lda, m, n
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
m = mat%mrows
n = mat%ncols
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
ALLOCATE(vmx(m))
ELSE
ALLOCATE(vmx(n))
END IF
!
CALL dgemv(trans_loc, m, n, one, mat%val, lda, x, 1, zero, vmx, 1)
END FUNCTION vmx_ge
!===========================================================================
FUNCTION vmx_gen(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(gemat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE PRECISION, ALLOCATABLE :: vmx(:,:)
DOUBLE PRECISION :: one=1.0d0, zero=0.0d0
INTEGER :: lda, ldb, m, n, k
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
ldb = SIZE(x,1)
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
m = mat%mrows
n = SIZE(x,2)
k = mat%ncols
ALLOCATE(vmx(m,n))
ELSE
m = mat%ncols
n = SIZE(x,2)
k = mat%mrows
ALLOCATE(vmx(m,n))
END IF
!
CALL dgemm(trans_loc, 'N', m, n, k, one, mat%val, lda, x, ldb, zero, vmx, &
& lda)
!
END FUNCTION vmx_gen
!===========================================================================
FUNCTION vmx_zge(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(zgemat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE COMPLEX, ALLOCATABLE :: vmx(:)
DOUBLE COMPLEX :: one=1.0d0, zero=0.0d0
INTEGER :: lda, m, n
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
m = mat%mrows
n = mat%ncols
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
ALLOCATE(vmx(m))
ELSE
ALLOCATE(vmx(n))
END IF
!
CALL zgemv(trans_loc, m, n, one, mat%val, lda, x, 1, zero, vmx, 1)
END FUNCTION vmx_zge
!===========================================================================
FUNCTION vmx_zgen(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(zgemat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:,:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE COMPLEX, ALLOCATABLE :: vmx(:,:)
DOUBLE COMPLEX :: one=1.0d0, zero=0.0d0
INTEGER :: lda, ldb, m, n, k
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
ldb = SIZE(x,1)
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
m = mat%mrows
n = SIZE(x,2)
k = mat%ncols
ALLOCATE(vmx(m,n))
ELSE
m = mat%ncols
n = SIZE(x,2)
k = mat%mrows
ALLOCATE(vmx(m,n))
END IF
!
CALL zgemm(trans_loc, 'N', m, n, k, one, mat%val, lda, x, ldb, zero, vmx, &
& lda)
!
END FUNCTION vmx_zgen
!===========================================================================
FUNCTION vmx_periodic(mat, x)
!
! Return product mat*x
!
TYPE(periodic_mat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: x
DOUBLE PRECISION, DIMENSION(SIZE(x)) :: vmx_periodic
INTEGER :: kl, ku, n, i, ii
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
vmx_periodic = vmx(mat%mat, x)
!
DO i=1,kl
vmx_periodic(i) = vmx_periodic(i) + &
& DOT_PRODUCT(mat%matvt(i,n-kl+1:n), x(n-kl+1:n))
END DO
!
DO i=n-ku+1,n
ii = i-n+ku+kl
vmx_periodic(i) = vmx_periodic(i) + &
& DOT_PRODUCT(mat%matvt(ii,1:ku), x(1:ku))
END DO
END FUNCTION vmx_periodic
!===========================================================================
FUNCTION vmx_zperiodic(mat, x)
!
! Return product mat*x
!
TYPE(zperiodic_mat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: x
DOUBLE COMPLEX, DIMENSION(SIZE(x)) :: vmx_zperiodic
INTEGER :: kl, ku, n, i, ii
!
kl = mat%mat%kl
ku = mat%mat%ku
n = mat%mat%rank
!
vmx_zperiodic = vmx(mat%mat, x)
!
DO i=1,kl
vmx_zperiodic(i) = vmx_zperiodic(i) + &
& DOT_PRODUCT(mat%matvt(i,n-kl+1:n), x(n-kl+1:n))
END DO
!
DO i=n-ku+1,n
ii = i-n+ku+kl
vmx_zperiodic(i) = vmx_zperiodic(i) + &
& DOT_PRODUCT(mat%matvt(ii,1:ku), x(1:ku))
END DO
END FUNCTION vmx_zperiodic
!===========================================================================
FUNCTION vmx_gbn(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(gbmat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE PRECISION, ALLOCATABLE :: vmx(:,:)
DOUBLE PRECISION :: one=1.0d0, zero=0.0d0
INTEGER :: lda, kl, ku, m, n, j, k, imin, imax, ibmin, ibmax
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
m = mat%mrows
n = mat%ncols
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
ALLOCATE(vmx(m,SIZE(x,2)))
ELSE
ALLOCATE(vmx(n,SIZE(x,2)))
END IF
!
DO k=1,SIZE(x,2)
CALL dgbmv(trans_loc, m, n, kl, ku, one, mat%val(kl+1,1), lda, &
& x(1,k), 1, zero, vmx(1,k), 1)
END DO
END FUNCTION vmx_gbn
!===========================================================================
FUNCTION vmx_pb(mat, x)
!
! Return product mat*x
!
TYPE(pbmat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: x
DOUBLE PRECISION, DIMENSION(SIZE(x)) :: vmx_pb
INTEGER :: lda, ku, n, i, j, imin, imax, ib, ibmin, ibmax
!
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
!
vmx_pb = 0.0d0
DO j=1,n
imin = MAX(1,j-ku); imax = j ! The column in the upper diagonal part
ibmin = ku+imin-j+1; ibmax = ku+imax-j+1
vmx_pb(imin:imax) = vmx_pb(imin:imax) + mat%val(ibmin:ibmax,j)*x(j)
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
vmx_pb(i) = vmx_pb(i) + mat%val(ib,i)*x(j)
END DO
END DO
END FUNCTION vmx_pb
!===========================================================================
FUNCTION vmx_pbn(mat, x)
!
! Return product mat*x
!
TYPE(pbmat), INTENT(in) :: mat
DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: x
DOUBLE PRECISION, DIMENSION(SIZE(x,1),SIZE(x,2)) :: vmx_pbn
INTEGER :: lda, ku, n, i, j, k, imin, imax, ib, ibmin, ibmax
!
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
!
vmx_pbn = 0.0d0
DO j=1,n
imin = MAX(1,j-ku); imax = j ! The column in the upper diagonal part
ibmin = ku+imin-j+1; ibmax = ku+imax-j+1
DO k=1,SIZE(x,2)
vmx_pbn(imin:imax,k) = vmx_pbn(imin:imax,k) + &
& mat%val(ibmin:ibmax,j)*x(j,k)
END DO
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
vmx_pbn(i,:) = vmx_pbn(i,:) + mat%val(ib,i)*x(j,:)
END DO
END DO
END FUNCTION vmx_pbn
!===========================================================================
FUNCTION vmx_zgb(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(zgbmat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE COMPLEX, ALLOCATABLE :: vmx(:)
DOUBLE COMPLEX :: one=1.0d0, zero=0.0d0
INTEGER :: lda, kl, ku, m, n, j, imin, imax, ibmin, ibmax
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
m = mat%mrows
n = mat%ncols
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
ALLOCATE(vmx(m))
ELSE
ALLOCATE(vmx(n))
END IF
!
CALL zgbmv(trans_loc, m, n, kl, ku, one, mat%val(kl+1,1), lda, x, 1, zero,&
& vmx, 1)
END FUNCTION vmx_zgb
!===========================================================================
FUNCTION vmx_zgbn(mat, x, trans) RESULT(vmx)
!
! Return product mat*x
!
TYPE(zgbmat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:,:), INTENT(in) :: x
CHARACTER(len=1), OPTIONAL :: trans
DOUBLE COMPLEX, ALLOCATABLE :: vmx(:,:)
DOUBLE COMPLEX :: one=1.0d0, zero=0.0d0
INTEGER :: lda, kl, ku, m, n, j, k, imin, imax, ibmin, ibmax
CHARACTER(len=1) :: trans_loc
!
lda = SIZE(mat%val, 1)
kl = mat%kl
ku = mat%ku
m = mat%mrows
n = mat%ncols
trans_loc = 'N'
IF(PRESENT(trans)) trans_loc = trans
!
IF(trans_loc.EQ.'N') THEN
ALLOCATE(vmx(m,SIZE(x,2)))
ELSE
ALLOCATE(vmx(n,SIZE(x,2)))
END IF
!
DO k=1,SIZE(x,2)
CALL zgbmv(trans_loc, m, n, kl, ku, one, mat%val(kl+1,1), lda, &
& x(1,k), 1, zero, vmx(1,k), 1)
END DO
END FUNCTION vmx_zgbn
!===========================================================================
FUNCTION vmx_zpb(mat, x)
!
! Return product mat*x
!
TYPE(zpbmat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:), INTENT(in) :: x
DOUBLE COMPLEX, DIMENSION(SIZE(x)) :: vmx_zpb
INTEGER :: lda, ku, n, i, j, imin, imax, ib, ibmin, ibmax
!
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
!
vmx_zpb = 0.0d0
DO j=1,n
imin = MAX(1,j-ku); imax = j ! The column in the upper diagonal part
ibmin = ku+imin-j+1; ibmax = ku+imax-j+1
vmx_zpb(imin:imax) = vmx_zpb(imin:imax) + mat%val(ibmin:ibmax,j)*x(j)
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
vmx_zpb(i) = vmx_zpb(i) + CONJG(mat%val(ib,i))*x(j)
END DO
END DO
END FUNCTION vmx_zpb
!===========================================================================
FUNCTION vmx_zpbn(mat, x)
!
! Return product mat*x
!
TYPE(zpbmat), INTENT(in) :: mat
DOUBLE COMPLEX, DIMENSION(:,:), INTENT(in) :: x
DOUBLE COMPLEX, DIMENSION(SIZE(x,1),SIZE(x,2)) :: vmx_zpbn
INTEGER :: lda, ku, n, i, j, k, imin, imax, ib, ibmin, ibmax
!
lda = SIZE(mat%val, 1)
ku = mat%ku
n = mat%rank
!
vmx_zpbn = 0.0d0
DO j=1,n
imin = MAX(1,j-ku); imax = j ! The column in the upper diagonal part
ibmin = ku+imin-j+1; ibmax = ku+imax-j+1
DO k=1,SIZE(x,2)
vmx_zpbn(imin:imax,k) = vmx_zpbn(imin:imax,k) + &
& mat%val(ibmin:ibmax,j)*x(j,k)
END DO
!
imin=j+1; imax = MIN(n,j+ku) ! The column in the lower diagonal part
DO i=imin,imax
ib = ku+1+j-i
vmx_zpbn(i,:) = vmx_zpbn(i,:) + CONJG(mat%val(ib,i))*x(j,:)
END DO
END DO
END FUNCTION vmx_zpbn
!===========================================================================
SUBROUTINE determinant_ge(mat, base, pow)
!
! Return the determinant of mat
!
TYPE(gemat) :: mat
INTEGER :: pow, i
DOUBLE PRECISION :: base
!
CALL factor(mat)
base = 1.0d0
pow = 0
DO i=1,mat%rank
IF( mat%piv(i) .NE. i) base = -base
base = mat%val(i,i)*base
IF( base .EQ. 0.0d0 ) THEN
WRITE(*,*) 'DETERMINANT_GE: matrix is singular'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO
IF( ABS(base) .GE. 1.0d0 ) EXIT
base = 10.0d0*base
pow = pow - 1
END DO
DO
IF( ABS(base) .LT. 10.0d0 ) EXIT
base = base/10.0d0
pow = pow + 1
END DO
END DO
END SUBROUTINE determinant_ge
!===========================================================================
SUBROUTINE determinant_gb(mat, base, pow)
!
! Return the determinant of mat
!
TYPE(gbmat) :: mat
INTEGER :: pow, i, ib
DOUBLE PRECISION :: base
!
CALL factor(mat)
base = 1.0d0
pow = 0
ib=mat%kl + mat%ku + 1
DO i=1,mat%rank
IF( mat%piv(i) .NE. i) base = -base
base = mat%val(ib,i)*base
IF( base .EQ. 0.0d0 ) THEN
WRITE(*,*) 'DETERMINANT_GB: matrix is singular'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO
IF( ABS(base) .GE. 1.0d0 ) EXIT
base = 10.0d0*base
pow = pow - 1
END DO
DO
IF( ABS(base) .LT. 10.0d0 ) EXIT
base = base/10.0d0
pow = pow + 1
END DO
END DO
END SUBROUTINE determinant_gb
!===========================================================================
SUBROUTINE determinant_pb(mat, base, pow)
!
! Return the determinant of mat
!
TYPE(pbmat) :: mat
INTEGER :: pow, i, ib
DOUBLE PRECISION :: base
!
CALL factor(mat)
base = 1.0d0
pow = 0
ib = mat%ku + 1
DO i=1,mat%rank
base = mat%val(ib,i)*base
IF( base .EQ. 0.0d0 ) THEN
WRITE(*,*) 'DETERMINANT_PB: matrix is singular'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO
IF( ABS(base) .GE. 1.0d0 ) EXIT
base = 10.0d0*base
pow = pow - 1
END DO
DO
IF( ABS(base) .LT. 10.0d0 ) EXIT
base = base/10.0d0
pow = pow + 1
END DO
END DO
base=base**2
pow=pow*2
END SUBROUTINE determinant_pb
!===========================================================================
SUBROUTINE determinant_zge(mat, base, pow)
!
! Return the determinant of mat
!
TYPE(zgemat) :: mat
INTEGER :: pow, i
DOUBLE COMPLEX :: base
!
CALL factor(mat)
base = 1.0d0
pow = 0
DO i=1,mat%rank
IF( mat%piv(i) .NE. i) base = -base
base = mat%val(i,i)*base
IF( base .EQ. 0.0d0 ) THEN
WRITE(*,*) 'DETERMINANT_ZGE: matrix is singular'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO
IF( ABS(base) .GE. 1.0d0 ) EXIT
base = 10.0d0*base
pow = pow - 1
END DO
DO
IF( ABS(base) .LT. 10.0d0 ) EXIT
base = base/10.0d0
pow = pow + 1
END DO
END DO
END SUBROUTINE determinant_zge
!===========================================================================
SUBROUTINE determinant_zgb(mat, base, pow)
!
! Return the determinant of mat
!
TYPE(zgbmat) :: mat
INTEGER :: pow, i, ib
DOUBLE COMPLEX :: base
!
CALL factor(mat)
base = 1.0d0
pow = 0
ib=mat%kl + mat%ku + 1
DO i=1,mat%rank
IF( mat%piv(i) .NE. i) base = -base
base = mat%val(ib,i)*base
IF( base .EQ. 0.0d0 ) THEN
WRITE(*,*) 'DETERMINANT_GB: matrix is singular'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO
IF( ABS(base) .GE. 1.0d0 ) EXIT
base = 10.0d0*base
pow = pow - 1
END DO
DO
IF( ABS(base) .LT. 10.0d0 ) EXIT
base = base/10.0d0
pow = pow + 1
END DO
END DO
END SUBROUTINE determinant_zgb
!===========================================================================
SUBROUTINE determinant_zpb(mat, base, pow)
!
! Return the determinant of mat
!
TYPE(zpbmat) :: mat
INTEGER :: pow, i, ib
DOUBLE COMPLEX :: base
!
CALL factor(mat)
base = 1.0d0
pow = 0
ib = mat%ku + 1
DO i=1,mat%rank
base = mat%val(ib,i)*base
IF( base .EQ. 0.0d0 ) THEN
WRITE(*,*) 'DETERMINANT_PB: matrix is singular'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO
IF( ABS(base) .GE. 1.0d0 ) EXIT
base = 10.0d0*base
pow = pow - 1
END DO
DO
IF( ABS(base) .LT. 10.0d0 ) EXIT
base = base/10.0d0
pow = pow + 1
END DO
END DO
base=base**2
pow=pow*2
END SUBROUTINE determinant_zpb
!===========================================================================
SUBROUTINE putmat_gb(fid, label, mat, str)
!
! Write GB matrix in hdf5 file
!
USE futils
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(gbmat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
IF(PRESENT(str)) THEN
CALL putarr(fid, label, mat%val, str)
ELSE
CALL putarr(fid, label, mat%val)
END IF
CALL attach(fid, label, 'KL', mat%kl)
CALL attach(fid, label, 'KU', mat%ku)
CALL attach(fid, label, 'RANK', mat%rank)
END SUBROUTINE putmat_gb
!===========================================================================
SUBROUTINE getmat_gb(fid, label, mat, str)
!
! Read in GB matrix from hdf5 file
!
USE futils
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(gbmat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
CALL getatt(fid, label, 'KL', mat%kl)
CALL getatt(fid, label, 'KU', mat%ku)
CALL getatt(fid, label, 'RANK', mat%rank)
CALL getarr(fid, label, mat%val)
END SUBROUTINE getmat_gb
!===========================================================================
SUBROUTINE kron_ge(mata, matb, matc)
!
! Krocnecker product of 2 dense matrices
!
TYPE(gemat), INTENT(in) :: mata, matb
TYPE(gemat), INTENT(out) :: matc
!
INTEGER :: i1, j1, i3, j3, m1, n1, m2, n2, m3, n3
!
m1 = mata%mrows
n1 = mata%ncols
m2 = matb%mrows
n2 = matb%ncols
m3 = m1*m2
n3 = n1*n2
!
CALL init(n3, 0, matc, mrows=m3)
DO i1=1,m1
i3 = (i1-1)*m2
DO j1=1,n1
j3 = (j1-1)*n2
matc%val(i3+1:i3+m2,j3+1:j3+n2) = mata%val(i1,j1)*matb%val(1:m2,1:n2)
END DO
END DO
END SUBROUTINE kron_ge
!===========================================================================
END MODULE matrix
c4science · Help