1       SUBROUTINE DQRT01P( M, N, A, AF, Q, R, LDA, TAU, WORK, LWORK,
2      $RWORK, RESULT ) 3 * 4 * -- LAPACK test routine (version 3.1) -- 5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 6 * June 2010 7 * 8 * .. Scalar Arguments .. 9 INTEGER LDA, LWORK, M, N 10 * .. 11 * .. Array Arguments .. 12 DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), Q( LDA, * ), 13$                   R( LDA, * ), RESULT* ), RWORK( * ), TAU( * ),
14      $WORK( LWORK ) 15 * .. 16 * 17 * Purpose 18 * ======= 19 * 20 * DQRT01P tests DGEQRFP, which computes the QR factorization of an m-by-n 21 * matrix A, and partially tests DORGQR which forms the m-by-m 22 * orthogonal matrix Q. 23 * 24 * DQRT01P compares R with Q'*A, and checks that Q is orthogonal. 25 * 26 * Arguments 27 * ========= 28 * 29 * M (input) INTEGER 30 * The number of rows of the matrix A. M >= 0. 31 * 32 * N (input) INTEGER 33 * The number of columns of the matrix A. N >= 0. 34 * 35 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 36 * The m-by-n matrix A. 37 * 38 * AF (output) DOUBLE PRECISION array, dimension (LDA,N) 39 * Details of the QR factorization of A, as returned by DGEQRFP. 40 * See DGEQRFP for further details. 41 * 42 * Q (output) DOUBLE PRECISION array, dimension (LDA,M) 43 * The m-by-m orthogonal matrix Q. 44 * 45 * R (workspace) DOUBLE PRECISION array, dimension (LDA,max(M,N)) 46 * 47 * LDA (input) INTEGER 48 * The leading dimension of the arrays A, AF, Q and R. 49 * LDA >= max(M,N). 50 * 51 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) 52 * The scalar factors of the elementary reflectors, as returned 53 * by DGEQRFP. 54 * 55 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 56 * 57 * LWORK (input) INTEGER 58 * The dimension of the array WORK. 59 * 60 * RWORK (workspace) DOUBLE PRECISION array, dimension (M) 61 * 62 * RESULT (output) DOUBLE PRECISION array, dimension (2) 63 * The test ratios: 64 * RESULT(1) = norm( R - Q'*A ) / ( M * norm(A) * EPS ) 65 * RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) 66 * 67 * ===================================================================== 68 * 69 * .. Parameters .. 70 DOUBLE PRECISION ZERO, ONE 71 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 72 DOUBLE PRECISION ROGUE 73 PARAMETER ( ROGUE = -1.0D+10 ) 74 * .. 75 * .. Local Scalars .. 76 INTEGER INFO, MINMN 77 DOUBLE PRECISION ANORM, EPS, RESID 78 * .. 79 * .. External Functions .. 80 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 81 EXTERNAL DLAMCH, DLANGE, DLANSY 82 * .. 83 * .. External Subroutines .. 84 EXTERNAL DGEMM, DGEQRFP, DLACPY, DLASET, DORGQR, DSYRK 85 * .. 86 * .. Intrinsic Functions .. 87 INTRINSIC DBLEMAXMIN 88 * .. 89 * .. Scalars in Common .. 90 CHARACTER*32 SRNAMT 91 * .. 92 * .. Common blocks .. 93 COMMON / SRNAMC / SRNAMT 94 * .. 95 * .. Executable Statements .. 96 * 97 MINMN = MIN( M, N ) 98 EPS = DLAMCH( 'Epsilon' ) 99 * 100 * Copy the matrix A to the array AF. 101 * 102 CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA ) 103 * 104 * Factorize the matrix A in the array AF. 105 * 106 SRNAMT = 'DGEQRFP' 107 CALL DGEQRFP( M, N, AF, LDA, TAU, WORK, LWORK, INFO ) 108 * 109 * Copy details of Q 110 * 111 CALL DLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) 112 CALL DLACPY( 'Lower', M-1, N, AF( 21 ), LDA, Q( 21 ), LDA ) 113 * 114 * Generate the m-by-m matrix Q 115 * 116 SRNAMT = 'DORGQR' 117 CALL DORGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO ) 118 * 119 * Copy R 120 * 121 CALL DLASET( 'Full', M, N, ZERO, ZERO, R, LDA ) 122 CALL DLACPY( 'Upper', M, N, AF, LDA, R, LDA ) 123 * 124 * Compute R - Q'*A 125 * 126 CALL DGEMM( 'Transpose''No transpose', M, N, M, -ONE, Q, LDA, A, 127$            LDA, ONE, R, LDA )
128 *
129 *     Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) .
130 *
131       ANORM = DLANGE( '1', M, N, A, LDA, RWORK )
132       RESID = DLANGE( '1', M, N, R, LDA, RWORK )
133       IF( ANORM.GT.ZERO ) THEN
134          RESULT1 ) = ( ( RESID / DBLEMAX1, M ) ) ) / ANORM ) / EPS
135       ELSE
136          RESULT1 ) = ZERO
137       END IF
138 *
139 *     Compute I - Q'*Q
140 *
141       CALL DLASET( 'Full', M, M, ZERO, ONE, R, LDA )
142       CALL DSYRK( 'Upper''Transpose', M, M, -ONE, Q, LDA, ONE, R,
143      \$            LDA )
144 *
145 *     Compute norm( I - Q'*Q ) / ( M * EPS ) .
146 *
147       RESID = DLANSY( '1''Upper', M, R, LDA, RWORK )
148 *
149       RESULT2 ) = ( RESID / DBLEMAX1, M ) ) ) / EPS
150 *
151       RETURN
152 *
153 *     End of DQRT01P
154 *
155       END