Move LAPACK trunk into position.

1 files changed, 221 insertions, 0 deletions

diff --git a/SRC/zlacn2.f b/SRC/zlacn2.f
new file mode 100644
index 00000000..99f7ae35
--- /dev/null
+++ b/SRC/zlacn2.f
@@ -0,0 +1,221 @@
+      SUBROUTINE ZLACN2( N, V, X, EST, KASE, ISAVE )
+*
+*  -- LAPACK auxiliary routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            KASE, N
+      DOUBLE PRECISION   EST
+*     ..
+*     .. Array Arguments ..
+      INTEGER            ISAVE( 3 )
+      COMPLEX*16         V( * ), X( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZLACN2 estimates the 1-norm of a square, complex matrix A.
+*  Reverse communication is used for evaluating matrix-vector products.
+*
+*  Arguments
+*  =========
+*
+*  N      (input) INTEGER
+*         The order of the matrix.  N >= 1.
+*
+*  V      (workspace) COMPLEX*16 array, dimension (N)
+*         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
+*         (W is not returned).
+*
+*  X      (input/output) COMPLEX*16 array, dimension (N)
+*         On an intermediate return, X should be overwritten by
+*               A * X,   if KASE=1,
+*               A' * X,  if KASE=2,
+*         where A' is the conjugate transpose of A, and ZLACN2 must be
+*         re-called with all the other parameters unchanged.
+*
+*  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 ZLACN2.
+*         On exit, EST is an estimate (a lower bound) for norm(A). 
+*
+*  KASE   (input/output) INTEGER
+*         On the initial call to ZLACN2, 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 ZLACN2, KASE will again be 0.
+*
+*  ISAVE  (input/output) INTEGER array, dimension (3)
+*         ISAVE is used to save variables between calls to ZLACN2
+*
+*  Further Details
+*  ======= =======
+*
+*  Contributed by Nick Higham, University of Manchester.
+*  Originally named CONEST, 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.
+*
+*  Last modified:  April, 1999
+*
+*  This is a thread safe version of ZLACON, which uses the array ISAVE
+*  in place of a SAVE statement, as follows:
+*
+*     ZLACON     ZLACN2
+*      JUMP     ISAVE(1)
+*      J        ISAVE(2)
+*      ITER     ISAVE(3)
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER              ITMAX
+      PARAMETER          ( ITMAX = 5 )
+      DOUBLE PRECISION     ONE,         TWO
+      PARAMETER          ( ONE = 1.0D0, TWO = 2.0D0 )
+      COMPLEX*16           CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
+     $                            CONE = ( 1.0D0, 0.0D0 ) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, JLAST
+      DOUBLE PRECISION   ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
+*     ..
+*     .. External Functions ..
+      INTEGER            IZMAX1
+      DOUBLE PRECISION   DLAMCH, DZSUM1
+      EXTERNAL           IZMAX1, DLAMCH, DZSUM1
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZCOPY
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DBLE, DCMPLX, DIMAG
+*     ..
+*     .. Executable Statements ..
+*
+      SAFMIN = DLAMCH( 'Safe minimum' )
+      IF( KASE.EQ.0 ) THEN
+         DO 10 I = 1, N
+            X( I ) = DCMPLX( ONE / DBLE( N ) )
+   10    CONTINUE
+         KASE = 1
+         ISAVE( 1 ) = 1
+         RETURN
+      END IF
+*
+      GO TO ( 20, 40, 70, 90, 120 )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 130
+      END IF
+      EST = DZSUM1( N, X, 1 )
+*
+      DO 30 I = 1, N
+         ABSXI = ABS( X( I ) )
+         IF( ABSXI.GT.SAFMIN ) THEN
+            X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
+     $               DIMAG( X( I ) ) / ABSXI )
+         ELSE
+            X( I ) = CONE
+         END IF
+   30 CONTINUE
+      KASE = 2
+      ISAVE( 1 ) = 2
+      RETURN
+*
+*     ................ ENTRY   (ISAVE( 1 ) = 2)
+*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
+*
+   40 CONTINUE
+      ISAVE( 2 ) = IZMAX1( N, X, 1 )
+      ISAVE( 3 ) = 2
+*
+*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
+*
+   50 CONTINUE
+      DO 60 I = 1, N
+         X( I ) = CZERO
+   60 CONTINUE
+      X( ISAVE( 2 ) ) = CONE
+      KASE = 1
+      ISAVE( 1 ) = 3
+      RETURN
+*
+*     ................ ENTRY   (ISAVE( 1 ) = 3)
+*     X HAS BEEN OVERWRITTEN BY A*X.
+*
+   70 CONTINUE
+      CALL ZCOPY( N, X, 1, V, 1 )
+      ESTOLD = EST
+      EST = DZSUM1( N, V, 1 )
+*
+*     TEST FOR CYCLING.
+      IF( EST.LE.ESTOLD )
+     $   GO TO 100
+*
+      DO 80 I = 1, N
+         ABSXI = ABS( X( I ) )
+         IF( ABSXI.GT.SAFMIN ) THEN
+            X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
+     $               DIMAG( X( I ) ) / ABSXI )
+         ELSE
+            X( I ) = CONE
+         END IF
+   80 CONTINUE
+      KASE = 2
+      ISAVE( 1 ) = 4
+      RETURN
+*
+*     ................ ENTRY   (ISAVE( 1 ) = 4)
+*     X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
+*
+   90 CONTINUE
+      JLAST = ISAVE( 2 )
+      ISAVE( 2 ) = IZMAX1( N, X, 1 )
+      IF( ( ABS( 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.
+*
+  100 CONTINUE
+      ALTSGN = ONE
+      DO 110 I = 1, N
+         X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
+         ALTSGN = -ALTSGN
+  110 CONTINUE
+      KASE = 1
+      ISAVE( 1 ) = 5
+      RETURN
+*
+*     ................ ENTRY   (ISAVE( 1 ) = 5)
+*     X HAS BEEN OVERWRITTEN BY A*X.
+*
+  120 CONTINUE
+      TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
+      IF( TEMP.GT.EST ) THEN
+         CALL ZCOPY( N, X, 1, V, 1 )
+         EST = TEMP
+      END IF
+*
+  130 CONTINUE
+      KASE = 0
+      RETURN
+*
+*     End of ZLACN2
+*
+      END