1 DOUBLE PRECISION FUNCTION ZQPT01( M, N, K, A, AF, LDA, TAU, JPVT,
2 $ WORK, LWORK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER K, LDA, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER JPVT( * )
13 COMPLEX*16 A( LDA, * ), AF( LDA, * ), TAU( * ),
14 $ WORK( LWORK )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZQPT01 tests the QR-factorization with pivoting of a matrix A. The
21 * array AF contains the (possibly partial) QR-factorization of A, where
22 * the upper triangle of AF(1:k,1:k) is a partial triangular factor,
23 * the entries below the diagonal in the first k columns are the
24 * Householder vectors, and the rest of AF contains a partially updated
25 * matrix.
26 *
27 * This function returns ||A*P - Q*R||/(||norm(A)||*eps*M)
28 *
29 * Arguments
30 * =========
31 *
32 * M (input) INTEGER
33 * The number of rows of the matrices A and AF.
34 *
35 * N (input) INTEGER
36 * The number of columns of the matrices A and AF.
37 *
38 * K (input) INTEGER
39 * The number of columns of AF that have been reduced
40 * to upper triangular form.
41 *
42 * A (input) COMPLEX*16 array, dimension (LDA, N)
43 * The original matrix A.
44 *
45 * AF (input) COMPLEX*16 array, dimension (LDA,N)
46 * The (possibly partial) output of ZGEQPF. The upper triangle
47 * of AF(1:k,1:k) is a partial triangular factor, the entries
48 * below the diagonal in the first k columns are the Householder
49 * vectors, and the rest of AF contains a partially updated
50 * matrix.
51 *
52 * LDA (input) INTEGER
53 * The leading dimension of the arrays A and AF.
54 *
55 * TAU (input) COMPLEX*16 array, dimension (K)
56 * Details of the Householder transformations as returned by
57 * ZGEQPF.
58 *
59 * JPVT (input) INTEGER array, dimension (N)
60 * Pivot information as returned by ZGEQPF.
61 *
62 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
63 *
64 * LWORK (input) INTEGER
65 * The length of the array WORK. LWORK >= M*N+N.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 DOUBLE PRECISION ZERO, ONE
71 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
72 * ..
73 * .. Local Scalars ..
74 INTEGER I, INFO, J
75 DOUBLE PRECISION NORMA
76 * ..
77 * .. Local Arrays ..
78 DOUBLE PRECISION RWORK( 1 )
79 * ..
80 * .. External Functions ..
81 DOUBLE PRECISION DLAMCH, ZLANGE
82 EXTERNAL DLAMCH, ZLANGE
83 * ..
84 * .. External Subroutines ..
85 EXTERNAL XERBLA, ZAXPY, ZCOPY, ZUNMQR
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC DBLE, DCMPLX, MAX, MIN
89 * ..
90 * .. Executable Statements ..
91 *
92 ZQPT01 = ZERO
93 *
94 * Test if there is enough workspace
95 *
96 IF( LWORK.LT.M*N+N ) THEN
97 CALL XERBLA( 'ZQPT01', 10 )
98 RETURN
99 END IF
100 *
101 * Quick return if possible
102 *
103 IF( M.LE.0 .OR. N.LE.0 )
104 $ RETURN
105 *
106 NORMA = ZLANGE( 'One-norm', M, N, A, LDA, RWORK )
107 *
108 DO 30 J = 1, K
109 DO 10 I = 1, MIN( J, M )
110 WORK( ( J-1 )*M+I ) = AF( I, J )
111 10 CONTINUE
112 DO 20 I = J + 1, M
113 WORK( ( J-1 )*M+I ) = ZERO
114 20 CONTINUE
115 30 CONTINUE
116 DO 40 J = K + 1, N
117 CALL ZCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
118 40 CONTINUE
119 *
120 CALL ZUNMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
121 $ M, WORK( M*N+1 ), LWORK-M*N, INFO )
122 *
123 DO 50 J = 1, N
124 *
125 * Compare i-th column of QR and jpvt(i)-th column of A
126 *
127 CALL ZAXPY( M, DCMPLX( -ONE ), A( 1, JPVT( J ) ), 1,
128 $ WORK( ( J-1 )*M+1 ), 1 )
129 50 CONTINUE
130 *
131 ZQPT01 = ZLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
132 $ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) )
133 IF( NORMA.NE.ZERO )
134 $ ZQPT01 = ZQPT01 / NORMA
135 *
136 RETURN
137 *
138 * End of ZQPT01
139 *
140 END
2 $ WORK, LWORK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER K, LDA, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER JPVT( * )
13 COMPLEX*16 A( LDA, * ), AF( LDA, * ), TAU( * ),
14 $ WORK( LWORK )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZQPT01 tests the QR-factorization with pivoting of a matrix A. The
21 * array AF contains the (possibly partial) QR-factorization of A, where
22 * the upper triangle of AF(1:k,1:k) is a partial triangular factor,
23 * the entries below the diagonal in the first k columns are the
24 * Householder vectors, and the rest of AF contains a partially updated
25 * matrix.
26 *
27 * This function returns ||A*P - Q*R||/(||norm(A)||*eps*M)
28 *
29 * Arguments
30 * =========
31 *
32 * M (input) INTEGER
33 * The number of rows of the matrices A and AF.
34 *
35 * N (input) INTEGER
36 * The number of columns of the matrices A and AF.
37 *
38 * K (input) INTEGER
39 * The number of columns of AF that have been reduced
40 * to upper triangular form.
41 *
42 * A (input) COMPLEX*16 array, dimension (LDA, N)
43 * The original matrix A.
44 *
45 * AF (input) COMPLEX*16 array, dimension (LDA,N)
46 * The (possibly partial) output of ZGEQPF. The upper triangle
47 * of AF(1:k,1:k) is a partial triangular factor, the entries
48 * below the diagonal in the first k columns are the Householder
49 * vectors, and the rest of AF contains a partially updated
50 * matrix.
51 *
52 * LDA (input) INTEGER
53 * The leading dimension of the arrays A and AF.
54 *
55 * TAU (input) COMPLEX*16 array, dimension (K)
56 * Details of the Householder transformations as returned by
57 * ZGEQPF.
58 *
59 * JPVT (input) INTEGER array, dimension (N)
60 * Pivot information as returned by ZGEQPF.
61 *
62 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
63 *
64 * LWORK (input) INTEGER
65 * The length of the array WORK. LWORK >= M*N+N.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 DOUBLE PRECISION ZERO, ONE
71 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
72 * ..
73 * .. Local Scalars ..
74 INTEGER I, INFO, J
75 DOUBLE PRECISION NORMA
76 * ..
77 * .. Local Arrays ..
78 DOUBLE PRECISION RWORK( 1 )
79 * ..
80 * .. External Functions ..
81 DOUBLE PRECISION DLAMCH, ZLANGE
82 EXTERNAL DLAMCH, ZLANGE
83 * ..
84 * .. External Subroutines ..
85 EXTERNAL XERBLA, ZAXPY, ZCOPY, ZUNMQR
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC DBLE, DCMPLX, MAX, MIN
89 * ..
90 * .. Executable Statements ..
91 *
92 ZQPT01 = ZERO
93 *
94 * Test if there is enough workspace
95 *
96 IF( LWORK.LT.M*N+N ) THEN
97 CALL XERBLA( 'ZQPT01', 10 )
98 RETURN
99 END IF
100 *
101 * Quick return if possible
102 *
103 IF( M.LE.0 .OR. N.LE.0 )
104 $ RETURN
105 *
106 NORMA = ZLANGE( 'One-norm', M, N, A, LDA, RWORK )
107 *
108 DO 30 J = 1, K
109 DO 10 I = 1, MIN( J, M )
110 WORK( ( J-1 )*M+I ) = AF( I, J )
111 10 CONTINUE
112 DO 20 I = J + 1, M
113 WORK( ( J-1 )*M+I ) = ZERO
114 20 CONTINUE
115 30 CONTINUE
116 DO 40 J = K + 1, N
117 CALL ZCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
118 40 CONTINUE
119 *
120 CALL ZUNMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
121 $ M, WORK( M*N+1 ), LWORK-M*N, INFO )
122 *
123 DO 50 J = 1, N
124 *
125 * Compare i-th column of QR and jpvt(i)-th column of A
126 *
127 CALL ZAXPY( M, DCMPLX( -ONE ), A( 1, JPVT( J ) ), 1,
128 $ WORK( ( J-1 )*M+1 ), 1 )
129 50 CONTINUE
130 *
131 ZQPT01 = ZLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
132 $ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) )
133 IF( NORMA.NE.ZERO )
134 $ ZQPT01 = ZQPT01 / NORMA
135 *
136 RETURN
137 *
138 * End of ZQPT01
139 *
140 END