LAPACK 3.3.0
|
00001 SUBROUTINE ZQRT02( 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 RESULT( * ), RWORK( * ) 00013 COMPLEX*16 A( LDA, * ), AF( LDA, * ), Q( LDA, * ), 00014 $ R( LDA, * ), TAU( * ), WORK( LWORK ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * ZQRT02 tests ZUNGQR, 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, ZQRT02 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) COMPLEX*16 array, dimension (LDA,N) 00044 * The m-by-n matrix A which was factorized by ZQRT01. 00045 * 00046 * AF (input) COMPLEX*16 array, dimension (LDA,N) 00047 * Details of the QR factorization of A, as returned by ZGEQRF. 00048 * See ZGEQRF for further details. 00049 * 00050 * Q (workspace) COMPLEX*16 array, dimension (LDA,N) 00051 * 00052 * R (workspace) COMPLEX*16 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) COMPLEX*16 array, dimension (N) 00058 * The scalar factors of the elementary reflectors corresponding 00059 * to the QR factorization in AF. 00060 * 00061 * WORK (workspace) COMPLEX*16 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 COMPLEX*16 ROGUE 00079 PARAMETER ( ROGUE = ( -1.0D+10, -1.0D+10 ) ) 00080 * .. 00081 * .. Local Scalars .. 00082 INTEGER INFO 00083 DOUBLE PRECISION ANORM, EPS, RESID 00084 * .. 00085 * .. External Functions .. 00086 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY 00087 EXTERNAL DLAMCH, ZLANGE, ZLANSY 00088 * .. 00089 * .. External Subroutines .. 00090 EXTERNAL ZGEMM, ZHERK, ZLACPY, ZLASET, ZUNGQR 00091 * .. 00092 * .. Intrinsic Functions .. 00093 INTRINSIC DBLE, DCMPLX, 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 ZLASET( 'Full', M, N, ROGUE, ROGUE, Q, LDA ) 00108 CALL ZLACPY( '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 = 'ZUNGQR' 00113 CALL ZUNGQR( M, N, K, Q, LDA, TAU, WORK, LWORK, INFO ) 00114 * 00115 * Copy R(1:n,1:k) 00116 * 00117 CALL ZLASET( 'Full', N, K, DCMPLX( ZERO ), DCMPLX( ZERO ), R, 00118 $ LDA ) 00119 CALL ZLACPY( 'Upper', N, K, AF, LDA, R, LDA ) 00120 * 00121 * Compute R(1:n,1:k) - Q(1:m,1:n)' * A(1:m,1:k) 00122 * 00123 CALL ZGEMM( 'Conjugate transpose', 'No transpose', N, K, M, 00124 $ DCMPLX( -ONE ), Q, LDA, A, LDA, DCMPLX( ONE ), R, 00125 $ LDA ) 00126 * 00127 * Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . 00128 * 00129 ANORM = ZLANGE( '1', M, K, A, LDA, RWORK ) 00130 RESID = ZLANGE( '1', N, K, R, LDA, RWORK ) 00131 IF( ANORM.GT.ZERO ) THEN 00132 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M ) ) ) / ANORM ) / EPS 00133 ELSE 00134 RESULT( 1 ) = ZERO 00135 END IF 00136 * 00137 * Compute I - Q'*Q 00138 * 00139 CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ), DCMPLX( ONE ), R, LDA ) 00140 CALL ZHERK( 'Upper', 'Conjugate transpose', N, M, -ONE, Q, LDA, 00141 $ ONE, R, LDA ) 00142 * 00143 * Compute norm( I - Q'*Q ) / ( M * EPS ) . 00144 * 00145 RESID = ZLANSY( '1', 'Upper', N, R, LDA, RWORK ) 00146 * 00147 RESULT( 2 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / EPS 00148 * 00149 RETURN 00150 * 00151 * End of ZQRT02 00152 * 00153 END