LAPACK: SRC/slacon.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLACON( N, V, X, ISGN, EST, KASE )
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       REAL               EST
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            ISGN( * )
00014       REAL               V( * ), X( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SLACON 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) REAL 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) REAL 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 SLACON must be re-called with all the other parameters
00038 *         unchanged.
00039 *
00040 *  ISGN   (workspace) INTEGER array, dimension (N)
00041 *
00042 *  EST    (input/output) REAL
00043 *         On entry with KASE = 1 or 2 and JUMP = 3, EST should be
00044 *         unchanged from the previous call to SLACON.
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 SLACON, 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 SLACON, KASE will again be 0.
00052 *
00053 *  Further Details
00054 *  ======= =======
00055 *
00056 *  Contributed by Nick Higham, University of Manchester.
00057 *  Originally named SONEST, dated March 16, 1988.
00058 *
00059 *  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
00060 *  a real or complex matrix, with applications to condition estimation",
00061 *  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
00062 *
00063 *  =====================================================================
00064 *
00065 *     .. Parameters ..
00066       INTEGER            ITMAX
00067       PARAMETER          ( ITMAX = 5 )
00068       REAL               ZERO, ONE, TWO
00069       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
00070 *     ..
00071 *     .. Local Scalars ..
00072       INTEGER            I, ITER, J, JLAST, JUMP
00073       REAL               ALTSGN, ESTOLD, TEMP
00074 *     ..
00075 *     .. External Functions ..
00076       INTEGER            ISAMAX
00077       REAL               SASUM
00078       EXTERNAL           ISAMAX, SASUM
00079 *     ..
00080 *     .. External Subroutines ..
00081       EXTERNAL           SCOPY
00082 *     ..
00083 *     .. Intrinsic Functions ..
00084       INTRINSIC          ABS, NINT, REAL, SIGN
00085 *     ..
00086 *     .. Save statement ..
00087       SAVE
00088 *     ..
00089 *     .. Executable Statements ..
00090 *
00091       IF( KASE.EQ.0 ) THEN
00092          DO 10 I = 1, N
00093             X( I ) = ONE / REAL( N )
00094    10    CONTINUE
00095          KASE = 1
00096          JUMP = 1
00097          RETURN
00098       END IF
00099 *
00100       GO TO ( 20, 40, 70, 110, 140 )JUMP
00101 *
00102 *     ................ ENTRY   (JUMP = 1)
00103 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
00104 *
00105    20 CONTINUE
00106       IF( N.EQ.1 ) THEN
00107          V( 1 ) = X( 1 )
00108          EST = ABS( V( 1 ) )
00109 *        ... QUIT
00110          GO TO 150
00111       END IF
00112       EST = SASUM( N, X, 1 )
00113 *
00114       DO 30 I = 1, N
00115          X( I ) = SIGN( ONE, X( I ) )
00116          ISGN( I ) = NINT( X( I ) )
00117    30 CONTINUE
00118       KASE = 2
00119       JUMP = 2
00120       RETURN
00121 *
00122 *     ................ ENTRY   (JUMP = 2)
00123 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
00124 *
00125    40 CONTINUE
00126       J = ISAMAX( N, X, 1 )
00127       ITER = 2
00128 *
00129 *     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
00130 *
00131    50 CONTINUE
00132       DO 60 I = 1, N
00133          X( I ) = ZERO
00134    60 CONTINUE
00135       X( J ) = ONE
00136       KASE = 1
00137       JUMP = 3
00138       RETURN
00139 *
00140 *     ................ ENTRY   (JUMP = 3)
00141 *     X HAS BEEN OVERWRITTEN BY A*X.
00142 *
00143    70 CONTINUE
00144       CALL SCOPY( N, X, 1, V, 1 )
00145       ESTOLD = EST
00146       EST = SASUM( N, V, 1 )
00147       DO 80 I = 1, N
00148          IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) )
00149      $      GO TO 90
00150    80 CONTINUE
00151 *     REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
00152       GO TO 120
00153 *
00154    90 CONTINUE
00155 *     TEST FOR CYCLING.
00156       IF( EST.LE.ESTOLD )
00157      $   GO TO 120
00158 *
00159       DO 100 I = 1, N
00160          X( I ) = SIGN( ONE, X( I ) )
00161          ISGN( I ) = NINT( X( I ) )
00162   100 CONTINUE
00163       KASE = 2
00164       JUMP = 4
00165       RETURN
00166 *
00167 *     ................ ENTRY   (JUMP = 4)
00168 *     X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
00169 *
00170   110 CONTINUE
00171       JLAST = J
00172       J = ISAMAX( N, X, 1 )
00173       IF( ( X( JLAST ).NE.ABS( X( J ) ) ) .AND. ( ITER.LT.ITMAX ) ) THEN
00174          ITER = ITER + 1
00175          GO TO 50
00176       END IF
00177 *
00178 *     ITERATION COMPLETE.  FINAL STAGE.
00179 *
00180   120 CONTINUE
00181       ALTSGN = ONE
00182       DO 130 I = 1, N
00183          X( I ) = ALTSGN*( ONE+REAL( I-1 ) / REAL( N-1 ) )
00184          ALTSGN = -ALTSGN
00185   130 CONTINUE
00186       KASE = 1
00187       JUMP = 5
00188       RETURN
00189 *
00190 *     ................ ENTRY   (JUMP = 5)
00191 *     X HAS BEEN OVERWRITTEN BY A*X.
00192 *
00193   140 CONTINUE
00194       TEMP = TWO*( SASUM( N, X, 1 ) / REAL( 3*N ) )
00195       IF( TEMP.GT.EST ) THEN
00196          CALL SCOPY( N, X, 1, V, 1 )
00197          EST = TEMP
00198       END IF
00199 *
00200   150 CONTINUE
00201       KASE = 0
00202       RETURN
00203 *
00204 *     End of SLACON
00205 *
00206       END