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