SRC/dlacon.f

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