LAPACK 3.3.0
|
00001 SUBROUTINE DQLT02( M, N, K, A, AF, Q, L, LDA, TAU, WORK, LWORK, 00002 $ 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 K, LDA, LWORK, M, N 00010 * .. 00011 * .. Array Arguments .. 00012 DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), L( LDA, * ), 00013 $ Q( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), 00014 $ WORK( LWORK ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DQLT02 tests DORGQL, which generates an m-by-n matrix Q with 00021 * orthonornmal columns that is defined as the product of k elementary 00022 * reflectors. 00023 * 00024 * Given the QL factorization of an m-by-n matrix A, DQLT02 generates 00025 * the orthogonal matrix Q defined by the factorization of the last k 00026 * columns of A; it compares L(m-n+1:m,n-k+1:n) with 00027 * Q(1:m,m-n+1:m)'*A(1:m,n-k+1:n), and checks that the columns of Q are 00028 * orthonormal. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * M (input) INTEGER 00034 * The number of rows of the matrix Q to be generated. M >= 0. 00035 * 00036 * N (input) INTEGER 00037 * The number of columns of the matrix Q to be generated. 00038 * M >= N >= 0. 00039 * 00040 * K (input) INTEGER 00041 * The number of elementary reflectors whose product defines the 00042 * matrix Q. N >= K >= 0. 00043 * 00044 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00045 * The m-by-n matrix A which was factorized by DQLT01. 00046 * 00047 * AF (input) DOUBLE PRECISION array, dimension (LDA,N) 00048 * Details of the QL factorization of A, as returned by DGEQLF. 00049 * See DGEQLF for further details. 00050 * 00051 * Q (workspace) DOUBLE PRECISION array, dimension (LDA,N) 00052 * 00053 * L (workspace) DOUBLE PRECISION array, dimension (LDA,N) 00054 * 00055 * LDA (input) INTEGER 00056 * The leading dimension of the arrays A, AF, Q and L. LDA >= M. 00057 * 00058 * TAU (input) DOUBLE PRECISION array, dimension (N) 00059 * The scalar factors of the elementary reflectors corresponding 00060 * to the QL factorization in AF. 00061 * 00062 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00063 * 00064 * LWORK (input) INTEGER 00065 * The dimension of the array WORK. 00066 * 00067 * RWORK (workspace) DOUBLE PRECISION array, dimension (M) 00068 * 00069 * RESULT (output) DOUBLE PRECISION array, dimension (2) 00070 * The test ratios: 00071 * RESULT(1) = norm( L - Q'*A ) / ( M * norm(A) * EPS ) 00072 * RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) 00073 * 00074 * ===================================================================== 00075 * 00076 * .. Parameters .. 00077 DOUBLE PRECISION ZERO, ONE 00078 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00079 DOUBLE PRECISION ROGUE 00080 PARAMETER ( ROGUE = -1.0D+10 ) 00081 * .. 00082 * .. Local Scalars .. 00083 INTEGER INFO 00084 DOUBLE PRECISION ANORM, EPS, RESID 00085 * .. 00086 * .. External Functions .. 00087 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00088 EXTERNAL DLAMCH, DLANGE, DLANSY 00089 * .. 00090 * .. External Subroutines .. 00091 EXTERNAL DGEMM, DLACPY, DLASET, DORGQL, DSYRK 00092 * .. 00093 * .. Intrinsic Functions .. 00094 INTRINSIC DBLE, MAX 00095 * .. 00096 * .. Scalars in Common .. 00097 CHARACTER*32 SRNAMT 00098 * .. 00099 * .. Common blocks .. 00100 COMMON / SRNAMC / SRNAMT 00101 * .. 00102 * .. Executable Statements .. 00103 * 00104 * Quick return if possible 00105 * 00106 IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN 00107 RESULT( 1 ) = ZERO 00108 RESULT( 2 ) = ZERO 00109 RETURN 00110 END IF 00111 * 00112 EPS = DLAMCH( 'Epsilon' ) 00113 * 00114 * Copy the last k columns of the factorization to the array Q 00115 * 00116 CALL DLASET( 'Full', M, N, ROGUE, ROGUE, Q, LDA ) 00117 IF( K.LT.M ) 00118 $ CALL DLACPY( 'Full', M-K, K, AF( 1, N-K+1 ), LDA, 00119 $ Q( 1, N-K+1 ), LDA ) 00120 IF( K.GT.1 ) 00121 $ CALL DLACPY( 'Upper', K-1, K-1, AF( M-K+1, N-K+2 ), LDA, 00122 $ Q( M-K+1, N-K+2 ), LDA ) 00123 * 00124 * Generate the last n columns of the matrix Q 00125 * 00126 SRNAMT = 'DORGQL' 00127 CALL DORGQL( M, N, K, Q, LDA, TAU( N-K+1 ), WORK, LWORK, INFO ) 00128 * 00129 * Copy L(m-n+1:m,n-k+1:n) 00130 * 00131 CALL DLASET( 'Full', N, K, ZERO, ZERO, L( M-N+1, N-K+1 ), LDA ) 00132 CALL DLACPY( 'Lower', K, K, AF( M-K+1, N-K+1 ), LDA, 00133 $ L( M-K+1, N-K+1 ), LDA ) 00134 * 00135 * Compute L(m-n+1:m,n-k+1:n) - Q(1:m,m-n+1:m)' * A(1:m,n-k+1:n) 00136 * 00137 CALL DGEMM( 'Transpose', 'No transpose', N, K, M, -ONE, Q, LDA, 00138 $ A( 1, N-K+1 ), LDA, ONE, L( M-N+1, N-K+1 ), LDA ) 00139 * 00140 * Compute norm( L - Q'*A ) / ( M * norm(A) * EPS ) . 00141 * 00142 ANORM = DLANGE( '1', M, K, A( 1, N-K+1 ), LDA, RWORK ) 00143 RESID = DLANGE( '1', N, K, L( M-N+1, N-K+1 ), LDA, RWORK ) 00144 IF( ANORM.GT.ZERO ) THEN 00145 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M ) ) ) / ANORM ) / EPS 00146 ELSE 00147 RESULT( 1 ) = ZERO 00148 END IF 00149 * 00150 * Compute I - Q'*Q 00151 * 00152 CALL DLASET( 'Full', N, N, ZERO, ONE, L, LDA ) 00153 CALL DSYRK( 'Upper', 'Transpose', N, M, -ONE, Q, LDA, ONE, L, 00154 $ LDA ) 00155 * 00156 * Compute norm( I - Q'*Q ) / ( M * EPS ) . 00157 * 00158 RESID = DLANSY( '1', 'Upper', N, L, LDA, RWORK ) 00159 * 00160 RESULT( 2 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / EPS 00161 * 00162 RETURN 00163 * 00164 * End of DQLT02 00165 * 00166 END