SRC/slacn2.f

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