Homec4science
Diffusion SPClibs (master) Browse

spclibs/examples/tmatrix_pb.f90
03540012fedemaster

Last Change
Raw File

tmatrix_pb.f90
View Options

!>
!> @file tmatrix_pb.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>
!>
PROGRAM main
!
! Test some routines of module matrix
!
USE matrix
IMPLICIT NONE
TYPE(pbmat) :: mata, matb
INTEGER, PARAMETER :: n=5, ku=3
DOUBLE PRECISION :: arr(n), fulla(n,n), fullb(n,n), base
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: p, pb
INTEGER :: i, j, info, pow
!
CALL init(ku, n, 0, mata)
CALL init(1, n, 0, matb)
CALL getvalp(mata, p)
CALL getvalp(matb, pb)
PRINT*, 'shape of A: ', SHAPE(p)
!
! Test updtmat
p = 0.0d0
DO i=1,n
DO j=i,n
arr(j) = 10*i + j
IF( ABS(j-i) .LT. ku+1 ) CALL updtmat(mata, i, j, arr(j))
END DO
END DO
CALL prntmat('Test of UPDTMAT', p)
!
! Test GETCOL
fulla = 0.0
DO j=1,n
CALL getcol(mata, j, fulla(:,j))
END DO
CALL prntmat('Full matrix from GETCOL', fulla)
!
! Test GETROW
fulla = 0.0
DO i=1,n
CALL getrow(mata, i, fulla(i,:))
END DO
CALL prntmat('Full matrix from GETROW', fulla)
!
! Test PUTCOL
p = 0.0d0
DO j=1,n
DO i=1,n
arr(i) = 10*i + j
END DO
CALL putcol(mata, j, arr)
END DO
CALL prntmat('Test of PUTCOL', p)
!
! Test PUTROW
p = 0.0d0
DO i=1,n
DO j=1,n
arr(j) = 10*i + j
END DO
CALL putrow(mata, i, arr)
END DO
CALL prntmat('Test of PUTROW', p)
!
arr = 0.0d0
arr(2) =1.0
CALL putrow(mata, 2, arr)
CALL prntmat('Away on i=2, j=2', p)
!
! Test GETCOL
fulla = 0.0
DO j=1,n
CALL getcol(mata, j, fulla(:,j))
END DO
CALL prntmat('Full matrix from GETCOL', fulla)
!
! Test GETROW
fulla = 0.0
DO i=1,n
CALL getrow(mata, i, fulla(i,:))
END DO
CALL prntmat('Full matrix from GETROW', fulla)
!
! Test GETELE
fulla = 0.0
DO i=1,n
DO j=1,n
IF(ABS(j-i).LT.ku+1) CALL getele(mata,i,j,fulla(i,j))
END DO
END DO
CALL prntmat('Full matrix from GETELE', fulla)
!
! Test of determinant
fullb = 0.0
DO i=1,n
fullb(i,i) = 2.0d0
IF(i.LT.n) fullb(i,i+1)=-1.0d0
IF(i.GT.1) fullb(i,i-1)=-1.0d0
END DO
DO j=1,n
CALL putcol(matb, j, fullb(:,j))
END DO
CALL prntmat('Mat. A (full)', fullb)
CALL prntmat('Mat. A (pb)', pb)
CALL determinant(matb, base, pow)
WRITE(*,'(a,1pe15.6,i6)') 'Determinant(bas,power) =', base, pow
!
CALL destroy(mata)
CALL destroy(matb)
CONTAINS
SUBROUTINE prntmat(str, a)
DOUBLE PRECISION, DIMENSION(:,:) :: a
CHARACTER(len=*) :: str
INTEGER :: i
WRITE(*,'(a)') TRIM(str)
DO i=1,SIZE(a,1)
WRITE(*,'(10f8.1)') a(i,:)
END DO
END SUBROUTINE prntmat
END PROGRAM main
c4science · Help