1       SUBROUTINE CQRT01P( 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 REAL RESULT* ), RWORK( * ) 13 COMPLEX A( LDA, * ), AF( LDA, * ), Q( LDA, * ), 14$                   R( LDA, * ), TAU( * ), WORK( LWORK )
15 *     ..
16 *
17 *  Purpose
18 *  =======
19 *
20 *  CQRT01P tests CGEQRFP, which computes the QR factorization of an m-by-n
21 *  matrix A, and partially tests CUNGQR which forms the m-by-m
22 *  orthogonal matrix Q.
23 *
24 *  CQRT01P 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) COMPLEX array, dimension (LDA,N)
36 *          The m-by-n matrix A.
37 *
38 *  AF      (output) COMPLEX array, dimension (LDA,N)
39 *          Details of the QR factorization of A, as returned by CGEQRFP.
40 *          See CGEQRFP for further details.
41 *
42 *  Q       (output) COMPLEX array, dimension (LDA,M)
43 *          The m-by-m orthogonal matrix Q.
44 *
45 *  R       (workspace) COMPLEX 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) COMPLEX array, dimension (min(M,N))
52 *          The scalar factors of the elementary reflectors, as returned
53 *          by CGEQRFP.
54 *
55 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
56 *
57 *  LWORK   (input) INTEGER
58 *          The dimension of the array WORK.
59 *
60 *  RWORK   (workspace) REAL array, dimension (M)
61 *
62 *  RESULT  (output) REAL 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       REAL               ZERO, ONE
71       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
72       COMPLEX            ROGUE
73       PARAMETER          ( ROGUE = ( -1.0E+10-1.0E+10 ) )
74 *     ..
75 *     .. Local Scalars ..
76       INTEGER            INFO, MINMN
77       REAL               ANORM, EPS, RESID
78 *     ..
79 *     .. External Functions ..
80       REAL               CLANGE, CLANSY, SLAMCH
81       EXTERNAL           CLANGE, CLANSY, SLAMCH
82 *     ..
83 *     .. External Subroutines ..
84       EXTERNAL           CGEMM, CGEQRFP, CHERK, CLACPY, CLASET, CUNGQR
85 *     ..
86 *     .. Intrinsic Functions ..
87       INTRINSIC          CMPLXMAXMIN, REAL
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 = SLAMCH( 'Epsilon' )
99 *
100 *     Copy the matrix A to the array AF.
101 *
102       CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
103 *
104 *     Factorize the matrix A in the array AF.
105 *
106       SRNAMT = 'CGEQRFP'
107       CALL CGEQRFP( M, N, AF, LDA, TAU, WORK, LWORK, INFO )
108 *
109 *     Copy details of Q
110 *
111       CALL CLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA )
112       CALL CLACPY( 'Lower', M-1, N, AF( 21 ), LDA, Q( 21 ), LDA )
113 *
114 *     Generate the m-by-m matrix Q
115 *
116       SRNAMT = 'CUNGQR'
117       CALL CUNGQR( M, M, MINMN, Q, LDA, TAU, WORK, LWORK, INFO )
118 *
119 *     Copy R
120 *
121       CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), R, LDA )
122       CALL CLACPY( 'Upper', M, N, AF, LDA, R, LDA )
123 *
124 *     Compute R - Q'*A
125 *
126       CALL CGEMM( 'Conjugate transpose''No transpose', M, N, M,
127      $CMPLX-ONE ), Q, LDA, A, LDA, CMPLX( ONE ), R, LDA ) 128 * 129 * Compute norm( R - Q'*A ) / ( M * norm(A) * EPS ) . 130 * 131 ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) 132 RESID = CLANGE( '1', M, N, R, LDA, RWORK ) 133 IF( ANORM.GT.ZERO ) THEN 134 RESULT1 ) = ( ( RESID / REALMAX1, M ) ) ) / ANORM ) / EPS 135 ELSE 136 RESULT1 ) = ZERO 137 END IF 138 * 139 * Compute I - Q'*Q 140 * 141 CALL CLASET( 'Full', M, M, CMPLX( ZERO ), CMPLX( ONE ), R, LDA ) 142 CALL CHERK( 'Upper''Conjugate transpose', M, M, -ONE, Q, LDA, 143$            ONE, R, LDA )
144 *
145 *     Compute norm( I - Q'*Q ) / ( M * EPS ) .
146 *
147       RESID = CLANSY( '1''Upper', M, R, LDA, RWORK )
148 *
149       RESULT2 ) = ( RESID / REALMAX1, M ) ) ) / EPS
150 *
151       RETURN
152 *
153 *     End of CQRT01P
154 *
155       END