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