Balances the matrix so that the rows with zero entries off the diagonal are isolated and the remaining columns and rows are resized to have one norm close to 1.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | n | Dimension of |
||
real(kind=double) | :: | mat(0:n-1,0:n-1) |
|
|||
real(kind=double) | :: | scal(0:n-1) | Contains isolated eigenvalue in the positions 0- |
|||
integer, | intent(out) | :: | low | |||
integer, | intent(out) | :: | high |
subroutine balance(n, mat, scal, low, high )
!+Balances the matrix so that the rows with zero entries
!+off the diagonal are isolated and the remaining columns and rows
!+are resized to have one norm close to 1.
integer, intent(in) :: n
!+ Dimension of `mat`
real(double) :: mat(0:n-1,0:n-1)
!+ `n`x`n` scaled matrix
real(double) :: scal(0:n-1)
!+ Contains isolated eigenvalue in the positions 0-`low` and `high`-`n`-1
!+ its other components contain the scaling factors for transforming `mat`
integer, intent(out) :: high
integer, intent(out) :: low
integer, parameter :: basis = 2
real(double) :: b2, r, c, f, g, s
integer :: m, k, i, j, iter
scal=0.d0
b2 = basis * basis
m = 0
k = n - 1
iter=1
do while(iter==1)
iter = 0
do j = k, 0, -1
r = ZERO
do i = 0, k
if (i.ne.j) r = r + DABS(mat(j,i))
enddo
if (r == ZERO) then
scal(k) = j
if (j.ne.k) then
do i = 0, k
call RSWAP(mat(i,j), mat(i,k))
enddo
do i = m, n-1
call RSWAP(mat(j,i), mat(k,i))
enddo
endif
k=k-1
iter = 1
endif
enddo !j loop
enddo !while iter=1
iter=1
do while (iter==1)
iter = 0
do j = m, k
c = ZERO
do i = m, k
if (i.ne.j) c = c + DABS(mat(i,j))
enddo
if (c == ZERO) then
scal(m) = j
if (j.ne.m) then
do i = 0, k
call RSWAP(mat(i,j), mat(i,m))
enddo
do i = m, n-1
call RSWAP(mat(j,i), mat(m,i))
enddo
endif
m = m + 1
iter = 1
endif
enddo !j loop
enddo !while iter=1
low = m
high = k
do i = m, k
scal(i) = ONE
enddo
iter=1
do while (iter==1)
iter = 0
do i = m, k
c=ZERO; r=ZERO
do j = m, k
if (j.ne.i) then
c = c + DABS(mat(j,i))
r = r + DABS(mat(i,j))
endif
enddo
g = r / basis
f = ONE
s = c + r
do while (c < g)
f = f * basis
c = c * b2
enddo
g = r * basis
do while (c >= g)
f = f / basis
c = c / b2
enddo
if ((c + r) / f < 0.95 * s) then
g = ONE / f
scal(i) = scal(i) * f
iter = 1
do j = m, n-1
mat(i,j) = mat(i,j) * g
enddo
do j = 0, k
mat(j,i) = mat(j,i) * f
enddo
endif
enddo !i loop
enddo !while iter=1
return
end subroutine balance