Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/poisson_mumps.f90
f1da25c60953feature/refactoring

Last Change
Raw File

poisson_mumps.f90
View Options

!>
!> @file poisson_mumps.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
USE mumps_bsplines
USE cds
!
IMPLICIT NONE
INCLUDE 'mpif.h'
TYPE(mumps_mat) :: amat
TYPE(cds_mat) :: amat_cds
DOUBLE PRECISION, ALLOCATABLE :: rhs(:), sol(:), arow(:)
INTEGER :: nx=5, ny=4
INTEGER :: n
INTEGER :: i, j, irow
INTEGER :: ierr, me
INTEGER, ALLOCATABLE :: dists(:)
DOUBLE PRECISION :: t0
!
CALL mpi_init(ierr)
CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
!
IF(me.EQ.0) THEN
WRITE(*,'(a)', advance='no') 'Enter nx, ny: '
READ(*,*) nx, ny
END IF
CALL mpi_bcast(nx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(ny, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
n = nx*ny ! Rank of the matrix
ALLOCATE(rhs(n))
ALLOCATE(sol(n))
ALLOCATE(arow(n))
!
WRITE(*,'(/a)') 'Mumps using CSR mat ...'
CALL init(n, 1, amat)
!
! Construct the matrix and RHS
!
t0 = mpi_wtime(0)
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
CALL putrow(amat, irow, arow)
rhs(irow) = SUM(arow) ! => the exact solution is 1
END DO
END DO
!
WRITE(*,'(a,i6)') 'Rank of matrix', n
WRITE(*,'(a,i6)') 'Number of non-zeros of matrix', get_count(amat)
WRITE(*,'(a,1pe12.3)') 'Matrix construction time (s)', mpi_wtime()-t0
!
! Factor the amat matrix (Reordering, symbolic and numerical factorization)
!
t0 = mpi_wtime(0)
CALL factor(amat, nlmetis=.TRUE.)
sol=rhs
CALL bsolve(amat, sol)
WRITE(*,'(a,1pe12.3)') 'Direct solve time (s)', mpi_wtime()-t0
!
IF(me.EQ.0) THEN
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
END IF
CALL destroy(amat)
!
! CDS matrix
!
WRITE(*,'(/a)') 'Mumps using CDS mat ...'
IF(ALLOCATED(dists)) DEALLOCATE(dists)
ALLOCATE(dists(-2:2))
dists = [-nx, -1, 0, 1, nx]
WRITE(*,'(a/(20i4))') 'dists used in INIT =', dists
CALL init(n, dists, 1, amat_cds)
!
t0 = mpi_wtime(0)
DO j=1,ny
DO i=1,nx
arow = 0.0d0
irow = numb(i,j)
amat_cds%val(irow,0) = 4.0d0
IF(i.GT.1) amat_cds%val(irow,-1) = -1.0d0
IF(i.LT.nx) amat_cds%val(irow,+1) = -1.0d0
IF(j.GT.1) amat_cds%val(irow,-2) = -1.0d0
IF(j.LT.ny) amat_cds%val(irow,+2) = -1.0d0
END DO
END DO
WRITE(*,'(a,1pe12.3)') 'Matrix construction time (s)', mpi_wtime()-t0
!
! Compute dists of amat
PRINT*, 'stat of mata%mat', ASSOCIATED(amat%mat)
PRINT*, 'rank of mata', amat%mat%rank
CALL mstruct(amat%mat, dists)
WRITE(*,'(A/(20i4))') 'dists from MSTRUCT=', dists
!
t0 = mpi_wtime(0)
CALL cds2mumps(amat_cds, amat)
CALL factor(amat, debug=.FALSE.)
sol = rhs
CALL bsolve(amat, sol)
WRITE(*,'(a,1pe12.3)') 'Direct solve time (s)', mpi_wtime()-t0
!
IF(me.EQ.0) THEN
WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
END IF
!
! Clean up
!
DEALLOCATE(rhs)
DEALLOCATE(sol)
DEALLOCATE(arow)
CALL destroy(amat)
CALL mpi_finalize(ierr)
CONTAINS
SUBROUTINE mstruct(mat, dists)
TYPE(spmat), INTENT(in) :: mat
INTEGER, ALLOCATABLE, INTENT(inout) :: dists(:)
TYPE(elt), POINTER :: t
INTEGER :: n, i, j0
j0 = LBOUND(dists,1)
n = mat%rank
PRINT*, 'rank of mat', n
DO i=1,n ! scan the matrix rows
t => mat%row(i)%row0
DO WHILE(ASSOCIATED(t)) ! walk thru the linked list row(i)
j = t%index
IF(ABS(t%val) .LE. EPSILON(0.0d0)) THEN
dists(j0) = t%index-i ! distance from main diag
j0 = j0+1
END IF
t => t%next
END DO
END DO
END SUBROUTINE mstruct
INTEGER FUNCTION numb(i,j)
!
! One-dimensional numbering
! Number first x then y
!
INTEGER, INTENT(in) :: i, j
numb = (j-1)*nx + i
END FUNCTION numb
END PROGRAM main
c4science · Help