LAPACK: SRC/clacn2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLACN2( N, V, X, 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       REAL               EST
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            ISAVE( 3 )
00014       COMPLEX            V( * ), X( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CLACN2 estimates the 1-norm of a square, complex 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) COMPLEX 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) COMPLEX 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 *         where A' is the conjugate transpose of A, and CLACN2 must be
00038 *         re-called with all the other parameters unchanged.
00039 *
00040 *  EST    (input/output) REAL
00041 *         On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be
00042 *         unchanged from the previous call to CLACN2.
00043 *         On exit, EST is an estimate (a lower bound) for norm(A). 
00044 *
00045 *  KASE   (input/output) INTEGER
00046 *         On the initial call to CLACN2, KASE should be 0.
00047 *         On an intermediate return, KASE will be 1 or 2, indicating
00048 *         whether X should be overwritten by A * X  or A' * X.
00049 *         On the final return from CLACN2, KASE will again be 0.
00050 *
00051 *  ISAVE  (input/output) INTEGER array, dimension (3)
00052 *         ISAVE is used to save variables between calls to SLACN2
00053 *
00054 *  Further Details
00055 *  ======= =======
00056 *
00057 *  Contributed by Nick Higham, University of Manchester.
00058 *  Originally named CONEST, dated March 16, 1988.
00059 *
00060 *  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
00061 *  a real or complex matrix, with applications to condition estimation",
00062 *  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
00063 *
00064 *  Last modified:  April, 1999
00065 *
00066 *  This is a thread safe version of CLACON, which uses the array ISAVE
00067 *  in place of a SAVE statement, as follows:
00068 *
00069 *     CLACON     CLACN2
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       REAL                 ONE,         TWO
00080       PARAMETER          ( ONE = 1.0E0, TWO = 2.0E0 )
00081       COMPLEX              CZERO, CONE
00082       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
00083      $                            CONE = ( 1.0E0, 0.0E0 ) )
00084 *     ..
00085 *     .. Local Scalars ..
00086       INTEGER            I, JLAST
00087       REAL               ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
00088 *     ..
00089 *     .. External Functions ..
00090       INTEGER            ICMAX1
00091       REAL               SCSUM1, SLAMCH
00092       EXTERNAL           ICMAX1, SCSUM1, SLAMCH
00093 *     ..
00094 *     .. External Subroutines ..
00095       EXTERNAL           CCOPY
00096 *     ..
00097 *     .. Intrinsic Functions ..
00098       INTRINSIC          ABS, AIMAG, CMPLX, REAL
00099 *     ..
00100 *     .. Executable Statements ..
00101 *
00102       SAFMIN = SLAMCH( 'Safe minimum' )
00103       IF( KASE.EQ.0 ) THEN
00104          DO 10 I = 1, N
00105             X( I ) = CMPLX( ONE / REAL( N ) )
00106    10    CONTINUE
00107          KASE = 1
00108          ISAVE( 1 ) = 1
00109          RETURN
00110       END IF
00111 *
00112       GO TO ( 20, 40, 70, 90, 120 )ISAVE( 1 )
00113 *
00114 *     ................ ENTRY   (ISAVE( 1 ) = 1)
00115 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
00116 *
00117    20 CONTINUE
00118       IF( N.EQ.1 ) THEN
00119          V( 1 ) = X( 1 )
00120          EST = ABS( V( 1 ) )
00121 *        ... QUIT
00122          GO TO 130
00123       END IF
00124       EST = SCSUM1( N, X, 1 )
00125 *
00126       DO 30 I = 1, N
00127          ABSXI = ABS( X( I ) )
00128          IF( ABSXI.GT.SAFMIN ) THEN
00129             X( I ) = CMPLX( REAL( X( I ) ) / ABSXI,
00130      $               AIMAG( X( I ) ) / ABSXI )
00131          ELSE
00132             X( I ) = CONE
00133          END IF
00134    30 CONTINUE
00135       KASE = 2
00136       ISAVE( 1 ) = 2
00137       RETURN
00138 *
00139 *     ................ ENTRY   (ISAVE( 1 ) = 2)
00140 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
00141 *
00142    40 CONTINUE
00143       ISAVE( 2 ) = ICMAX1( N, X, 1 )
00144       ISAVE( 3 ) = 2
00145 *
00146 *     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
00147 *
00148    50 CONTINUE
00149       DO 60 I = 1, N
00150          X( I ) = CZERO
00151    60 CONTINUE
00152       X( ISAVE( 2 ) ) = CONE
00153       KASE = 1
00154       ISAVE( 1 ) = 3
00155       RETURN
00156 *
00157 *     ................ ENTRY   (ISAVE( 1 ) = 3)
00158 *     X HAS BEEN OVERWRITTEN BY A*X.
00159 *
00160    70 CONTINUE
00161       CALL CCOPY( N, X, 1, V, 1 )
00162       ESTOLD = EST
00163       EST = SCSUM1( N, V, 1 )
00164 *
00165 *     TEST FOR CYCLING.
00166       IF( EST.LE.ESTOLD )
00167      $   GO TO 100
00168 *
00169       DO 80 I = 1, N
00170          ABSXI = ABS( X( I ) )
00171          IF( ABSXI.GT.SAFMIN ) THEN
00172             X( I ) = CMPLX( REAL( X( I ) ) / ABSXI,
00173      $               AIMAG( X( I ) ) / ABSXI )
00174          ELSE
00175             X( I ) = CONE
00176          END IF
00177    80 CONTINUE
00178       KASE = 2
00179       ISAVE( 1 ) = 4
00180       RETURN
00181 *
00182 *     ................ ENTRY   (ISAVE( 1 ) = 4)
00183 *     X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
00184 *
00185    90 CONTINUE
00186       JLAST = ISAVE( 2 )
00187       ISAVE( 2 ) = ICMAX1( N, X, 1 )
00188       IF( ( ABS( X( JLAST ) ).NE.ABS( X( ISAVE( 2 ) ) ) ) .AND.
00189      $    ( ISAVE( 3 ).LT.ITMAX ) ) THEN
00190          ISAVE( 3 ) = ISAVE( 3 ) + 1
00191          GO TO 50
00192       END IF
00193 *
00194 *     ITERATION COMPLETE.  FINAL STAGE.
00195 *
00196   100 CONTINUE
00197       ALTSGN = ONE
00198       DO 110 I = 1, N
00199          X( I ) = CMPLX( ALTSGN*( ONE + REAL( I-1 ) / REAL( N-1 ) ) )
00200          ALTSGN = -ALTSGN
00201   110 CONTINUE
00202       KASE = 1
00203       ISAVE( 1 ) = 5
00204       RETURN
00205 *
00206 *     ................ ENTRY   (ISAVE( 1 ) = 5)
00207 *     X HAS BEEN OVERWRITTEN BY A*X.
00208 *
00209   120 CONTINUE
00210       TEMP = TWO*( SCSUM1( N, X, 1 ) / REAL( 3*N ) )
00211       IF( TEMP.GT.EST ) THEN
00212          CALL CCOPY( N, X, 1, V, 1 )
00213          EST = TEMP
00214       END IF
00215 *
00216   130 CONTINUE
00217       KASE = 0
00218       RETURN
00219 *
00220 *     End of CLACN2
00221 *
00222       END