home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / lapack / zgebal.f < prev    next >
Text File  |  1996-09-28  |  9KB  |  329 lines
  1.       SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          JOB
  10.       INTEGER            IHI, ILO, INFO, LDA, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   SCALE( * )
  14.       COMPLEX*16         A( LDA, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZGEBAL balances a general complex matrix A.  This involves, first,
  21. *  permuting A by a similarity transformation to isolate eigenvalues
  22. *  in the first 1 to ILO-1 and last IHI+1 to N elements on the
  23. *  diagonal; and second, applying a diagonal similarity transformation
  24. *  to rows and columns ILO to IHI to make the rows and columns as
  25. *  close in norm as possible.  Both steps are optional.
  26. *
  27. *  Balancing may reduce the 1-norm of the matrix, and improve the
  28. *  accuracy of the computed eigenvalues and/or eigenvectors.
  29. *
  30. *  Arguments
  31. *  =========
  32. *
  33. *  JOB     (input) CHARACTER*1
  34. *          Specifies the operations to be performed on A:
  35. *          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
  36. *                  for i = 1,...,N;
  37. *          = 'P':  permute only;
  38. *          = 'S':  scale only;
  39. *          = 'B':  both permute and scale.
  40. *
  41. *  N       (input) INTEGER
  42. *          The order of the matrix A.  N >= 0.
  43. *
  44. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  45. *          On entry, the input matrix A.
  46. *          On exit,  A is overwritten by the balanced matrix.
  47. *          If JOB = 'N', A is not referenced.
  48. *          See Further Details.
  49. *
  50. *  LDA     (input) INTEGER
  51. *          The leading dimension of the array A.  LDA >= max(1,N).
  52. *
  53. *  ILO     (output) INTEGER
  54. *  IHI     (output) INTEGER
  55. *          ILO and IHI are set to integers such that on exit
  56. *          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
  57. *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  58. *
  59. *  SCALE   (output) DOUBLE PRECISION array, dimension (N)
  60. *          Details of the permutations and scaling factors applied to
  61. *          A.  If P(j) is the index of the row and column interchanged
  62. *          with row and column j and D(j) is the scaling factor
  63. *          applied to row and column j, then
  64. *          SCALE(j) = P(j)    for j = 1,...,ILO-1
  65. *                   = D(j)    for j = ILO,...,IHI
  66. *                   = P(j)    for j = IHI+1,...,N.
  67. *          The order in which the interchanges are made is N to IHI+1,
  68. *          then 1 to ILO-1.
  69. *
  70. *  INFO    (output) INTEGER
  71. *          = 0:  successful exit.
  72. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  73. *
  74. *  Further Details
  75. *  ===============
  76. *
  77. *  The permutations consist of row and column interchanges which put
  78. *  the matrix in the form
  79. *
  80. *             ( T1   X   Y  )
  81. *     P A P = (  0   B   Z  )
  82. *             (  0   0   T2 )
  83. *
  84. *  where T1 and T2 are upper triangular matrices whose eigenvalues lie
  85. *  along the diagonal.  The column indices ILO and IHI mark the starting
  86. *  and ending columns of the submatrix B. Balancing consists of applying
  87. *  a diagonal similarity transformation inv(D) * B * D to make the
  88. *  1-norms of each row of B and its corresponding column nearly equal.
  89. *  The output matrix is
  90. *
  91. *     ( T1     X*D          Y    )
  92. *     (  0  inv(D)*B*D  inv(D)*Z ).
  93. *     (  0      0           T2   )
  94. *
  95. *  Information about the permutations P and the diagonal matrix D is
  96. *  returned in the vector SCALE.
  97. *
  98. *  This subroutine is based on the EISPACK routine CBAL.
  99. *
  100. *  =====================================================================
  101. *
  102. *     .. Parameters ..
  103.       DOUBLE PRECISION   ZERO, ONE
  104.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  105.       DOUBLE PRECISION   SCLFAC
  106.       PARAMETER          ( SCLFAC = 1.0D+1 )
  107.       DOUBLE PRECISION   FACTOR
  108.       PARAMETER          ( FACTOR = 0.95D+0 )
  109. *     ..
  110. *     .. Local Scalars ..
  111.       LOGICAL            NOCONV
  112.       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
  113.       DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
  114.      $                   SFMIN2
  115.       COMPLEX*16         CDUM
  116. *     ..
  117. *     .. External Functions ..
  118.       LOGICAL            LSAME
  119.       INTEGER            IZAMAX
  120.       DOUBLE PRECISION   DLAMCH
  121.       EXTERNAL           LSAME, IZAMAX, DLAMCH
  122. *     ..
  123. *     .. External Subroutines ..
  124.       EXTERNAL           XERBLA, ZDSCAL, ZSWAP
  125. *     ..
  126. *     .. Intrinsic Functions ..
  127.       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
  128. *     ..
  129. *     .. Statement Functions ..
  130.       DOUBLE PRECISION   CABS1
  131. *     ..
  132. *     .. Statement Function definitions ..
  133.       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  134. *     ..
  135. *     .. Executable Statements ..
  136. *
  137. *     Test the input parameters
  138. *
  139.       INFO = 0
  140.       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
  141.      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
  142.          INFO = -1
  143.       ELSE IF( N.LT.0 ) THEN
  144.          INFO = -2
  145.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  146.          INFO = -4
  147.       END IF
  148.       IF( INFO.NE.0 ) THEN
  149.          CALL XERBLA( 'ZGEBAL', -INFO )
  150.          RETURN
  151.       END IF
  152. *
  153.       K = 1
  154.       L = N
  155. *
  156.       IF( N.EQ.0 )
  157.      $   GO TO 210
  158. *
  159.       IF( LSAME( JOB, 'N' ) ) THEN
  160.          DO 10 I = 1, N
  161.             SCALE( I ) = ONE
  162.    10    CONTINUE
  163.          GO TO 210
  164.       END IF
  165. *
  166.       IF( LSAME( JOB, 'S' ) )
  167.      $   GO TO 120
  168. *
  169. *     Permutation to isolate eigenvalues if possible
  170. *
  171.       GO TO 50
  172. *
  173. *     Row and column exchange.
  174. *
  175.    20 CONTINUE
  176.       SCALE( M ) = J
  177.       IF( J.EQ.M )
  178.      $   GO TO 30
  179. *
  180.       CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
  181.       CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
  182. *
  183.    30 CONTINUE
  184.       GO TO ( 40, 80 )IEXC
  185. *
  186. *     Search for rows isolating an eigenvalue and push them down.
  187. *
  188.    40 CONTINUE
  189.       IF( L.EQ.1 )
  190.      $   GO TO 210
  191.       L = L - 1
  192. *
  193.    50 CONTINUE
  194.       DO 70 J = L, 1, -1
  195. *
  196.          DO 60 I = 1, L
  197.             IF( I.EQ.J )
  198.      $         GO TO 60
  199.             IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
  200.      $          ZERO )GO TO 70
  201.    60    CONTINUE
  202. *
  203.          M = L
  204.          IEXC = 1
  205.          GO TO 20
  206.    70 CONTINUE
  207. *
  208.       GO TO 90
  209. *
  210. *     Search for columns isolating an eigenvalue and push them left.
  211. *
  212.    80 CONTINUE
  213.       K = K + 1
  214. *
  215.    90 CONTINUE
  216.       DO 110 J = K, L
  217. *
  218.          DO 100 I = K, L
  219.             IF( I.EQ.J )
  220.      $         GO TO 100
  221.             IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
  222.      $          ZERO )GO TO 110
  223.   100    CONTINUE
  224. *
  225.          M = K
  226.          IEXC = 2
  227.          GO TO 20
  228.   110 CONTINUE
  229. *
  230.   120 CONTINUE
  231.       DO 130 I = K, L
  232.          SCALE( I ) = ONE
  233.   130 CONTINUE
  234. *
  235.       IF( LSAME( JOB, 'P' ) )
  236.      $   GO TO 210
  237. *
  238. *     Balance the submatrix in rows K to L.
  239. *
  240. *     Iterative loop for norm reduction
  241. *
  242.       SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
  243.       SFMAX1 = ONE / SFMIN1
  244.       SFMIN2 = SFMIN1*SCLFAC
  245.       SFMAX2 = ONE / SFMIN2
  246.   140 CONTINUE
  247.       NOCONV = .FALSE.
  248. *
  249.       DO 200 I = K, L
  250.          C = ZERO
  251.          R = ZERO
  252. *
  253.          DO 150 J = K, L
  254.             IF( J.EQ.I )
  255.      $         GO TO 150
  256.             C = C + CABS1( A( J, I ) )
  257.             R = R + CABS1( A( I, J ) )
  258.   150    CONTINUE
  259.          ICA = IZAMAX( L, A( 1, I ), 1 )
  260.          CA = ABS( A( ICA, I ) )
  261.          IRA = IZAMAX( N-K+1, A( I, K ), LDA )
  262.          RA = ABS( A( I, IRA+K-1 ) )
  263. *
  264. *        Guard against zero C or R due to underflow.
  265. *
  266.          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
  267.      $      GO TO 200
  268.          G = R / SCLFAC
  269.          F = ONE
  270.          S = C + R
  271.   160    CONTINUE
  272.          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
  273.      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
  274.          F = F*SCLFAC
  275.          C = C*SCLFAC
  276.          CA = CA*SCLFAC
  277.          R = R / SCLFAC
  278.          G = G / SCLFAC
  279.          RA = RA / SCLFAC
  280.          GO TO 160
  281. *
  282.   170    CONTINUE
  283.          G = C / SCLFAC
  284.   180    CONTINUE
  285.          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
  286.      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
  287.          F = F / SCLFAC
  288.          C = C / SCLFAC
  289.          G = G / SCLFAC
  290.          CA = CA / SCLFAC
  291.          R = R*SCLFAC
  292.          RA = RA*SCLFAC
  293.          GO TO 180
  294. *
  295. *        Now balance.
  296. *
  297.   190    CONTINUE
  298.          IF( ( C+R ).GE.FACTOR*S )
  299.      $      GO TO 200
  300.          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
  301.             IF( F*SCALE( I ).LE.SFMIN1 )
  302.      $         GO TO 200
  303.          END IF
  304.          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
  305.             IF( SCALE( I ).GE.SFMAX1 / F )
  306.      $         GO TO 200
  307.          END IF
  308.          G = ONE / F
  309.          SCALE( I ) = SCALE( I )*F
  310.          NOCONV = .TRUE.
  311. *
  312.          CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
  313.          CALL ZDSCAL( L, F, A( 1, I ), 1 )
  314. *
  315.   200 CONTINUE
  316. *
  317.       IF( NOCONV )
  318.      $   GO TO 140
  319. *
  320.   210 CONTINUE
  321.       ILO = K
  322.       IHI = L
  323. *
  324.       RETURN
  325. *
  326. *     End of ZGEBAL
  327. *
  328.       END
  329.