LAPACK: SRC/dlacn2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            KASE, N
00010       DOUBLE PRECISION   EST
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            ISGN( * ), ISAVE( 3 )
00014       DOUBLE PRECISION   V( * ), X( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLACN2 estimates the 1-norm of a square, real matrix A.
00021 *  Reverse communication is used for evaluating matrix-vector products.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  N      (input) INTEGER
00027 *         The order of the matrix.  N >= 1.
00028 *
00029 *  V      (workspace) DOUBLE PRECISION array, dimension (N)
00030 *         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
00031 *         (W is not returned).
00032 *
00033 *  X      (input/output) DOUBLE PRECISION array, dimension (N)
00034 *         On an intermediate return, X should be overwritten by
00035 *               A * X,   if KASE=1,
00036 *               A' * X,  if KASE=2,
00037 *         and DLACN2 must be re-called with all the other parameters
00038 *         unchanged.
00039 *
00040 *  ISGN   (workspace) INTEGER array, dimension (N)
00041 *
00042 *  EST    (input/output) DOUBLE PRECISION
00043 *         On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be
00044 *         unchanged from the previous call to DLACN2.
00045 *         On exit, EST is an estimate (a lower bound) for norm(A). 
00046 *
00047 *  KASE   (input/output) INTEGER
00048 *         On the initial call to DLACN2, KASE should be 0.
00049 *         On an intermediate return, KASE will be 1 or 2, indicating
00050 *         whether X should be overwritten by A * X  or A' * X.
00051 *         On the final return from DLACN2, KASE will again be 0.
00052 *
00053 *  ISAVE  (input/output) INTEGER array, dimension (3)
00054 *         ISAVE is used to save variables between calls to DLACN2
00055 *
00056 *  Further Details
00057 *  ======= =======
00058 *
00059 *  Contributed by Nick Higham, University of Manchester.
00060 *  Originally named SONEST, dated March 16, 1988.
00061 *
00062 *  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
00063 *  a real or complex matrix, with applications to condition estimation",
00064 *  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
00065 *
00066 *  This is a thread safe version of DLACON, which uses the array ISAVE
00067 *  in place of a SAVE statement, as follows:
00068 *
00069 *     DLACON     DLACN2
00070 *      JUMP     ISAVE(1)
00071 *      J        ISAVE(2)
00072 *      ITER     ISAVE(3)
00073 *
00074 *  =====================================================================
00075 *
00076 *     .. Parameters ..
00077       INTEGER            ITMAX
00078       PARAMETER          ( ITMAX = 5 )
00079       DOUBLE PRECISION   ZERO, ONE, TWO
00080       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
00081 *     ..
00082 *     .. Local Scalars ..
00083       INTEGER            I, JLAST
00084       DOUBLE PRECISION   ALTSGN, ESTOLD, TEMP
00085 *     ..
00086 *     .. External Functions ..
00087       INTEGER            IDAMAX
00088       DOUBLE PRECISION   DASUM
00089       EXTERNAL           IDAMAX, DASUM
00090 *     ..
00091 *     .. External Subroutines ..
00092       EXTERNAL           DCOPY
00093 *     ..
00094 *     .. Intrinsic Functions ..
00095       INTRINSIC          ABS, DBLE, NINT, SIGN
00096 *     ..
00097 *     .. Executable Statements ..
00098 *
00099       IF( KASE.EQ.0 ) THEN
00100          DO 10 I = 1, N
00101             X( I ) = ONE / DBLE( N )
00102    10    CONTINUE
00103          KASE = 1
00104          ISAVE( 1 ) = 1
00105          RETURN
00106       END IF
00107 *
00108       GO TO ( 20, 40, 70, 110, 140 )ISAVE( 1 )
00109 *
00110 *     ................ ENTRY   (ISAVE( 1 ) = 1)
00111 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
00112 *
00113    20 CONTINUE
00114       IF( N.EQ.1 ) THEN
00115          V( 1 ) = X( 1 )
00116          EST = ABS( V( 1 ) )
00117 *        ... QUIT
00118          GO TO 150
00119       END IF
00120       EST = DASUM( N, X, 1 )
00121 *
00122       DO 30 I = 1, N
00123          X( I ) = SIGN( ONE, X( I ) )
00124          ISGN( I ) = NINT( X( I ) )
00125    30 CONTINUE
00126       KASE = 2
00127       ISAVE( 1 ) = 2
00128       RETURN
00129 *
00130 *     ................ ENTRY   (ISAVE( 1 ) = 2)
00131 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
00132 *
00133    40 CONTINUE
00134       ISAVE( 2 ) = IDAMAX( N, X, 1 )
00135       ISAVE( 3 ) = 2
00136 *
00137 *     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
00138 *
00139    50 CONTINUE
00140       DO 60 I = 1, N
00141          X( I ) = ZERO
00142    60 CONTINUE
00143       X( ISAVE( 2 ) ) = ONE
00144       KASE = 1
00145       ISAVE( 1 ) = 3
00146       RETURN
00147 *
00148 *     ................ ENTRY   (ISAVE( 1 ) = 3)
00149 *     X HAS BEEN OVERWRITTEN BY A*X.
00150 *
00151    70 CONTINUE
00152       CALL DCOPY( N, X, 1, V, 1 )
00153       ESTOLD = EST
00154       EST = DASUM( N, V, 1 )
00155       DO 80 I = 1, N
00156          IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) )
00157      $      GO TO 90
00158    80 CONTINUE
00159 *     REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
00160       GO TO 120
00161 *
00162    90 CONTINUE
00163 *     TEST FOR CYCLING.
00164       IF( EST.LE.ESTOLD )
00165      $   GO TO 120
00166 *
00167       DO 100 I = 1, N
00168          X( I ) = SIGN( ONE, X( I ) )
00169          ISGN( I ) = NINT( X( I ) )
00170   100 CONTINUE
00171       KASE = 2
00172       ISAVE( 1 ) = 4
00173       RETURN
00174 *
00175 *     ................ ENTRY   (ISAVE( 1 ) = 4)
00176 *     X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
00177 *
00178   110 CONTINUE
00179       JLAST = ISAVE( 2 )
00180       ISAVE( 2 ) = IDAMAX( N, X, 1 )
00181       IF( ( X( JLAST ).NE.ABS( X( ISAVE( 2 ) ) ) ) .AND.
00182      $    ( ISAVE( 3 ).LT.ITMAX ) ) THEN
00183          ISAVE( 3 ) = ISAVE( 3 ) + 1
00184          GO TO 50
00185       END IF
00186 *
00187 *     ITERATION COMPLETE.  FINAL STAGE.
00188 *
00189   120 CONTINUE
00190       ALTSGN = ONE
00191       DO 130 I = 1, N
00192          X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) )
00193          ALTSGN = -ALTSGN
00194   130 CONTINUE
00195       KASE = 1
00196       ISAVE( 1 ) = 5
00197       RETURN
00198 *
00199 *     ................ ENTRY   (ISAVE( 1 ) = 5)
00200 *     X HAS BEEN OVERWRITTEN BY A*X.
00201 *
00202   140 CONTINUE
00203       TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) )
00204       IF( TEMP.GT.EST ) THEN
00205          CALL DCOPY( N, X, 1, V, 1 )
00206          EST = TEMP
00207       END IF
00208 *
00209   150 CONTINUE
00210       KASE = 0
00211       RETURN
00212 *
00213 *     End of DLACN2
00214 *
00215       END