Move LAPACK trunk into position.

1 files changed, 172 insertions, 0 deletions

diff --git a/TESTING/LIN/zgbt01.f b/TESTING/LIN/zgbt01.f
new file mode 100644
index 00000000..f792850a
--- /dev/null
+++ b/TESTING/LIN/zgbt01.f
@@ -0,0 +1,172 @@
+      SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
+     $                   RESID )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            KL, KU, LDA, LDAFAC, M, N
+      DOUBLE PRECISION   RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZGBT01 reconstructs a band matrix  A  from its L*U factorization and
+*  computes the residual:
+*     norm(L*U - A) / ( N * norm(A) * EPS ),
+*  where EPS is the machine epsilon.
+*
+*  The expression L*U - A is computed one column at a time, so A and
+*  AFAC are not modified.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix A.  M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix A.  N >= 0.
+*
+*  KL      (input) INTEGER
+*          The number of subdiagonals within the band of A.  KL >= 0.
+*
+*  KU      (input) INTEGER
+*          The number of superdiagonals within the band of A.  KU >= 0.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          The original matrix A in band storage, stored in rows 1 to
+*          KL+KU+1.
+*
+*  LDA     (input) INTEGER.
+*          The leading dimension of the array A.  LDA >= max(1,KL+KU+1).
+*
+*  AFAC    (input) COMPLEX*16 array, dimension (LDAFAC,N)
+*          The factored form of the matrix A.  AFAC contains the banded
+*          factors L and U from the L*U factorization, as computed by
+*          ZGBTRF.  U is stored as an upper triangular band matrix with
+*          KL+KU superdiagonals in rows 1 to KL+KU+1, and the
+*          multipliers used during the factorization are stored in rows
+*          KL+KU+2 to 2*KL+KU+1.  See ZGBTRF for further details.
+*
+*  LDAFAC  (input) INTEGER
+*          The leading dimension of the array AFAC.
+*          LDAFAC >= max(1,2*KL*KU+1).
+*
+*  IPIV    (input) INTEGER array, dimension (min(M,N))
+*          The pivot indices from ZGBTRF.
+*
+*  WORK    (workspace) COMPLEX*16 array, dimension (2*KL+KU+1)
+*
+*  RESID   (output) DOUBLE PRECISION
+*          norm(L*U - A) / ( N * norm(A) * EPS )
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
+      DOUBLE PRECISION   ANORM, EPS
+      COMPLEX*16         T
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMCH, DZASUM
+      EXTERNAL           DLAMCH, DZASUM
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZAXPY, ZCOPY
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          DBLE, DCMPLX, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick exit if M = 0 or N = 0.
+*
+      RESID = ZERO
+      IF( M.LE.0 .OR. N.LE.0 )
+     $   RETURN
+*
+*     Determine EPS and the norm of A.
+*
+      EPS = DLAMCH( 'Epsilon' )
+      KD = KU + 1
+      ANORM = ZERO
+      DO 10 J = 1, N
+         I1 = MAX( KD+1-J, 1 )
+         I2 = MIN( KD+M-J, KL+KD )
+         IF( I2.GE.I1 )
+     $      ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
+   10 CONTINUE
+*
+*     Compute one column at a time of L*U - A.
+*
+      KD = KL + KU + 1
+      DO 40 J = 1, N
+*
+*        Copy the J-th column of U to WORK.
+*
+         JU = MIN( KL+KU, J-1 )
+         JL = MIN( KL, M-J )
+         LENJ = MIN( M, J ) - J + JU + 1
+         IF( LENJ.GT.0 ) THEN
+            CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
+            DO 20 I = LENJ + 1, JU + JL + 1
+               WORK( I ) = ZERO
+   20       CONTINUE
+*
+*           Multiply by the unit lower triangular matrix L.  Note that L
+*           is stored as a product of transformations and permutations.
+*
+            DO 30 I = MIN( M-1, J ), J - JU, -1
+               IL = MIN( KL, M-I )
+               IF( IL.GT.0 ) THEN
+                  IW = I - J + JU + 1
+                  T = WORK( IW )
+                  CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
+     $                        1 )
+                  IP = IPIV( I )
+                  IF( I.NE.IP ) THEN
+                     IP = IP - J + JU + 1
+                     WORK( IW ) = WORK( IP )
+                     WORK( IP ) = T
+                  END IF
+               END IF
+   30       CONTINUE
+*
+*           Subtract the corresponding column of A.
+*
+            JUA = MIN( JU, KU )
+            IF( JUA+JL+1.GT.0 )
+     $         CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ),
+     $                     1, WORK( JU+1-JUA ), 1 )
+*
+*           Compute the 1-norm of the column.
+*
+            RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) )
+         END IF
+   40 CONTINUE
+*
+*     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
+*
+      IF( ANORM.LE.ZERO ) THEN
+         IF( RESID.NE.ZERO )
+     $      RESID = ONE / EPS
+      ELSE
+         RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
+      END IF
+*
+      RETURN
+*
+*     End of ZGBT01
+*
+      END