Move LAPACK trunk into position.

1 files changed, 142 insertions, 0 deletions

diff --git a/SRC/zlaqsy.f b/SRC/zlaqsy.f
new file mode 100644
index 00000000..2a08bed3
--- /dev/null
+++ b/SRC/zlaqsy.f
@@ -0,0 +1,142 @@
+      SUBROUTINE ZLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
+*
+*  -- LAPACK auxiliary routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          EQUED, UPLO
+      INTEGER            LDA, N
+      DOUBLE PRECISION   AMAX, SCOND
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   S( * )
+      COMPLEX*16         A( LDA, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZLAQSY equilibrates a symmetric matrix A using the scaling factors
+*  in the vector S.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the upper or lower triangular part of the
+*          symmetric matrix A is stored.
+*          = 'U':  Upper triangular
+*          = 'L':  Lower triangular
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
+*          n by n upper triangular part of A contains the upper
+*          triangular part of the matrix A, and the strictly lower
+*          triangular part of A is not referenced.  If UPLO = 'L', the
+*          leading n by n lower triangular part of A contains the lower
+*          triangular part of the matrix A, and the strictly upper
+*          triangular part of A is not referenced.
+*
+*          On exit, if EQUED = 'Y', the equilibrated matrix:
+*          diag(S) * A * diag(S).
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(N,1).
+*
+*  S       (input) DOUBLE PRECISION array, dimension (N)
+*          The scale factors for A.
+*
+*  SCOND   (input) DOUBLE PRECISION
+*          Ratio of the smallest S(i) to the largest S(i).
+*
+*  AMAX    (input) DOUBLE PRECISION
+*          Absolute value of largest matrix entry.
+*
+*  EQUED   (output) CHARACTER*1
+*          Specifies whether or not equilibration was done.
+*          = 'N':  No equilibration.
+*          = 'Y':  Equilibration was done, i.e., A has been replaced by
+*                  diag(S) * A * diag(S).
+*
+*  Internal Parameters
+*  ===================
+*
+*  THRESH is a threshold value used to decide if scaling should be done
+*  based on the ratio of the scaling factors.  If SCOND < THRESH,
+*  scaling is done.
+*
+*  LARGE and SMALL are threshold values used to decide if scaling should
+*  be done based on the absolute size of the largest matrix element.
+*  If AMAX > LARGE or AMAX < SMALL, scaling is done.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, THRESH
+      PARAMETER          ( ONE = 1.0D+0, THRESH = 0.1D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J
+      DOUBLE PRECISION   CJ, LARGE, SMALL
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      DOUBLE PRECISION   DLAMCH
+      EXTERNAL           LSAME, DLAMCH
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( N.LE.0 ) THEN
+         EQUED = 'N'
+         RETURN
+      END IF
+*
+*     Initialize LARGE and SMALL.
+*
+      SMALL = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
+      LARGE = ONE / SMALL
+*
+      IF( SCOND.GE.THRESH .AND. AMAX.GE.SMALL .AND. AMAX.LE.LARGE ) THEN
+*
+*        No equilibration
+*
+         EQUED = 'N'
+      ELSE
+*
+*        Replace A by diag(S) * A * diag(S).
+*
+         IF( LSAME( UPLO, 'U' ) ) THEN
+*
+*           Upper triangle of A is stored.
+*
+            DO 20 J = 1, N
+               CJ = S( J )
+               DO 10 I = 1, J
+                  A( I, J ) = CJ*S( I )*A( I, J )
+   10          CONTINUE
+   20       CONTINUE
+         ELSE
+*
+*           Lower triangle of A is stored.
+*
+            DO 40 J = 1, N
+               CJ = S( J )
+               DO 30 I = J, N
+                  A( I, J ) = CJ*S( I )*A( I, J )
+   30          CONTINUE
+   40       CONTINUE
+         END IF
+         EQUED = 'Y'
+      END IF
+*
+      RETURN
+*
+*     End of ZLAQSY
+*
+      END