Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/zssmp_ex1.f
03540012fedefeature/fft-split

Last Change
Raw File

zssmp_ex1.f
View Options

!>
!> @file zssmp_ex1.f
!>
!> @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 wsmp_bsplines
USE pardiso_bsplines
IMPLICIT NONE
c
INTEGER :: n=9
INTEGER ia(10)
INTEGER ja(29)
COMPLEX*16 avals(29)
COMPLEX*16 b(9), sol(9), arow(9)
c
type(zwsmp_mat) :: mat
integer :: i, k
c
DATA ia /1,5,9,13,17,21,25,27,29,30/
data ja
1 /1, 3, 7, 8,
2 2, 3, 8, 9,
3 3, 7, 8, 9,
4 4, 6, 7, 8,
5 5, 6, 8, 9,
6 6, 7, 8, 9,
7 7, 8,
8 8, 9,
9 9/
data avals
1 /(14.d0,0.d0), (-1.d0,-5.d-2), (-1.d0,-5.0d-2), (-3.d0,-1.5d-1),
2 (14.d0,0.d0), (-1.d0,-5.d-2), (-3.d0,-1.5d-1), (-1.d0,-5.0d-2),
3 (16.d0,0.d0), (-2.d0,-1.d-1), (-4.d0,-2.0d-1), (-2.d0,-1.0d-1),
4 (14.d0,0.d0), (-1.d0,-5.d-2), (-1.d0,-5.0d-2), (-3.d0,-1.5d-1),
5 (14.d0,0.d0), (-1.d0,-5.d-2), (-3.d0,-1.5d-1), (-1.d0,-5.0d-2),
6 (16.d0,0.d0), (-2.d0,-1.d-1), (-4.d0,-2.0d-1), (-2.d0,-1.0d-1),
7 (16.d0,0.d0), (-4.d0,-2.d-1),
8 (71.d0,0.d0), (-4.d0,-2.d-1),
9 (16.d0,0.d0)/
c
call init(n, 1, mat, nlherm=.true., nlpos=.true.)
do i=1,n
do k=ia(i),ia(i+1)-1
call putele(mat, i, ja(k), avals(k))
end do
end do
c
print*, 'The RHS before tomat:'
do i = 1, n
call getrow(mat,i, arow)
b(i) = sum(arow)
print *,i,' : ',b(i)
end do
c
call factor(mat)
c
write(*,'(a/(20f6.2))') 'avals', avals
write(*,'(a/(20f6.2))') 'mat%val', mat%val
print*, 'diff of val', cnorm2(avals-mat%val)
print*, 'diff of ia', ia-mat%irow
print*,' diff ja', ja-mat%cols
c
print*, 'Check getrow'
do i = 1, n
call getrow(mat,i, arow)
write(*,'(i3,": ",(20f6.2))') i, arow(i:n)
end do
c
print*, 'The RHS:'
do i = 1, n
call getrow(mat,i, arow)
b(i) = sum(arow)
print *,i,' : ',b(i)
end do
call bsolve(mat,b,sol)
print *,'Norm of Residual = ',mat%p%dparm(7)
print *,'The solution of the system is as follows:'
do i = 1, n
print *,i,' : ',sol(i)
end do
print*, 'Residue =', cnorm2(vmx(mat,sol)-b)
c
contains
FUNCTION cnorm2(x)
DOUBLE COMPLEX, INTENT(in) :: x(:)
DOUBLE PRECISION :: cnorm2
cnorm2 = SQRT(DOT_PRODUCT(x,x))
END FUNCTION cnorm2
c
END PROGRAM main
c4science · Help