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 / zlacon.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  206 lines
  1.       SUBROUTINE ZLACON( N, V, X, EST, KASE )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            KASE, N
  10.       DOUBLE PRECISION   EST
  11. *     ..
  12. *     .. Array Arguments ..
  13.       COMPLEX*16         V( N ), X( N )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  ZLACON estimates the 1-norm of a square, complex matrix A.
  20. *  Reverse communication is used for evaluating matrix-vector products.
  21. *
  22. *  Arguments
  23. *  =========
  24. *
  25. *  N      (input) INTEGER
  26. *         The order of the matrix.  N >= 1.
  27. *
  28. *  V      (workspace) COMPLEX*16 array, dimension (N)
  29. *         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
  30. *         (W is not returned).
  31. *
  32. *  X      (input/output) COMPLEX*16 array, dimension (N)
  33. *         On an intermediate return, X should be overwritten by
  34. *               A * X,   if KASE=1,
  35. *               A' * X,  if KASE=2,
  36. *         where A' is the conjugate transpose of A, and ZLACON must be
  37. *         re-called with all the other parameters unchanged.
  38. *
  39. *  EST    (output) DOUBLE PRECISION
  40. *         An estimate (a lower bound) for norm(A).
  41. *
  42. *  KASE   (input/output) INTEGER
  43. *         On the initial call to ZLACON, KASE should be 0.
  44. *         On an intermediate return, KASE will be 1 or 2, indicating
  45. *         whether X should be overwritten by A * X  or A' * X.
  46. *         On the final return from ZLACON, KASE will again be 0.
  47. *
  48. *  Further Details
  49. *  ======= =======
  50. *
  51. *  Contributed by Nick Higham, University of Manchester.
  52. *  Originally named CONEST, dated March 16, 1988.
  53. *
  54. *  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
  55. *  a real or complex matrix, with applications to condition estimation",
  56. *  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
  57. *
  58. *  =====================================================================
  59. *
  60. *     .. Parameters ..
  61.       INTEGER            ITMAX
  62.       PARAMETER          ( ITMAX = 5 )
  63.       DOUBLE PRECISION   ONE, TWO
  64.       PARAMETER          ( ONE = 1.0D0, TWO = 2.0D0 )
  65.       COMPLEX*16         CZERO, CONE
  66.       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
  67.      $                   CONE = ( 1.0D0, 0.0D0 ) )
  68. *     ..
  69. *     .. Local Scalars ..
  70.       INTEGER            I, ITER, J, JLAST, JUMP
  71.       DOUBLE PRECISION   ALTSGN, ESTOLD, SAFMIN, TEMP
  72. *     ..
  73. *     .. External Functions ..
  74.       INTEGER            IZMAX1
  75.       DOUBLE PRECISION   DLAMCH, DZSUM1
  76.       EXTERNAL           IZMAX1, DLAMCH, DZSUM1
  77. *     ..
  78. *     .. External Subroutines ..
  79.       EXTERNAL           ZCOPY
  80. *     ..
  81. *     .. Intrinsic Functions ..
  82.       INTRINSIC          ABS, DBLE, DCMPLX
  83. *     ..
  84. *     .. Save statement ..
  85.       SAVE
  86. *     ..
  87. *     .. Executable Statements ..
  88. *
  89.       SAFMIN = DLAMCH( 'Safe minimum' )
  90.       IF( KASE.EQ.0 ) THEN
  91.          DO 10 I = 1, N
  92.             X( I ) = DCMPLX( ONE / DBLE( N ) )
  93.    10    CONTINUE
  94.          KASE = 1
  95.          JUMP = 1
  96.          RETURN
  97.       END IF
  98. *
  99.       GO TO ( 20, 40, 70, 90, 120 )JUMP
  100. *
  101. *     ................ ENTRY   (JUMP = 1)
  102. *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
  103. *
  104.    20 CONTINUE
  105.       IF( N.EQ.1 ) THEN
  106.          V( 1 ) = X( 1 )
  107.          EST = ABS( V( 1 ) )
  108. *        ... QUIT
  109.          GO TO 130
  110.       END IF
  111.       EST = DZSUM1( N, X, 1 )
  112. *
  113.       DO 30 I = 1, N
  114.          IF( ABS( X( I ) ).GT.SAFMIN ) THEN
  115.             X( I ) = X( I ) / DCMPLX( ABS( X( I ) ) )
  116.          ELSE
  117.             X( I ) = CONE
  118.          END IF
  119.    30 CONTINUE
  120.       KASE = 2
  121.       JUMP = 2
  122.       RETURN
  123. *
  124. *     ................ ENTRY   (JUMP = 2)
  125. *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY ZTRANS(A)*X.
  126. *
  127.    40 CONTINUE
  128.       J = IZMAX1( N, X, 1 )
  129.       ITER = 2
  130. *
  131. *     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
  132. *
  133.    50 CONTINUE
  134.       DO 60 I = 1, N
  135.          X( I ) = CZERO
  136.    60 CONTINUE
  137.       X( J ) = CONE
  138.       KASE = 1
  139.       JUMP = 3
  140.       RETURN
  141. *
  142. *     ................ ENTRY   (JUMP = 3)
  143. *     X HAS BEEN OVERWRITTEN BY A*X.
  144. *
  145.    70 CONTINUE
  146.       CALL ZCOPY( N, X, 1, V, 1 )
  147.       ESTOLD = EST
  148.       EST = DZSUM1( N, V, 1 )
  149. *
  150. *     TEST FOR CYCLING.
  151.       IF( EST.LE.ESTOLD )
  152.      $   GO TO 100
  153. *
  154.       DO 80 I = 1, N
  155.          IF( ABS( X( I ) ).GT.SAFMIN ) THEN
  156.             X( I ) = X( I ) / DCMPLX( ABS( X( I ) ) )
  157.          ELSE
  158.             X( I ) = CONE
  159.          END IF
  160.    80 CONTINUE
  161.       KASE = 2
  162.       JUMP = 4
  163.       RETURN
  164. *
  165. *     ................ ENTRY   (JUMP = 4)
  166. *     X HAS BEEN OVERWRITTEN BY ZTRANS(A)*X.
  167. *
  168.    90 CONTINUE
  169.       JLAST = J
  170.       J = IZMAX1( N, X, 1 )
  171.       IF( ( DBLE( X( JLAST ) ).NE.ABS( DBLE( X( J ) ) ) ) .AND.
  172.      $    ( ITER.LT.ITMAX ) ) THEN
  173.          ITER = ITER + 1
  174.          GO TO 50
  175.       END IF
  176. *
  177. *     ITERATION COMPLETE.  FINAL STAGE.
  178. *
  179.   100 CONTINUE
  180.       ALTSGN = ONE
  181.       DO 110 I = 1, N
  182.          X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
  183.          ALTSGN = -ALTSGN
  184.   110 CONTINUE
  185.       KASE = 1
  186.       JUMP = 5
  187.       RETURN
  188. *
  189. *     ................ ENTRY   (JUMP = 5)
  190. *     X HAS BEEN OVERWRITTEN BY A*X.
  191. *
  192.   120 CONTINUE
  193.       TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
  194.       IF( TEMP.GT.EST ) THEN
  195.          CALL ZCOPY( N, X, 1, V, 1 )
  196.          EST = TEMP
  197.       END IF
  198. *
  199.   130 CONTINUE
  200.       KASE = 0
  201.       RETURN
  202. *
  203. *     End of ZLACON
  204. *
  205.       END
  206.