File:  [local] / rpl / lapack / lapack / dlacn2.f
Revision 1.12: download - view: text, annotated - select for diffs - revision graph
Fri Dec 14 12:30:21 2012 UTC (11 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Mise à jour de Lapack vers la version 3.4.2 et des scripts de compilation
pour rplcas. En particulier, le Makefile.am de giac a été modifié pour ne
compiler que le répertoire src.

    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: *> \date September 2012
  105: *
  106: *> \ingroup doubleOTHERauxiliary
  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 DLACON, which uses the array ISAVE
  116: *>  in place of a SAVE statement, as follows:
  117: *>
  118: *>     DLACON     DLACN2
  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 DLACN2( 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:       DOUBLE PRECISION   EST
  147: *     ..
  148: *     .. Array Arguments ..
  149:       INTEGER            ISGN( * ), ISAVE( 3 )
  150:       DOUBLE PRECISION   V( * ), X( * )
  151: *     ..
  152: *
  153: *  =====================================================================
  154: *
  155: *     .. Parameters ..
  156:       INTEGER            ITMAX
  157:       PARAMETER          ( ITMAX = 5 )
  158:       DOUBLE PRECISION   ZERO, ONE, TWO
  159:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
  160: *     ..
  161: *     .. Local Scalars ..
  162:       INTEGER            I, JLAST
  163:       DOUBLE PRECISION   ALTSGN, ESTOLD, TEMP
  164: *     ..
  165: *     .. External Functions ..
  166:       INTEGER            IDAMAX
  167:       DOUBLE PRECISION   DASUM
  168:       EXTERNAL           IDAMAX, DASUM
  169: *     ..
  170: *     .. External Subroutines ..
  171:       EXTERNAL           DCOPY
  172: *     ..
  173: *     .. Intrinsic Functions ..
  174:       INTRINSIC          ABS, DBLE, NINT, SIGN
  175: *     ..
  176: *     .. Executable Statements ..
  177: *
  178:       IF( KASE.EQ.0 ) THEN
  179:          DO 10 I = 1, N
  180:             X( I ) = ONE / DBLE( 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 = DASUM( 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 ) = IDAMAX( 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 DCOPY( N, X, 1, V, 1 )
  232:       ESTOLD = EST
  233:       EST = DASUM( 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 ) = IDAMAX( 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+DBLE( I-1 ) / DBLE( 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*( DASUM( N, X, 1 ) / DBLE( 3*N ) )
  283:       IF( TEMP.GT.EST ) THEN
  284:          CALL DCOPY( N, X, 1, V, 1 )
  285:          EST = TEMP
  286:       END IF
  287: *
  288:   150 CONTINUE
  289:       KASE = 0
  290:       RETURN
  291: *
  292: *     End of DLACN2
  293: *
  294:       END

CVSweb interface <joel.bertrand@systella.fr>