balance Subroutine

  1. eigensystem.f90
  2. eigensystem
  3. balance

public 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.

Arguments

Type IntentOptional AttributesName
integer, intent(in) :: n

Dimension of mat
real(kind=double) :: mat(0:n-1,0:n-1)

nxn scaled matrix
real(kind=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) :: low
integer, intent(out) :: high

Called by

proc~~balance~~CalledByGraph proc~balance balance proc~eigen eigen proc~eigen->proc~balance proc~colrad colrad proc~colrad->proc~eigen proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~colrad proc~read_equilibrium read_equilibrium proc~read_equilibrium->proc~colrad proc~attenuate attenuate proc~attenuate->proc~colrad proc~pnpa_f pnpa_f proc~pnpa_f->proc~attenuate proc~pnpa_mc pnpa_mc proc~pnpa_mc->proc~attenuate proc~npa_f npa_f proc~npa_f->proc~attenuate program~fidasim fidasim program~fidasim->proc~fida_weights_los program~fidasim->proc~read_equilibrium program~fidasim->proc~pnpa_f program~fidasim->proc~pnpa_mc program~fidasim->proc~npa_f proc~npa_mc npa_mc program~fidasim->proc~npa_mc proc~npa_weights npa_weights program~fidasim->proc~npa_weights proc~npa_mc->proc~attenuate proc~npa_weights->proc~attenuate

Contents

Source Code

balance

Source Code

  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