Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/eqblok.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

eqblok.f90
View Options

!>
!> @file eqblok.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 eqblok ( t, n, kpm, work1, work2, bloks, lenblk, integs, &
nbloks, b )
!*************************************************************************
!
!! EQBLOK is to be called in COLLOC.
!
! Method:
!
! Each breakpoint interval gives rise to a block in the linear system.
! this block is determined by the K collocation equations in the interval
! with the side conditions (if any) in the interval interspersed ap-
! propriately, and involves the kpm b-splines having the interval in
! their support. correspondingly, such a block has nrow=k+isidel
! rows, with isidel=number of side conditions in this and the prev-
! ious intervals, and ncol=kpm columns.
!
! Further, because the interior knots have multiplicity k, we can
! carry out (in slvblk) k elimination steps in a block before pivot-
! ing might involve an equation from the next block. in the last block,
! of course, all kpm elimination steps will be carried out (in slvblk).
!
! see the detailed comments in the solveblok package for further in-
! formation about the almost block diagonal form used here.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! input
!
! Input, real ( kind = 8 ) T(N+KPM), the knot sequence.
!
! Input, integer N, the dimension of the approximating spline space,
! that is, the order of the linear system to be constructed.
!
! Input, integer KPM, = K + M, the order of the approximating spline.
!
! Input, integer LENBLK, the maximum length of the array BLOKS,
! as allowed by the dimension statement in COLLOC.
!
! work areas
!
! work1 used in putit, of size (kpm,kpm)
! work2 used in putit, of size (kpm,m+1)
!
! output
!
! bloks the coefficient matrix of the linear system, stored in al-
! most block diagonal form, of size
! kpm*sum(integs(1,i) , i=1,...,nbloks)
!
! integs an integer array, of size (3,nbloks), describing the block
! structure.
! integs(1,i) = number of rows in block i
! integs(2,i) = number of columns in block i
! integs(3,i) = number of elimination steps which can be
! carried out in block i before pivoting might
! bring in an equation from the next block.
!
! nbloks number of blocks, equals number of polynomial pieces
!
! b the right side of the linear system, stored corresponding to the
! almost block diagonal form, of size sum(integs(1,i) , i=1,...,
! nbloks).
!
implicit none
integer kpm
integer n
real ( kind = 8 ) b(*)
real ( kind = 8 ) bloks(*)
integer i
integer index
integer indexb
integer integs(3,*)
integer iside
integer isidel
integer itermx
integer k
integer left
integer lenblk
integer m
integer nbloks
integer nrow
real ( kind = 8 ) rho
real ( kind = 8 ) t(n+kpm)
real ( kind = 8 ) work1(kpm,kpm)
real ( kind = 8 ) work2(kpm,*)
real ( kind = 8 ) xside
common /side/ m,iside,xside(10)
common /other/ itermx,k,rho(19)
index = 1
indexb = 1
i = 0
iside = 1
do left = kpm, n, k
i = i + 1
!
! Determine integs(.,i)
!
integs(2,i) = kpm
if ( n <= left ) then
integs(3,i) = kpm
isidel = m
go to 30
end if
integs(3,i) = k
!
! At this point, iside-1 gives the number of side conditions
! incorporated so far. adding to this the side conditions in the
! current interval gives the number isidel .
!
isidel = iside - 1
do
if ( isidel == m ) then
exit
end if
if ( t(left+1) <= xside(isidel+1) ) then
exit
end if
isidel = isidel + 1
end do
30 continue
nrow = k + isidel
integs(1,i) = nrow
!
! The detailed equations for this block are generated and put
! together in PUTIT.
!
if ( lenblk < index + nrow * kpm - 1 ) then
write ( *, * ) ' '
write ( *, * ) 'EQBLOK - Fatal error!'
write ( *, * ) ' The dimension of BLOKS is too small.'
write ( *, * ) ' LENBLK = ', lenblk
stop
end if
call putit ( t, kpm, left, work1, work2, bloks(index), nrow, b(indexb) )
index = index + nrow * kpm
indexb = indexb + nrow
end do
nbloks = i
return
end
c4science · Help