Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/pppack/slvblk.f90
a1c89766d1dafeature/fft-split

Last Change
Raw File

slvblk.f90
View Options

!>
!> @file slvblk.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 slvblk ( bloks, integs, nbloks, b, ipivot, x, iflag )
!*************************************************************************
!
!! SLVBLK solves the almost block diagonal linear system A*x=b.
!
! Discussion:
!
! Such almost block diagonal matrices arise naturally in piecewise
! polynomial interpolation or approximation and in finite element
! methods for two-point boundary value problems. The PLU factorization
! method is implemented here to take advantage of the special structure
! of such systems for savings in computing time and storage requirements.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! bloks a one-dimenional array, of length
! sum( integs(1,i)*integs(2,i) ; i=1,nbloks )
! on input, contains the blocks of the almost block diagonal
! matrix a . the array integs (see below and the example)
! describes the block structure.
! on output, contains correspondingly the plu factorization
! of a (if iflag /= 0). certain of the entries into bloks
! are arbitrary (where the blocks overlap).
!
! integs integer array description of the block structure of a .
! integs(1,i)=no. of rows of block i = nrow
! integs(2,i)=no. of colums of block i = ncol
! integs(3,i)=no. of elim. steps in block i = last
! i =1,2,...,nbloks
! the linear system is of order
! n = sum ( integs(3,i) , i=1,...,nbloks ),
! but the total number of rows in the blocks is
! nbrows=sum( integs(1,i) ; i = 1,...,nbloks)
!
! nbloks number of blocks
! b right side of the linear system, array of length nbrows.
! certain of the entries are arbitrary, corresponding to
! rows of the blocks which overlap (see block structure and
! the example below).
! ipivot on output, integer array containing the pivoting sequence
! used. length is nbrows
! x on output, contains the computed solution (if iflag /= 0)
! length is n.
! iflag on output, integer
! =(-1)**(no. of interchanges during factorization)
! if a is invertible
! =0 if a is singular
!
! auxiliary programs
! fcblok (bloks,integs,nbloks,ipivot,scrtch,iflag) factors the matrix
! a , and is used for this purpose in slvblk. its arguments
! are as in slvblk, except for
! scrtch=a work array of length max(integs(1,i)).
!
! sbblok (bloks,integs,nbloks,ipivot,b,x) solves the system a*x=b
! once a is factored. this is done automatically by slvblk
! for one right side b, but subsequent solutions may be
! obtained for additional b-vectors. the arguments are all
! as in slvblk.
!
! dtblok (bloks,integs,nbloks,ipivot,iflag,detsgn,detlog) computes the
! determinant of a once slvblk or fcblok has done the fact-
! orization.the first five arguments are as in slvblk.
! detsgn =sign of the determinant
! detlog =natural log of the determinant
!
! block structure of a
! the nbloks blocks are stored consecutively in the array bloks .
! the first block has its (1,1)-entry at bloks(1), and, if the i-th
! block has its (1,1)-entry at bloks(index(i)), then
! index(i+1)=index(i) + nrow(i)*ncol(i) .
! the blocks are pieced together to give the interesting part of a
! as follows. for i=1,2,...,nbloks-1, the (1,1)-entry of the next
! block (the (i+1)st block ) corresponds to the (last+1,last+1)-entry
! of the current i-th block. recall last=integs(3,i) and note that
! this means that
! a. every block starts on the diagonal of a .
! b. the blocks overlap (usually). the rows of the (i+1)st block
! which are overlapped by the i-th block may be arbitrarily de-
! fined initially. they are overwritten during elimination.
! the right side for the equations in the i-th block are stored cor-
! respondingly as the last entries of a piece of b of length nrow
! (= integs(1,i)) and following immediately in b the corresponding
! piece for the right side of the preceding block, with the right side
! for the first block starting at b(1) . in this, the right side for
! an equation need only be specified once on input, in the first block
! in which the equation appears.
!
! example and test driver
! the test driver for this package contains an example, a linear
! system of order 11, whose nonzero entries are indicated in the fol-
! lowing schema by their row and column index modulo 10. next to it
! are the contents of the integs arrray when the matrix is taken to
! be almost block diagonal with nbloks=5, and below it are the five
! blocks.
!
! nrow1=3, ncol1 = 4
! 11 12 13 14
! 21 22 23 24 nrow2=3, ncol2 = 3
! 31 32 33 34
! last1=2 43 44 45
! 53 54 55 nrow3=3, ncol3 = 4
! last2=3 66 67 68 69 nrow4 = 3, ncol4 = 4
! 76 77 78 79 nrow5=4, ncol5 = 4
! 86 87 88 89
! last3=1 97 98 99 90
! last4=1 08 09 00 01
! 18 19 10 11
! last5=4
!
! actual input to bloks shown by rows of blocks of a .
! (the ** items are arbitrary, this storage is used by slvblk)
!
! 11 12 13 14 / ** ** ** / 66 67 68 69 / ** ** ** ** / ** ** ** **
! 21 22 23 24 / 43 44 45 / 76 77 78 79 / ** ** ** ** / ** ** ** **
! 31 32 33 34/ 53 54 55/ 86 87 88 89/ 97 98 99 90/ 08 09 00 01
! 18 19 10 11
!
! index=1 index = 13 index = 22 index = 34 index = 46
!
! actual right side values with ** for arbitrary values
! b1 b2 b3 ** b4 b5 b6 b7 b8 ** ** b9 ** ** b10 b11
!
! (it would have been more efficient to combine block 3 with block 4)
!
implicit none
integer nbloks
real ( kind = 8 ) b(*)
real ( kind = 8 ) bloks(*)
integer iflag
integer integs(3,nbloks)
integer ipivot(*)
real ( kind = 8 ) x(*)
!
! In the call to FCBLOK, X is used for temporary storage.
!
call fcblok ( bloks, integs, nbloks, ipivot, x, iflag )
if ( iflag == 0 ) then
return
end if
call sbblok ( bloks, integs, nbloks, ipivot, b, x )
return
end
c4science · Help