File:  [local] / rpl / lapack / lapack / dlacon.f
Revision 1.19: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:53 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

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

CVSweb interface <joel.bertrand@systella.fr>