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

CVSweb interface <joel.bertrand@systella.fr>