dlacn2.f

Go to the documentation of this file.

      SUBROUTINE dlacn2( N, V, X, ISGN, EST, KASE, ISAVE )
*
*  -- LAPACK auxiliary routine (version 3.2) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            kase, n
      DOUBLE PRECISION   est
*     ..
*     .. Array Arguments ..
      INTEGER            isgn( * ), isave( 3 )
      DOUBLE PRECISION   v( * ), x( * )
*     ..
*
*  Purpose
*  =======
*
*  DLACN2 estimates the 1-norm of a square, real matrix A.
*  Reverse communication is used for evaluating matrix-vector products.
*
*  Arguments
*  =========
*
*  N      (input) INTEGER
*         The order of the matrix.  N >= 1.
*
*  V      (workspace) DOUBLE PRECISION array, dimension (N)
*         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
*         (W is not returned).
*
*  X      (input/output) DOUBLE PRECISION array, dimension (N)
*         On an intermediate return, X should be overwritten by
*               A * X,   if KASE=1,
*               A' * X,  if KASE=2,
*         and DLACN2 must be re-called with all the other parameters
*         unchanged.
*
*  ISGN   (workspace) INTEGER array, dimension (N)
*
*  EST    (input/output) DOUBLE PRECISION
*         On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be
*         unchanged from the previous call to DLACN2.
*         On exit, EST is an estimate (a lower bound) for norm(A). 
*
*  KASE   (input/output) INTEGER
*         On the initial call to DLACN2, KASE should be 0.
*         On an intermediate return, KASE will be 1 or 2, indicating
*         whether X should be overwritten by A * X  or A' * X.
*         On the final return from DLACN2, KASE will again be 0.
*
*  ISAVE  (input/output) INTEGER array, dimension (3)
*         ISAVE is used to save variables between calls to DLACN2
*
*  Further Details
*  ======= =======
*
*  Contributed by Nick Higham, University of Manchester.
*  Originally named SONEST, dated March 16, 1988.
*
*  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
*  a real or complex matrix, with applications to condition estimation",
*  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
*
*  This is a thread safe version of DLACON, which uses the array ISAVE
*  in place of a SAVE statement, as follows:
*
*     DLACON     DLACN2
*      JUMP     ISAVE(1)
*      J        ISAVE(2)
*      ITER     ISAVE(3)
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            itmax
      parameter( itmax = 5 )
      DOUBLE PRECISION   zero, one, two
      parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            i, jlast
      DOUBLE PRECISION   altsgn, estold, temp
*     ..
*     .. External Functions ..
      INTEGER            idamax
      DOUBLE PRECISION   dasum
      EXTERNAL           idamax, dasum
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, nint, sign
*     ..
*     .. Executable Statements ..
*
      IF( kase.EQ.0 ) THEN
         DO 10 i = 1, n
            x( i ) = one / dble( n )
   10    continue
         kase = 1
         isave( 1 ) = 1
         return
      END IF
*
      go to( 20, 40, 70, 110, 140 )isave( 1 )
*
*     ................ ENTRY   (ISAVE( 1 ) = 1)
*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
*
   20 continue
      IF( n.EQ.1 ) THEN
         v( 1 ) = x( 1 )
         est = abs( v( 1 ) )
*        ... QUIT
         go to 150
      END IF
      est = dasum( n, x, 1 )
*
      DO 30 i = 1, n
         x( i ) = sign( one, x( i ) )
         isgn( i ) = nint( x( i ) )
   30 continue
      kase = 2
      isave( 1 ) = 2
      return
*
*     ................ ENTRY   (ISAVE( 1 ) = 2)
*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
*
   40 continue
      isave( 2 ) = idamax( n, x, 1 )
      isave( 3 ) = 2
*
*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
*
   50 continue
      DO 60 i = 1, n
         x( i ) = zero
   60 continue
      x( isave( 2 ) ) = one
      kase = 1
      isave( 1 ) = 3
      return
*
*     ................ ENTRY   (ISAVE( 1 ) = 3)
*     X HAS BEEN OVERWRITTEN BY A*X.
*
   70 continue
      CALL dcopy( n, x, 1, v, 1 )
      estold = est
      est = dasum( n, v, 1 )
      DO 80 i = 1, n
         IF( nint( sign( one, x( i ) ) ).NE.isgn( i ) )
     $      go to 90
   80 continue
*     REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
      go to 120
*
   90 continue
*     TEST FOR CYCLING.
      IF( est.LE.estold )
     $   go to 120
*
      DO 100 i = 1, n
         x( i ) = sign( one, x( i ) )
         isgn( i ) = nint( x( i ) )
  100 continue
      kase = 2
      isave( 1 ) = 4
      return
*
*     ................ ENTRY   (ISAVE( 1 ) = 4)
*     X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
*
  110 continue
      jlast = isave( 2 )
      isave( 2 ) = idamax( n, x, 1 )
      IF( ( x( jlast ).NE.abs( x( isave( 2 ) ) ) ) .AND.
     $    ( isave( 3 ).LT.itmax ) ) THEN
         isave( 3 ) = isave( 3 ) + 1
         go to 50
      END IF
*
*     ITERATION COMPLETE.  FINAL STAGE.
*
  120 continue
      altsgn = one
      DO 130 i = 1, n
         x( i ) = altsgn*( one+dble( i-1 ) / dble( n-1 ) )
         altsgn = -altsgn
  130 continue
      kase = 1
      isave( 1 ) = 5
      return
*
*     ................ ENTRY   (ISAVE( 1 ) = 5)
*     X HAS BEEN OVERWRITTEN BY A*X.
*
  140 continue
      temp = two*( dasum( n, x, 1 ) / dble( 3*n ) )
      IF( temp.GT.est ) THEN
         CALL dcopy( n, x, 1, v, 1 )
         est = temp
      END IF
*
  150 continue
      kase = 0
      return
*
*     End of DLACN2
*
      END