Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/test_kron.f90
3324689e437cfeature/refactoring

Last Change
Raw File

test_kron.f90
View Options

!>
!> @file test_kron.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/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
! Test Kronecker product for both GE and CSR versions.
!
USE iso_fortran_env, ONLY : rkind => real64
USE matrix, ONLY : gemat, init, kron, vmx
USE csr, ONLY : csr_mat, init, full_to_csr, kron, vmx
IMPLICIT NONE
!
TYPE(gemat) :: mata, matb, matc
TYPE(gemat) :: matu1
TYPE(csr_mat) :: mata_csr, matb_csr, matc_csr
INTEGER, PARAMETER :: m=3, n=2
REAL(rkind) :: a(m*m)=[ 1.1, 0.0, 0.0, 0.2, 1.5, 1.0, 0.5, 0.0, 1.0 ]
REAL(rkind) :: b(n*n)=[ 2.0, 1.0, 0.0, 3.0 ]
REAL(rkind), TARGET :: u(m,n), uu(m,n), v(m,n)
REAL(rkind), POINTER :: u1d(:), uu1d(:)
INTEGER :: i, s, e
!
CALL init(m, 0, mata)
CALL init(n, 0, matb)
mata%val(1:m,1:m) = RESHAPE(a, [m,m])
matb%val(1:n,1:n) = RESHAPE(b, [n,n])
!
CALL printmat_ge('Matrix A', mata)
CALL printmat_ge('Matrix B', matb)
!
u1d(1:m*n) => u ! u1d = vec(u)
u1d = [ (REAL(i,rkind), i=1,m*n) ]
CALL printmat('Array U', u)
!
! Compute (A.U).B^T
!
CALL init(n, 0, matu1, mrows=m)
v = TRANSPOSE(vmx(matb, TRANSPOSE(vmx(mata, u))))
CALL printmat('(A.U).B^T', v)
!
! Compute (BxA).vec(U)
!
CALL kron(matb, mata, matc)
uu1d(1:m*n) => uu
uu1d = vmx(matc, u1d)
CALL printmat_ge('Matrix C', matc)
CALL printmat('(BxA).vec(U)', uu)
!----------------------------------------------------------------------
!
! Using CSR matrices
!
CALL full_to_csr(mata%val, mata_csr)
CALL full_to_csr(matb%val, matb_csr)
CALL printmat_csr('Matrix A', mata_csr)
CALL printmat_csr('Matrix B', matb_csr)
!
CALL kron(matb_csr, mata_csr, matc_csr)
uu1d = vmx(matc_csr, u1d)
!
CALL printmat_csr('Matrix C', matc_csr)
CALL printmat('(BxA).vec(U)', uu)
!
CONTAINS
SUBROUTINE printmat(str, a)
CHARACTER(len=*) :: str
REAL(rkind), INTENT(in) :: a(:,:)
INTEGER :: i,m,n
WRITE(*,'(/a)') TRIM(str)
m=SIZE(a,1)
n=SIZE(a,2)
DO i=1,m
WRITE(*,'(12f8.3)') a(i,:)
END DO
END SUBROUTINE printmat
SUBROUTINE printmat_ge(str, a)
CHARACTER(len=*) :: str
TYPE(gemat) :: a
INTEGER :: i
WRITE(*,'(/a)') TRIM(str)
DO i=1,a%mrows
WRITE(*,'(12f8.3)') a%val(i,:)
END DO
END SUBROUTINE printmat_ge
SUBROUTINE printmat_csr(str, a)
CHARACTER(len=*) :: str
TYPE(csr_mat) :: a
INTEGER :: i, s, e
REAL(rkind) :: arow(a%ncols)
WRITE(*,'(/a,a,3i4)') TRIM(str), ': m, n, nnz', a%mrows, a%ncols, a%nnz
DO i=1,a%mrows
arow = 0.0_rkind
s = a%irow(i)
e = a%irow(i+1)-1
arow(a%cols(s:e)) = a%val(s:e)
WRITE(*,'(12f8.3)') arow
END DO
IF(SIZE(a%idiag) .GT. 0) THEN
WRITE(*,'(a,(20i4))') 'idiag =', a%idiag
WRITE(*,'(a,12f8.3)') 'diag = ', a%val(a%idiag)
END IF
!
END SUBROUTINE printmat_csr
END PROGRAM main
c4science · Help