SRC/clacn2.f

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