Move LAPACK trunk into position.

1 files changed, 144 insertions, 0 deletions

diff --git a/TESTING/EIG/cglmts.f b/TESTING/EIG/cglmts.f
new file mode 100644
index 00000000..84d6b8df
--- /dev/null
+++ b/TESTING/EIG/cglmts.f
@@ -0,0 +1,144 @@
+      SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
+     $                   X, U, WORK, LWORK, RWORK, RESULT )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            LDA, LDB, LWORK, M, P, N
+      REAL               RESULT
+*     ..
+*     .. Array Arguments ..
+      REAL               RWORK( * )
+      COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
+     $                   BF( LDB, * ), D( * ), DF( * ), U( * ),
+     $                   WORK( LWORK ), X( * )
+*
+*  Purpose
+*  =======
+*
+*  CGLMTS tests CGGGLM - a subroutine for solving the generalized
+*  linear model problem.
+*
+*  Arguments
+*  =========
+*
+*  N       (input) INTEGER
+*          The number of rows of the matrices A and B.  N >= 0.
+*
+*  M       (input) INTEGER
+*          The number of columns of the matrix A.  M >= 0.
+*
+*  P       (input) INTEGER
+*          The number of columns of the matrix B.  P >= 0.
+*
+*  A       (input) COMPLEX array, dimension (LDA,M)
+*          The N-by-M matrix A.
+*
+*  AF      (workspace) COMPLEX array, dimension (LDA,M)
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the arrays A, AF. LDA >= max(M,N).
+*
+*  B       (input) COMPLEX array, dimension (LDB,P)
+*          The N-by-P matrix A.
+*
+*  BF      (workspace) COMPLEX array, dimension (LDB,P)
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the arrays B, BF. LDB >= max(P,N).
+*
+*  D       (input) COMPLEX array, dimension( N )
+*          On input, the left hand side of the GLM.
+*
+*  DF      (workspace) COMPLEX array, dimension( N )
+*
+*  X       (output) COMPLEX array, dimension( M )
+*          solution vector X in the GLM problem.
+*
+*  U       (output) COMPLEX array, dimension( P )
+*          solution vector U in the GLM problem.
+*
+*  WORK    (workspace) COMPLEX array, dimension (LWORK)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK.
+*
+*  RWORK   (workspace) REAL array, dimension (M)
+*
+*  RESULT   (output) REAL
+*          The test ratio:
+*                           norm( d - A*x - B*u )
+*            RESULT = -----------------------------------------
+*                     (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
+*
+*  ====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E+0 )
+      COMPLEX            CONE
+      PARAMETER          ( CONE = 1.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            INFO
+      REAL               ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL
+*     ..
+*     .. External Functions ..
+      REAL               SCASUM, SLAMCH, CLANGE
+      EXTERNAL           SCASUM, SLAMCH, CLANGE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLACPY
+*
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+      EPS = SLAMCH( 'Epsilon' )
+      UNFL = SLAMCH( 'Safe minimum' )
+      ANORM = MAX( CLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
+      BNORM = MAX( CLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
+*
+*     Copy the matrices A and B to the arrays AF and BF,
+*     and the vector D the array DF.
+*
+      CALL CLACPY( 'Full', N, M, A, LDA, AF, LDA )
+      CALL CLACPY( 'Full', N, P, B, LDB, BF, LDB )
+      CALL CCOPY( N, D, 1, DF, 1 )
+*
+*     Solve GLM problem
+*
+      CALL CGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
+     $             INFO )
+*
+*     Test the residual for the solution of LSE
+*
+*                       norm( d - A*x - B*u )
+*       RESULT = -----------------------------------------
+*                (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
+*
+      CALL CCOPY( N, D, 1, DF, 1 )
+      CALL CGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE,
+     $            DF, 1 )
+*
+      CALL CGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE,
+     $            DF, 1 )
+*
+      DNORM = SCASUM( N, DF, 1 )
+      XNORM = SCASUM( M, X, 1 ) + SCASUM( P, U, 1 )
+      YNORM = ANORM + BNORM
+*
+      IF( XNORM.LE.ZERO ) THEN
+         RESULT = ZERO
+      ELSE
+         RESULT =  ( ( DNORM / YNORM ) / XNORM ) /EPS
+      END IF
+*
+      RETURN
+*
+*     End of CGLMTS
+*
+      END