1       SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
  2      $                   X, U, WORK, LWORK, RWORK, RESULT )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            LDA, LDB, LWORK, M, P, N
 10       REAL               RESULT
 11 *     ..
 12 *     .. Array Arguments ..
 13       REAL               RWORK( * )
 14       COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
 15      $                   BF( LDB, * ), D( * ), DF( * ), U( * ),
 16      $                   WORK( LWORK ), X( * )
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  CGLMTS tests CGGGLM - a subroutine for solving the generalized
 22 *  linear model problem.
 23 *
 24 *  Arguments
 25 *  =========
 26 *
 27 *  N       (input) INTEGER
 28 *          The number of rows of the matrices A and B.  N >= 0.
 29 *
 30 *  M       (input) INTEGER
 31 *          The number of columns of the matrix A.  M >= 0.
 32 *
 33 *  P       (input) INTEGER
 34 *          The number of columns of the matrix B.  P >= 0.
 35 *
 36 *  A       (input) COMPLEX array, dimension (LDA,M)
 37 *          The N-by-M matrix A.
 38 *
 39 *  AF      (workspace) COMPLEX array, dimension (LDA,M)
 40 *
 41 *  LDA     (input) INTEGER
 42 *          The leading dimension of the arrays A, AF. LDA >= max(M,N).
 43 *
 44 *  B       (input) COMPLEX array, dimension (LDB,P)
 45 *          The N-by-P matrix A.
 46 *
 47 *  BF      (workspace) COMPLEX array, dimension (LDB,P)
 48 *
 49 *  LDB     (input) INTEGER
 50 *          The leading dimension of the arrays B, BF. LDB >= max(P,N).
 51 *
 52 *  D       (input) COMPLEX array, dimension( N )
 53 *          On input, the left hand side of the GLM.
 54 *
 55 *  DF      (workspace) COMPLEX array, dimension( N )
 56 *
 57 *  X       (output) COMPLEX array, dimension( M )
 58 *          solution vector X in the GLM problem.
 59 *
 60 *  U       (output) COMPLEX array, dimension( P )
 61 *          solution vector U in the GLM problem.
 62 *
 63 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
 64 *
 65 *  LWORK   (input) INTEGER
 66 *          The dimension of the array WORK.
 67 *
 68 *  RWORK   (workspace) REAL array, dimension (M)
 69 *
 70 *  RESULT   (output) REAL
 71 *          The test ratio:
 72 *                           norm( d - A*x - B*u )
 73 *            RESULT = -----------------------------------------
 74 *                     (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
 75 *
 76 *  ====================================================================
 77 *
 78 *     .. Parameters ..
 79       REAL               ZERO
 80       PARAMETER          ( ZERO = 0.0E+0 )
 81       COMPLEX            CONE
 82       PARAMETER          ( CONE = 1.0E+0 )
 83 *     ..
 84 *     .. Local Scalars ..
 85       INTEGER            INFO
 86       REAL               ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL
 87 *     ..
 88 *     .. External Functions ..
 89       REAL               SCASUM, SLAMCH, CLANGE
 90       EXTERNAL           SCASUM, SLAMCH, CLANGE
 91 *     ..
 92 *     .. External Subroutines ..
 93       EXTERNAL           CLACPY
 94 *
 95 *     .. Intrinsic Functions ..
 96       INTRINSIC          MAX
 97 *     ..
 98 *     .. Executable Statements ..
 99 *
100       EPS = SLAMCH( 'Epsilon' )
101       UNFL = SLAMCH( 'Safe minimum' )
102       ANORM = MAX( CLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
103       BNORM = MAX( CLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
104 *
105 *     Copy the matrices A and B to the arrays AF and BF,
106 *     and the vector D the array DF.
107 *
108       CALL CLACPY( 'Full', N, M, A, LDA, AF, LDA )
109       CALL CLACPY( 'Full', N, P, B, LDB, BF, LDB )
110       CALL CCOPY( N, D, 1, DF, 1 )
111 *
112 *     Solve GLM problem
113 *
114       CALL CGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
115      $             INFO )
116 *
117 *     Test the residual for the solution of LSE
118 *
119 *                       norm( d - A*x - B*u )
120 *       RESULT = -----------------------------------------
121 *                (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
122 *
123       CALL CCOPY( N, D, 1, DF, 1 )
124       CALL CGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE,
125      $            DF, 1 )
126 *
127       CALL CGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE,
128      $            DF, 1 )
129 *
130       DNORM = SCASUM( N, DF, 1 )
131       XNORM = SCASUM( M, X, 1 ) + SCASUM( P, U, 1 )
132       YNORM = ANORM + BNORM
133 *
134       IF( XNORM.LE.ZERO ) THEN
135          RESULT = ZERO
136       ELSE
137          RESULT =  ( ( DNORM / YNORM ) / XNORM ) /EPS
138       END IF
139 *
140       RETURN
141 *
142 *     End of CGLMTS
143 *
144       END