Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/chol1d.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

chol1d.f90
View Options

!>
!> @file chol1d.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>
!>
subroutine chol1d ( p, v, qty, npoint, u, qu )
!*************************************************************************
!
!! CHOL1D sets up and solves linear systems needed by SMOOTH.
!
! Discussion:
!
! This routine constructs the upper three diagonals of
!
! V(I,J), I=2 to NPOINT-1, J=1,3,
!
! of the matrix
!
! 6*(1-P)*Q-transpose*(D**2)*Q + P*R.
!
! It then computes its L*L' decomposition and stores it also
! in V, then applies forward and backsubstitution to the right side
!
! Q'*Y
!
! in QTY to obtain the solution in U.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) P, ?
!
! ?put, real ( kind = 8 ) V(NPOINT,7), ?
!
! ?put, real ( kind = 8 ) QTY(NPOINT), ?
!
! Input, integer NPOINT, ?
!
! Output, real ( kind = 8 ) U(NPOINT), the solution.
!
! Output, real ( kind = 8 ) QU(NPOINT), the value of Q * U.
!
implicit none
integer npoint
integer i
real ( kind = 8 ) p
real ( kind = 8 ) qty(npoint)
real ( kind = 8 ) qu(npoint)
real ( kind = 8 ) u(npoint)
real ( kind = 8 ) v(npoint,7)
real ( kind = 8 ) prev
real ( kind = 8 ) ratio
real ( kind = 8 ) six1mp
real ( kind = 8 ) twop
!
! Construct 6*(1-p)*q'*(d**2)*q + p*r
!
six1mp = 6.0D+00 * ( 1.0D+00 - p )
twop = 2.0D+00 * p
do i = 2, npoint-1
v(i,1) = six1mp * v(i,5)+twop*(v(i-1,4)+v(i,4))
v(i,2) = six1mp * v(i,6)+p*v(i,4)
v(i,3) = six1mp * v(i,7)
end do
if ( npoint < 4 ) then
u(1) = 0.0D+00
u(2) = qty(2) / v(2,1)
u(3) = 0.0D+00
!
! Factorization
!
else
do i = 2, npoint-2
ratio = v(i,2)/v(i,1)
v(i+1,1) = v(i+1,1)-ratio*v(i,2)
v(i+1,2) = v(i+1,2)-ratio*v(i,3)
v(i,2) = ratio
ratio = v(i,3)/v(i,1)
v(i+2,1) = v(i+2,1)-ratio*v(i,3)
v(i,3) = ratio
end do
!
! Forward substitution
!
u(1) = 0.0D+00
v(1,3) = 0.0D+00
u(2) = qty(2)
do i = 2, npoint-2
u(i+1) = qty(i+1)-v(i,2)*u(i)-v(i-1,3)*u(i-1)
end do
!
! Back substitution.
!
u(npoint) = 0.0D+00
u(npoint-1) = u(npoint-1) / v(npoint-1,1)
do i = npoint-2, 2, -1
u(i) = u(i)/v(i,1)-u(i+1)*v(i,2)-u(i+2)*v(i,3)
end do
end if
!
! Construct Q*U.
!
prev = 0.0D+00
do i = 2, npoint
qu(i) = (u(i)-u(i-1))/v(i-1,4)
qu(i-1) = qu(i)-prev
prev = qu(i)
end do
qu(npoint) = -qu(npoint)
return
end
c4science · Help