 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
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
1.7.3 
