Diffusion SPClibs (feature/fft-split) Browse

fft-split

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

spclibs/pppack/eqblok.f90c723281d03b6feature/fft-split

eqblok.f90View Options

spclibs/pppack/eqblok.f90
c723281d03b6feature/fft-split

eqblok.f90
View Options