LAPACK 3.3.0

zglmts.f

Go to the documentation of this file.
00001       SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
00002      $                   WORK, LWORK, RWORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            LDA, LDB, LWORK, M, N, P
00010       DOUBLE PRECISION   RESULT
00011 *     ..
00012 *     .. Array Arguments ..
00013 *
00014 *  Purpose
00015 *  =======
00016 *
00017 *  ZGLMTS tests ZGGGLM - a subroutine for solving the generalized
00018 *  linear model problem.
00019 *
00020 *  Arguments
00021 *  =========
00022 *
00023 *  N       (input) INTEGER
00024 *          The number of rows of the matrices A and B.  N >= 0.
00025 *
00026 *  M       (input) INTEGER
00027 *          The number of columns of the matrix A.  M >= 0.
00028 *
00029 *  P       (input) INTEGER
00030 *          The number of columns of the matrix B.  P >= 0.
00031 *
00032 *  A       (input) COMPLEX*16 array, dimension (LDA,M)
00033 *          The N-by-M matrix A.
00034 *
00035 *  AF      (workspace) COMPLEX*16 array, dimension (LDA,M)
00036 *
00037 *  LDA     (input) INTEGER
00038 *          The leading dimension of the arrays A, AF. LDA >= max(M,N).
00039 *
00040 *  B       (input) COMPLEX*16 array, dimension (LDB,P)
00041 *          The N-by-P matrix A.
00042 *
00043 *  BF      (workspace) COMPLEX*16 array, dimension (LDB,P)
00044 *
00045 *  LDB     (input) INTEGER
00046 *          The leading dimension of the arrays B, BF. LDB >= max(P,N).
00047 *
00048 *  D       (input) COMPLEX*16 array, dimension( N )
00049 *          On input, the left hand side of the GLM.
00050 *
00051 *  DF      (workspace) COMPLEX*16 array, dimension( N )
00052 *
00053 *  X       (output) COMPLEX*16 array, dimension( M )
00054 *          solution vector X in the GLM problem.
00055 *
00056 *  U       (output) COMPLEX*16 array, dimension( P )
00057 *          solution vector U in the GLM problem.
00058 *
00059 *  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
00060 *
00061 *  LWORK   (input) INTEGER
00062 *          The dimension of the array WORK.
00063 *
00064 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (M)
00065 *
00066 *  RESULT   (output) DOUBLE PRECISION
00067 *          The test ratio:
00068 *                           norm( d - A*x - B*u )
00069 *            RESULT = -----------------------------------------
00070 *                     (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
00071 *
00072 *  ====================================================================
00073 *
00074       DOUBLE PRECISION   RWORK( * )
00075       COMPLEX*16         A( LDA, * ), AF( LDA, * ), B( LDB, * ),
00076      $                   BF( LDB, * ), D( * ), DF( * ), U( * ),
00077      $                   WORK( LWORK ), X( * )
00078 *     ..
00079 *     .. Parameters ..
00080       DOUBLE PRECISION   ZERO
00081       PARAMETER          ( ZERO = 0.0D+0 )
00082       COMPLEX*16         CONE
00083       PARAMETER          ( CONE = 1.0D+0 )
00084 *     ..
00085 *     .. Local Scalars ..
00086       INTEGER            INFO
00087       DOUBLE PRECISION   ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
00088 *     ..
00089 *     .. External Functions ..
00090       DOUBLE PRECISION   DLAMCH, DZASUM, ZLANGE
00091       EXTERNAL           DLAMCH, DZASUM, ZLANGE
00092 *     ..
00093 *     .. External Subroutines ..
00094 *
00095       EXTERNAL           ZCOPY, ZGEMV, ZGGGLM, ZLACPY
00096 *     ..
00097 *     .. Intrinsic Functions ..
00098       INTRINSIC          MAX
00099 *     ..
00100 *     .. Executable Statements ..
00101 *
00102       EPS = DLAMCH( 'Epsilon' )
00103       UNFL = DLAMCH( 'Safe minimum' )
00104       ANORM = MAX( ZLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
00105       BNORM = MAX( ZLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
00106 *
00107 *     Copy the matrices A and B to the arrays AF and BF,
00108 *     and the vector D the array DF.
00109 *
00110       CALL ZLACPY( 'Full', N, M, A, LDA, AF, LDA )
00111       CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB )
00112       CALL ZCOPY( N, D, 1, DF, 1 )
00113 *
00114 *     Solve GLM problem
00115 *
00116       CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
00117      $             INFO )
00118 *
00119 *     Test the residual for the solution of LSE
00120 *
00121 *                       norm( d - A*x - B*u )
00122 *       RESULT = -----------------------------------------
00123 *                (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
00124 *
00125       CALL ZCOPY( N, D, 1, DF, 1 )
00126       CALL ZGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE, DF,
00127      $            1 )
00128 *
00129       CALL ZGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE, DF,
00130      $            1 )
00131 *
00132       DNORM = DZASUM( N, DF, 1 )
00133       XNORM = DZASUM( M, X, 1 ) + DZASUM( P, U, 1 )
00134       YNORM = ANORM + BNORM
00135 *
00136       IF( XNORM.LE.ZERO ) THEN
00137          RESULT = ZERO
00138       ELSE
00139          RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
00140       END IF
00141 *
00142       RETURN
00143 *
00144 *     End of ZGLMTS
00145 *
00146       END
 All Files Functions