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