SRC/zlacon.f

   1 *> \brief \b ZLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLACON + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlacon.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlacon.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlacon.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLACON( N, V, X, EST, KASE )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            KASE, N
  25 *       DOUBLE PRECISION   EST
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX*16         V( N ), X( N )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZLACON estimates the 1-norm of a square, complex matrix A.
  38 *> Reverse communication is used for evaluating matrix-vector products.
  39 *> \endverbatim
  40 *
  41 *  Arguments:
  42 *  ==========
  43 *
  44 *> \param[in] N
  45 *> \verbatim
  46 *>          N is INTEGER
  47 *>         The order of the matrix.  N >= 1.
  48 *> \endverbatim
  49 *>
  50 *> \param[out] V
  51 *> \verbatim
  52 *>          V is COMPLEX*16 array, dimension (N)
  53 *>         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
  54 *>         (W is not returned).
  55 *> \endverbatim
  56 *>
  57 *> \param[in,out] X
  58 *> \verbatim
  59 *>          X is COMPLEX*16 array, dimension (N)
  60 *>         On an intermediate return, X should be overwritten by
  61 *>               A * X,   if KASE=1,
  62 *>               A**H * X,  if KASE=2,
  63 *>         where A**H is the conjugate transpose of A, and ZLACON must be
  64 *>         re-called with all the other parameters unchanged.
  65 *> \endverbatim
  66 *>
  67 *> \param[in,out] EST
  68 *> \verbatim
  69 *>          EST is DOUBLE PRECISION
  70 *>         On entry with KASE = 1 or 2 and JUMP = 3, EST should be
  71 *>         unchanged from the previous call to ZLACON.
  72 *>         On exit, EST is an estimate (a lower bound) for norm(A).
  73 *> \endverbatim
  74 *>
  75 *> \param[in,out] KASE
  76 *> \verbatim
  77 *>          KASE is INTEGER
  78 *>         On the initial call to ZLACON, KASE should be 0.
  79 *>         On an intermediate return, KASE will be 1 or 2, indicating
  80 *>         whether X should be overwritten by A * X  or A**H * X.
  81 *>         On the final return from ZLACON, KASE will again be 0.
  82 *> \endverbatim
  83 *
  84 *  Authors:
  85 *  ========
  86 *
  87 *> \author Univ. of Tennessee
  88 *> \author Univ. of California Berkeley
  89 *> \author Univ. of Colorado Denver
  90 *> \author NAG Ltd.
  91 *
  92 *> \date September 2012
  93 *
  94 *> \ingroup complex16OTHERauxiliary
  95 *
  96 *> \par Further Details:
  97 *  =====================
  98 *>
  99 *>  Originally named CONEST, dated March 16, 1988. \n
 100 *>  Last modified:  April, 1999
 101 *
 102 *> \par Contributors:
 103 *  ==================
 104 *>
 105 *>     Nick Higham, University of Manchester
 106 *
 107 *> \par References:
 108 *  ================
 109 *>
 110 *>  N.J. Higham, "FORTRAN codes for estimating the one-norm of
 111 *>  a real or complex matrix, with applications to condition estimation",
 112 *>  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
 113 *>
 114 *  =====================================================================
 115       SUBROUTINE ZLACON( N, V, X, EST, KASE )
 116 *
 117 *  -- LAPACK auxiliary routine (version 3.4.2) --
 118 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 119 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 120 *     September 2012
 121 *
 122 *     .. Scalar Arguments ..
 123       INTEGER            KASE, N
 124       DOUBLE PRECISION   EST
 125 *     ..
 126 *     .. Array Arguments ..
 127       COMPLEX*16         V( N ), X( N )
 128 *     ..
 129 *
 130 *  =====================================================================
 131 *
 132 *     .. Parameters ..
 133       INTEGER            ITMAX
 134       PARAMETER          ( ITMAX = 5 )
 135       DOUBLE PRECISION   ONE, TWO
 136       PARAMETER          ( ONE = 1.0D0, TWO = 2.0D0 )
 137       COMPLEX*16         CZERO, CONE
 138       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
 139      $                   CONE = ( 1.0D0, 0.0D0 ) )
 140 *     ..
 141 *     .. Local Scalars ..
 142       INTEGER            I, ITER, J, JLAST, JUMP
 143       DOUBLE PRECISION   ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
 144 *     ..
 145 *     .. External Functions ..
 146       INTEGER            IZMAX1
 147       DOUBLE PRECISION   DLAMCH, DZSUM1
 148       EXTERNAL           IZMAX1, DLAMCH, DZSUM1
 149 *     ..
 150 *     .. External Subroutines ..
 151       EXTERNAL           ZCOPY
 152 *     ..
 153 *     .. Intrinsic Functions ..
 154       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG
 155 *     ..
 156 *     .. Save statement ..
 157       SAVE
 158 *     ..
 159 *     .. Executable Statements ..
 160 *
 161       SAFMIN = DLAMCH( 'Safe minimum' )
 162       IF( KASE.EQ.0 ) THEN
 163          DO 10 I = 1, N
 164             X( I ) = DCMPLX( ONE / DBLE( N ) )
 165    10    CONTINUE
 166          KASE = 1
 167          JUMP = 1
 168          RETURN
 169       END IF
 170 *
 171       GO TO ( 20, 40, 70, 90, 120 )JUMP
 172 *
 173 *     ................ ENTRY   (JUMP = 1)
 174 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
 175 *
 176    20 CONTINUE
 177       IF( N.EQ.1 ) THEN
 178          V( 1 ) = X( 1 )
 179          EST = ABS( V( 1 ) )
 180 *        ... QUIT
 181          GO TO 130
 182       END IF
 183       EST = DZSUM1( N, X, 1 )
 184 *
 185       DO 30 I = 1, N
 186          ABSXI = ABS( X( I ) )
 187          IF( ABSXI.GT.SAFMIN ) THEN
 188             X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
 189      $               DIMAG( X( I ) ) / ABSXI )
 190          ELSE
 191             X( I ) = CONE
 192          END IF
 193    30 CONTINUE
 194       KASE = 2
 195       JUMP = 2
 196       RETURN
 197 *
 198 *     ................ ENTRY   (JUMP = 2)
 199 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
 200 *
 201    40 CONTINUE
 202       J = IZMAX1( N, X, 1 )
 203       ITER = 2
 204 *
 205 *     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
 206 *
 207    50 CONTINUE
 208       DO 60 I = 1, N
 209          X( I ) = CZERO
 210    60 CONTINUE
 211       X( J ) = CONE
 212       KASE = 1
 213       JUMP = 3
 214       RETURN
 215 *
 216 *     ................ ENTRY   (JUMP = 3)
 217 *     X HAS BEEN OVERWRITTEN BY A*X.
 218 *
 219    70 CONTINUE
 220       CALL ZCOPY( N, X, 1, V, 1 )
 221       ESTOLD = EST
 222       EST = DZSUM1( N, V, 1 )
 223 *
 224 *     TEST FOR CYCLING.
 225       IF( EST.LE.ESTOLD )
 226      $   GO TO 100
 227 *
 228       DO 80 I = 1, N
 229          ABSXI = ABS( X( I ) )
 230          IF( ABSXI.GT.SAFMIN ) THEN
 231             X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
 232      $               DIMAG( X( I ) ) / ABSXI )
 233          ELSE
 234             X( I ) = CONE
 235          END IF
 236    80 CONTINUE
 237       KASE = 2
 238       JUMP = 4
 239       RETURN
 240 *
 241 *     ................ ENTRY   (JUMP = 4)
 242 *     X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
 243 *
 244    90 CONTINUE
 245       JLAST = J
 246       J = IZMAX1( N, X, 1 )
 247       IF( ( ABS( X( JLAST ) ).NE.ABS( X( J ) ) ) .AND.
 248      $    ( ITER.LT.ITMAX ) ) THEN
 249          ITER = ITER + 1
 250          GO TO 50
 251       END IF
 252 *
 253 *     ITERATION COMPLETE.  FINAL STAGE.
 254 *
 255   100 CONTINUE
 256       ALTSGN = ONE
 257       DO 110 I = 1, N
 258          X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
 259          ALTSGN = -ALTSGN
 260   110 CONTINUE
 261       KASE = 1
 262       JUMP = 5
 263       RETURN
 264 *
 265 *     ................ ENTRY   (JUMP = 5)
 266 *     X HAS BEEN OVERWRITTEN BY A*X.
 267 *
 268   120 CONTINUE
 269       TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
 270       IF( TEMP.GT.EST ) THEN
 271          CALL ZCOPY( N, X, 1, V, 1 )
 272          EST = TEMP
 273       END IF
 274 *
 275   130 CONTINUE
 276       KASE = 0
 277       RETURN
 278 *
 279 *     End of ZLACON
 280 *
 281       END