d5/d96/dglmts_8f_source.html

00001       SUBROUTINE DGLMTS( 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 *  DGLMTS tests DGGGLM - 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) DOUBLE PRECISION array, dimension (LDA,M)
00033 *          The N-by-M matrix A.
00034 *
00035 *  AF      (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDB,P)
00041 *          The N-by-P matrix A.
00042 *
00043 *  BF      (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension( N )
00049 *          On input, the left hand side of the GLM.
00050 *
00051 *  DF      (workspace) DOUBLE PRECISION array, dimension( N )
00052 *
00053 *  X       (output) DOUBLE PRECISION array, dimension( M )
00054 *          solution vector X in the GLM problem.
00055 *
00056 *  U       (output) DOUBLE PRECISION array, dimension( P )
00057 *          solution vector U in the GLM problem.
00058 *
00059 *  WORK    (workspace) DOUBLE PRECISION 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   A( LDA, * ), AF( LDA, * ), B( LDB, * ),
00075      $                   BF( LDB, * ), D( * ), DF( * ), RWORK( * ),
00076      $                   U( * ), WORK( LWORK ), X( * )
00077 *     ..
00078 *     .. Parameters ..
00079       DOUBLE PRECISION   ZERO, ONE
00080       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00081 *     ..
00082 *     .. Local Scalars ..
00083       INTEGER            INFO
00084       DOUBLE PRECISION   ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
00085 *     ..
00086 *     .. External Functions ..
00087       DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
00088       EXTERNAL           DASUM, DLAMCH, DLANGE
00089 *     ..
00090 *     .. External Subroutines ..
00091 *
00092       EXTERNAL           DCOPY, DGEMV, DGGGLM, DLACPY
00093 *     ..
00094 *     .. Intrinsic Functions ..
00095       INTRINSIC          MAX
00096 *     ..
00097 *     .. Executable Statements ..
00098 *
00099       EPS = DLAMCH( 'Epsilon' )
00100       UNFL = DLAMCH( 'Safe minimum' )
00101       ANORM = MAX( DLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
00102       BNORM = MAX( DLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
00103 *
00104 *     Copy the matrices A and B to the arrays AF and BF,
00105 *     and the vector D the array DF.
00106 *
00107       CALL DLACPY( 'Full', N, M, A, LDA, AF, LDA )
00108       CALL DLACPY( 'Full', N, P, B, LDB, BF, LDB )
00109       CALL DCOPY( N, D, 1, DF, 1 )
00110 *
00111 *     Solve GLM problem
00112 *
00113       CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
00114      $             INFO )
00115 *
00116 *     Test the residual for the solution of LSE
00117 *
00118 *                       norm( d - A*x - B*u )
00119 *       RESULT = -----------------------------------------
00120 *                (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
00121 *
00122       CALL DCOPY( N, D, 1, DF, 1 )
00123       CALL DGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE, DF, 1 )
00124 *
00125       CALL DGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE, DF, 1 )
00126 *
00127       DNORM = DASUM( N, DF, 1 )
00128       XNORM = DASUM( M, X, 1 ) + DASUM( P, U, 1 )
00129       YNORM = ANORM + BNORM
00130 *
00131       IF( XNORM.LE.ZERO ) THEN
00132          RESULT = ZERO
00133       ELSE
00134          RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
00135       END IF
00136 *
00137       RETURN
00138 *
00139 *     End of DGLMTS
00140 *
00141       END