Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/factrb.f90
3324689e437cfeature/refactoring

Last Change
Raw File

factrb.f90
View Options

!>
!> @file factrb.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 factrb ( w, ipivot, d, nrow, ncol, last, iflag )
!*************************************************************************
!
!! FACTRB constructs a partial PLU factorization.
!
! Discussion:
!
! This factorization corresponds to steps 1 through LAST in Gauss
! elimination for the matrix W of order ( NROW, NCOL ), using
! pivoting of scaled rows.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input/output, real ( kind = 8 ) W(NROW,NCOL); on input, contains the
! matrix to be partially factored; on output, the partial factorization.
!
! Output, integer IPIVOT(NROW), contains a record of the pivoting
! strategy used; row IPIVOT(I) is used during the I-th elimination step,
! for I = 1, ..., LAST.
!
! Workspace, real ( kind = 8 ) D(NROW), used to store the maximum entry
! in each row.
!
! Input, integer NROW, the number of rows of W.
!
! Input, integer NCOL, the number of columns of W.
!
! Input, integer LAST, the number of elimination steps to be carried out.
!
! Input/output, integer IFLAG. On output, equals the input value
! times (-1)**(number of row interchanges during the factorization
! process), in case no zero pivot was encountered.
! Otherwise, iflag=0 on output.
!
implicit none
integer ncol
integer nrow
real ( kind = 8 ) awikdi
real ( kind = 8 ) colmax
real ( kind = 8 ) d(nrow)
integer i
integer iflag
integer ipivi
integer ipivk
integer ipivot(nrow)
integer j
integer k
integer kp1
integer last
real ( kind = 8 ) ratio
real ( kind = 8 ) rowmax
real ( kind = 8 ) w(nrow,ncol)
!
! Initialize IPIVOT and D.
!
do i = 1, nrow
ipivot(i) = i
end do
do i = 1, nrow
rowmax = 0.0D+00
do j = 1, ncol
rowmax = max ( rowmax, abs ( w(i,j) ) )
end do
if ( rowmax == 0.0D+00 ) then
iflag = 0
return
end if
d(i) = rowmax
end do
!
! Gauss elimination with pivoting of scaled rows, loop over k=1,.,last
!
k = 1
!
! As pivot row for k-th step, pick among the rows not yet used,
! i.e., from rows ipivot(k),...,ipivot(nrow), the one whose k-th
! entry (compared to the row size) is largest. then, if this row
! does not turn out to be row ipivot(k), redefine ipivot(k) ap-
! propriately and record this interchange by changing the sign
! of IFLAG.
!
30 continue
ipivk = ipivot(k)
if ( k == nrow ) then
if ( abs(w(ipivk,nrow))+d(ipivk) <= d(ipivk) ) then
iflag = 0
end if
return
end if
j = k
kp1 = k+1
colmax = abs(w(ipivk,k))/d(ipivk)
!
! Find the largest pivot
!
do i = kp1, nrow
ipivi = ipivot(i)
awikdi = abs(w(ipivi,k)) / d(ipivi)
if ( colmax < awikdi ) then
colmax = awikdi
j = i
end if
end do
if ( j /= k ) then
ipivk = ipivot(j)
ipivot(j) = ipivot(k)
ipivot(k) = ipivk
iflag = -iflag
end if
!
! If pivot element is too small in absolute value, declare
! matrix to be noninvertible and quit.
!
if ( abs(w(ipivk,k))+d(ipivk) <= d(ipivk) ) then
iflag = 0
return
end if
!
! Otherwise, subtract the appropriate multiple of the pivot
! row from remaining rows, i.e., the rows ipivot(k+1),...,
! ipivot(nrow), to make k-th entry zero. save the multiplier in
! its place.
!
do i = kp1, nrow
ipivi = ipivot(i)
w(ipivi,k) = w(ipivi,k)/w(ipivk,k)
ratio = -w(ipivi,k)
do j = kp1, ncol
w(ipivi,j) = ratio*w(ipivk,j)+w(ipivi,j)
end do
end do
k = kp1
!
! Check for having reached the next block.
!
if ( k <= last ) then
go to 30
end if
return
end
c4science · Help