psgebal.f

Go to the documentation of this file.

      SUBROUTINE psgebal( JOB, N, A, DESCA, ILO, IHI, SCALE, INFO )
*
*     Contribution from the Department of Computing Science and HPC2N,
*     Umea University, Sweden
*
*  -- ScaLAPACK computational routine (version 2.0.1) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     Univ. of Colorado Denver and University of California, Berkeley.
*     January, 2012
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          JOB
      INTEGER            IHI, ILO, INFO, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      REAL               A( * ), SCALE( * )
*     ..
*
*  Purpose
*  =======
*
*  PSGEBAL balances a general real matrix A.  This involves, first,
*  permuting A by a similarity transformation to isolate eigenvalues
*  in the first 1 to ILO-1 and last IHI+1 to N elements on the
*  diagonal; and second, applying a diagonal similarity transformation
*  to rows and columns ILO to IHI to make the rows and columns as
*  close in norm as possible.  Both steps are optional.
*
*  Balancing may reduce the 1-norm of the matrix, and improve the
*  accuracy of the computed eigenvalues and/or eigenvectors.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*
*  Arguments
*  =========
*
*  JOB     (global input) CHARACTER*1
*          Specifies the operations to be performed on A:
*          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
*                  for i = 1,...,N;
*          = 'P':  permute only;
*          = 'S':  scale only;
*          = 'B':  both permute and scale.
*
*  N       (global input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  A       (local input/output) REAL             array, dimension
*          (DESCA(LLD_,LOCc(N))
*          On entry, the input matrix A.
*          On exit,  A is overwritten by the balanced matrix.
*          If JOB = 'N', A is not referenced.
*          See Further Details.
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  ILO     (global output) INTEGER
*  IHI     (global output) INTEGER
*          ILO and IHI are set to integers such that on exit
*          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
*          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
*
*  SCALE   (global output) REAL             array, dimension (N)
*          Details of the permutations and scaling factors applied to
*          A.  If P(j) is the index of the row and column interchanged
*          with row and column j and D(j) is the scaling factor
*          applied to row and column j, then
*          SCALE(j) = P(j)    for j = 1,...,ILO-1
*                   = D(j)    for j = ILO,...,IHI
*                   = P(j)    for j = IHI+1,...,N.
*          The order in which the interchanges are made is N to IHI+1,
*          then 1 to ILO-1.
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*
*  Further Details
*  ===============
*
*  The permutations consist of row and column interchanges which put
*  the matrix in the form
*
*             ( T1   X   Y  )
*     P A P = (  0   B   Z  )
*             (  0   0   T2 )
*
*  where T1 and T2 are upper triangular matrices whose eigenvalues lie
*  along the diagonal.  The column indices ILO and IHI mark the starting
*  and ending columns of the submatrix B. Balancing consists of applying
*  a diagonal similarity transformation inv(D) * B * D to make the
*  1-norms of each row of B and its corresponding column nearly equal.
*  The output matrix is
*
*     ( T1     X*D          Y    )
*     (  0  inv(D)*B*D  inv(D)*Z ).
*     (  0      0           T2   )
*
*  Information about the permutations P and the diagonal matrix D is
*  returned in the vector SCALE.
*
*  This subroutine is based on the EISPACK routine BALANC. In principle,
*  the parallelism is extracted by using PBLAS and BLACS routines for
*  the permutation and balancing.
*
*  Modified by Tzu-Yi Chen, Computer Science Division, University of
*    California at Berkeley, USA
*
*  Parallel version by Robert Granat and Meiyue Shao, Department of
*    Computing Science and HPC2N, Umea University, Sweden
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
     $                     ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                     rsrc_ = 7, csrc_ = 8, lld_ = 9 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
      REAL               SCLFAC
      parameter( sclfac = 2.0e+0 )
      REAL               FACTOR
      parameter( factor = 0.95e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOCONV
      INTEGER            I, ICA, IEXC, IRA, J, K, L, M, LLDA,
     $                   ICTXT, NPROW, NPCOL, MYROW, MYCOL, II, JJ,
     $                   ARSRC, ACSRC
      REAL               C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
     $                   SFMIN2, ELEM
*     ..
*     .. Local Arrays ..
      REAL               CR( 2 )
*     ..
*     .. External Functions ..
      LOGICAL            SISNAN, LSAME
      INTEGER            IDAMAX
      REAL               SLAMCH
      EXTERNAL           sisnan, lsame, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           psscal, psswap, psamax, pxerbla,
     $                   blacs_gridinfo, chk1mat, sgsum2d,
     $                   infog2l, pselget
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
      info = 0
      ictxt = desca( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
*     Test the input parameters.
*
      IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
     $    .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE
         CALL chk1mat( n, 2, n, 2, 1, 1, desca, 4, info )
      END IF
      IF( info.NE.0 ) THEN
         CALL pxerbla( ictxt, 'PSGEBAL', -info )
         RETURN
      END IF
*
*     Extract local leading dimension of A.
*
      llda = desca( lld_ )
*
      k = 1
      l = n
*
      IF( n.EQ.0 )
     $   GO TO 210
*
      IF( lsame( job, 'N' ) ) THEN
         DO 10 i = 1, n
            scale( i ) = one
   10    CONTINUE
         GO TO 210
      END IF
*
      IF( lsame( job, 'S' ) )
     $   GO TO 120
*
*     Permutation to isolate eigenvalues if possible.
*
      GO TO 50
*
*     Row and column exchange.
*
   20 CONTINUE
      scale( m ) = j
      IF( j.EQ.m )
     $   GO TO 30
*
      CALL psswap( l, a, 1, j, desca, 1, a, 1, m, desca, 1 )
      CALL psswap( n-k+1, a, j, k, desca, desca(m_), a, m, k, desca,
     $             desca(m_) )
*
   30 CONTINUE
      GO TO ( 40, 80 )iexc
*
*     Search for rows isolating an eigenvalue and push them down.
*
   40 CONTINUE
      IF( l.EQ.1 )
     $   GO TO 210
      l = l - 1
*
   50 CONTINUE
      DO 70 j = l, 1, -1
*
         DO 60 i = 1, l
            IF( i.EQ.j )
     $         GO TO 60
*
*           All processors need the information to make correct decisions.
*
            CALL pselget( 'All', '1-Tree', elem, a, j, i, desca )
            IF( elem.NE.zero )
     $         GO TO 70
   60    CONTINUE
*
         m = l
         iexc = 1
         GO TO 20
   70 CONTINUE
*
      GO TO 90
*
*     Search for columns isolating an eigenvalue and push them left.
*
   80 CONTINUE
      k = k + 1
*
   90 CONTINUE
      DO 110 j = k, l
*
         DO 100 i = k, l
            IF( i.EQ.j )
     $         GO TO 100
*
*           All processors need the information to make correct decisions.
*
            CALL pselget( 'All', '1-Tree', elem, a, i, j, desca )
            IF( elem.NE.zero )
     $         GO TO 110
  100    CONTINUE
*
         m = k
         iexc = 2
         GO TO 20
  110 CONTINUE
*
  120 CONTINUE
      DO 130 i = k, l
         scale( i ) = one
  130 CONTINUE
*
      IF( lsame( job, 'P' ) )
     $   GO TO 210
*
*     Balance the submatrix in rows K to L.
*
*     Iterative loop for norm reduction.
*
      sfmin1 = slamch( 'S' ) / slamch( 'P' )
      sfmax1 = one / sfmin1
      sfmin2 = sfmin1*sclfac
      sfmax2 = one / sfmin2
  140 CONTINUE
      noconv = .false.
*
      DO 200 i = k, l
         c = zero
         r = zero
*
*        Compute local partial values of R and C in parallel and combine
*        with a call to the BLACS global summation routine distributing
*        information to all processors.
*
         DO 150 j = k, l
            IF( j.EQ.i )
     $         GO TO 150
            CALL infog2l( j, i, desca, nprow, npcol, myrow,
     $                    mycol, ii, jj, arsrc, acsrc )
            IF( myrow.EQ.arsrc .AND. mycol.EQ.acsrc ) THEN
               c = c + abs( a( ii + (jj-1)*llda ) )
            END IF
            CALL infog2l( i, j, desca, nprow, npcol, myrow,
     $                    mycol, ii, jj, arsrc, acsrc )
            IF( myrow.EQ.arsrc .AND. mycol.EQ.acsrc ) THEN
               r = r + abs( a( ii + (jj-1)*llda ) )
            END IF
  150    CONTINUE
         cr( 1 ) = c
         cr( 2 ) = r
         CALL sgsum2d( ictxt, 'All', '1-Tree', 2, 1, cr, 2, -1, -1 )
         c = cr( 1 )
         r = cr( 2 )
*
*        Find global maximum absolute values and indices in parallel.
*
         CALL psamax( l, ca, ica, a, 1, i, desca, 1 )
         CALL psamax( n-k+1, ra, ira, a, i, k, desca, desca(m_) )
*
*        Guard against zero C or R due to underflow.
*
         IF( c.EQ.zero .OR. r.EQ.zero )
     $      GO TO 200
         g = r / sclfac
         f = one
         s = c + r
  160    CONTINUE
         IF( c.GE.g .OR. max( f, c, ca ).GE.sfmax2 .OR.
     $       min( r, g, ra ).LE.sfmin2 )GO TO 170
         IF( sisnan( c+f+ca+r+g+ra ) ) THEN
*
*           Exit if NaN to avoid infinite loop
*
            info = -3
            CALL pxerbla( ictxt, 'PDGEBAL', -info )
            RETURN
         END IF
         f = f*sclfac
         c = c*sclfac
         ca = ca*sclfac
         r = r / sclfac
         g = g / sclfac
         ra = ra / sclfac
         GO TO 160
*
  170    CONTINUE
         g = c / sclfac
  180    CONTINUE
         IF( g.LT.r .OR. max( r, ra ).GE.sfmax2 .OR.
     $       min( f, c, g, ca ).LE.sfmin2 )GO TO 190
         f = f / sclfac
         c = c / sclfac
         g = g / sclfac
         ca = ca / sclfac
         r = r*sclfac
         ra = ra*sclfac
         GO TO 180
*
*        Now balance.
*
  190    CONTINUE
         IF( ( c+r ).GE.factor*s )
     $      GO TO 200
         IF( f.LT.one .AND. scale( i ).LT.one ) THEN
            IF( f*scale( i ).LE.sfmin1 )
     $         GO TO 200
         END IF
         IF( f.GT.one .AND. scale( i ).GT.one ) THEN
            IF( scale( i ).GE.sfmax1 / f )
     $         GO TO 200
         END IF
         g = one / f
         scale( i ) = scale( i )*f
         noconv = .true.
*
         CALL psscal( n-k+1, g, a, i, k, desca, desca(m_) )
         CALL psscal( l, f, a, 1, i, desca, 1 )
*
  200 CONTINUE
*
      IF( noconv )
     $   GO TO 140
*
  210 CONTINUE
      ilo = k
      ihi = l
*
      RETURN
*
*     End of PSGEBAL
*
      END