Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/tsparse2.f90
03540012fedefeature/refactoring

Last Change
Raw File

tsparse2.f90
View Options

!>
!> @file tsparse2.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
!
! Simple 2D Poisson using 5 points FD
!
USE pardiso_bsplines
IMPLICIT NONE
!
TYPE(pardiso_mat) :: amat, bmat
DOUBLE PRECISION, ALLOCATABLE :: arow(:), rhs(:), sol(:)
DOUBLE PRECISION, ALLOCATABLE :: arown(:,:), rhsn(:,:), soln(:,:)
INTEGER :: nx=5, ny=4, n
INTEGER :: i, j, irow, jcol
!
DOUBLE PRECISION :: mem, seconds
!
n = nx*ny
CALL init(n, 1, amat) ! Non-symmetric matrix
CALL init(n, 1, bmat, nlsym=.TRUE.) ! Symmetric matrix
ALLOCATE(arow(n), arown(n,2))
ALLOCATE(rhs(n), rhsn(n,2))
ALLOCATE(sol(n), soln(n,2))
!
! Construct the FD matrix amat, using sparse rows (linked lists)
!
DO j=1,ny
DO i=1,nx
arow = 0.0d0
irow = numb(i,j)
arow(irow) = 4.0d0
IF(i.GT.1) arow(numb(i-1,j)) = -1.0d0
IF(i.LT.nx) arow(numb(i+1,j)) = -1.0d0
IF(j.GT.1) arow(numb(i,j-1)) = -1.0d0
IF(j.LT.ny) arow(numb(i,j+1)) = -1.0d0
rhs(irow) = SUM(arow)
CALL putrow(amat, irow, arow) ! General matrix
CALL putrow(bmat, irow, arow) ! Symmetric matrix
END DO
END DO
!
! Print the matrices
!
WRITE(*,'(/a)') 'Matrix A'
DO i=1,n
CALL getrow(amat, i, arow)
WRITE(*,'(30f4.0)') arow
END DO
PRINT*, 'nnz from get_count', get_count(amat)
!
WRITE(*,'(/a)') 'Matrix B'
DO i=1,n
CALL getrow(bmat, i, arow)
WRITE(*,'(30f4.0)') arow
END DO
PRINT*, 'nnz from get_count', get_count(bmat)
!
! Factor the matrix using Pardiso
!
CALL factor(amat, nlmetis=.TRUE.)
WRITE(*,'(/a,i5)') 'Number of nonzeros in factors of A = ',amat%p%iparm(18)
WRITE(*,'(a,i5)') 'Number of factorization MFLOPS = ',amat%p%iparm(19)
!
CALL factor(bmat, nlmetis=.TRUE.)
WRITE(*,'(/a,i5)') 'Number of nonzeros in factors of B = ',bmat%p%iparm(18)
WRITE(*,'(a,i5)') 'Number of factorization MFLOPS = ',bmat%p%iparm(19)
!
WRITE(*,'(/a/(10f8.4))') 'rhs', rhs
!
! Backsolve Ax = b, using Pardiso
!
sol = rhs
CALL bsolve(amat, sol, debug=.FALSE.)
WRITE(*,'(/a/(10f8.4))') 'sol (non-sym)', sol
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
!
! Backsolve Bx = b, using Pardiso
!
sol = rhs
CALL bsolve(bmat, sol, debug=.FALSE.)
WRITE(*,'(/a/(10f8.4))') 'sol (sym)', sol
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
!
arow = vmx(amat, sol)
WRITE(*,'(/a/(10f8.4))') 'A*x (non-sym)', arow
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arow-rhs))
!
arow = vmx(bmat, sol)
WRITE(*,'(/a/(10f8.4))') 'B*x (sym)', arow
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arow-rhs))
!
! Multiple RHS
!
rhsn(:,1) = -rhs(:)
rhsn(:,2) = 2.d0*rhs(:)
CALL bsolve(amat, rhsn, soln)
WRITE(*,'(/a/(10f8.4))') 'soln (non-sym)', soln
WRITE(*,'(a,2(1pe12.3))') 'Error', MAXVAL(ABS(soln(:,1)+1.0d0)), &
& MAXVAL(ABS(soln(:,2)-2.0d0))
!
CALL bsolve(bmat, rhsn, soln)
WRITE(*,'(/a/(10f8.4))') 'soln (sym)', soln
WRITE(*,'(a,2(1pe12.3))') 'Error', MAXVAL(ABS(soln(:,1)+1.0d0)), &
& MAXVAL(ABS(soln(:,2)-2.0d0))
!
arown = vmx(amat, soln)
WRITE(*,'(/a/(10f8.4))') 'A*x (non-sym)', arown
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arown-rhsn))
!
arown = vmx(bmat, soln)
WRITE(*,'(/a/(10f8.4))') 'A*x (sym)', arown
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arown-rhsn))
!
! Clean up
!
DEALLOCATE(arow,arown)
DEALLOCATE(rhs,rhsn)
DEALLOCATE(sol,soln)
CALL destroy(amat)
CALL destroy(bmat)
CONTAINS
INTEGER FUNCTION numb(i,j)
INTEGER, INTENT(in) :: i, j
numb = (j-1)*nx + i
END FUNCTION numb
!
END PROGRAM main
c4science · Help