1 SUBROUTINE DGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGEQRFP computes a QR factorization of a real M-by-N matrix A:
19 * A = Q * R.
20 *
21 * Arguments
22 * =========
23 *
24 * M (input) INTEGER
25 * The number of rows of the matrix A. M >= 0.
26 *
27 * N (input) INTEGER
28 * The number of columns of the matrix A. N >= 0.
29 *
30 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
31 * On entry, the M-by-N matrix A.
32 * On exit, the elements on and above the diagonal of the array
33 * contain the min(M,N)-by-N upper trapezoidal matrix R (R is
34 * upper triangular if m >= n); the elements below the diagonal,
35 * with the array TAU, represent the orthogonal matrix Q as a
36 * product of min(m,n) elementary reflectors (see Further
37 * Details).
38 *
39 * LDA (input) INTEGER
40 * The leading dimension of the array A. LDA >= max(1,M).
41 *
42 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
43 * The scalar factors of the elementary reflectors (see Further
44 * Details).
45 *
46 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
47 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
48 *
49 * LWORK (input) INTEGER
50 * The dimension of the array WORK. LWORK >= max(1,N).
51 * For optimum performance LWORK >= N*NB, where NB is
52 * the optimal blocksize.
53 *
54 * If LWORK = -1, then a workspace query is assumed; the routine
55 * only calculates the optimal size of the WORK array, returns
56 * this value as the first entry of the WORK array, and no error
57 * message related to LWORK is issued by XERBLA.
58 *
59 * INFO (output) INTEGER
60 * = 0: successful exit
61 * < 0: if INFO = -i, the i-th argument had an illegal value
62 *
63 * Further Details
64 * ===============
65 *
66 * The matrix Q is represented as a product of elementary reflectors
67 *
68 * Q = H(1) H(2) . . . H(k), where k = min(m,n).
69 *
70 * Each H(i) has the form
71 *
72 * H(i) = I - tau * v * v**T
73 *
74 * where tau is a real scalar, and v is a real vector with
75 * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
76 * and tau in TAU(i).
77 *
78 * =====================================================================
79 *
80 * .. Local Scalars ..
81 LOGICAL LQUERY
82 INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB,
83 $ NBMIN, NX
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL DGEQR2P, DLARFB, DLARFT, XERBLA
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC MAX, MIN
90 * ..
91 * .. External Functions ..
92 INTEGER ILAENV
93 EXTERNAL ILAENV
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input arguments
98 *
99 INFO = 0
100 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
101 LWKOPT = N*NB
102 WORK( 1 ) = LWKOPT
103 LQUERY = ( LWORK.EQ.-1 )
104 IF( M.LT.0 ) THEN
105 INFO = -1
106 ELSE IF( N.LT.0 ) THEN
107 INFO = -2
108 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
109 INFO = -4
110 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
111 INFO = -7
112 END IF
113 IF( INFO.NE.0 ) THEN
114 CALL XERBLA( 'DGEQRFP', -INFO )
115 RETURN
116 ELSE IF( LQUERY ) THEN
117 RETURN
118 END IF
119 *
120 * Quick return if possible
121 *
122 K = MIN( M, N )
123 IF( K.EQ.0 ) THEN
124 WORK( 1 ) = 1
125 RETURN
126 END IF
127 *
128 NBMIN = 2
129 NX = 0
130 IWS = N
131 IF( NB.GT.1 .AND. NB.LT.K ) THEN
132 *
133 * Determine when to cross over from blocked to unblocked code.
134 *
135 NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )
136 IF( NX.LT.K ) THEN
137 *
138 * Determine if workspace is large enough for blocked code.
139 *
140 LDWORK = N
141 IWS = LDWORK*NB
142 IF( LWORK.LT.IWS ) THEN
143 *
144 * Not enough workspace to use optimal NB: reduce NB and
145 * determine the minimum value of NB.
146 *
147 NB = LWORK / LDWORK
148 NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1,
149 $ -1 ) )
150 END IF
151 END IF
152 END IF
153 *
154 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
155 *
156 * Use blocked code initially
157 *
158 DO 10 I = 1, K - NX, NB
159 IB = MIN( K-I+1, NB )
160 *
161 * Compute the QR factorization of the current block
162 * A(i:m,i:i+ib-1)
163 *
164 CALL DGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
165 $ IINFO )
166 IF( I+IB.LE.N ) THEN
167 *
168 * Form the triangular factor of the block reflector
169 * H = H(i) H(i+1) . . . H(i+ib-1)
170 *
171 CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB,
172 $ A( I, I ), LDA, TAU( I ), WORK, LDWORK )
173 *
174 * Apply H**T to A(i:m,i+ib:n) from the left
175 *
176 CALL DLARFB( 'Left', 'Transpose', 'Forward',
177 $ 'Columnwise', M-I+1, N-I-IB+1, IB,
178 $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ),
179 $ LDA, WORK( IB+1 ), LDWORK )
180 END IF
181 10 CONTINUE
182 ELSE
183 I = 1
184 END IF
185 *
186 * Use unblocked code to factor the last or only block.
187 *
188 IF( I.LE.K )
189 $ CALL DGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK,
190 $ IINFO )
191 *
192 WORK( 1 ) = IWS
193 RETURN
194 *
195 * End of DGEQRFP
196 *
197 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGEQRFP computes a QR factorization of a real M-by-N matrix A:
19 * A = Q * R.
20 *
21 * Arguments
22 * =========
23 *
24 * M (input) INTEGER
25 * The number of rows of the matrix A. M >= 0.
26 *
27 * N (input) INTEGER
28 * The number of columns of the matrix A. N >= 0.
29 *
30 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
31 * On entry, the M-by-N matrix A.
32 * On exit, the elements on and above the diagonal of the array
33 * contain the min(M,N)-by-N upper trapezoidal matrix R (R is
34 * upper triangular if m >= n); the elements below the diagonal,
35 * with the array TAU, represent the orthogonal matrix Q as a
36 * product of min(m,n) elementary reflectors (see Further
37 * Details).
38 *
39 * LDA (input) INTEGER
40 * The leading dimension of the array A. LDA >= max(1,M).
41 *
42 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
43 * The scalar factors of the elementary reflectors (see Further
44 * Details).
45 *
46 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
47 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
48 *
49 * LWORK (input) INTEGER
50 * The dimension of the array WORK. LWORK >= max(1,N).
51 * For optimum performance LWORK >= N*NB, where NB is
52 * the optimal blocksize.
53 *
54 * If LWORK = -1, then a workspace query is assumed; the routine
55 * only calculates the optimal size of the WORK array, returns
56 * this value as the first entry of the WORK array, and no error
57 * message related to LWORK is issued by XERBLA.
58 *
59 * INFO (output) INTEGER
60 * = 0: successful exit
61 * < 0: if INFO = -i, the i-th argument had an illegal value
62 *
63 * Further Details
64 * ===============
65 *
66 * The matrix Q is represented as a product of elementary reflectors
67 *
68 * Q = H(1) H(2) . . . H(k), where k = min(m,n).
69 *
70 * Each H(i) has the form
71 *
72 * H(i) = I - tau * v * v**T
73 *
74 * where tau is a real scalar, and v is a real vector with
75 * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
76 * and tau in TAU(i).
77 *
78 * =====================================================================
79 *
80 * .. Local Scalars ..
81 LOGICAL LQUERY
82 INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB,
83 $ NBMIN, NX
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL DGEQR2P, DLARFB, DLARFT, XERBLA
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC MAX, MIN
90 * ..
91 * .. External Functions ..
92 INTEGER ILAENV
93 EXTERNAL ILAENV
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input arguments
98 *
99 INFO = 0
100 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
101 LWKOPT = N*NB
102 WORK( 1 ) = LWKOPT
103 LQUERY = ( LWORK.EQ.-1 )
104 IF( M.LT.0 ) THEN
105 INFO = -1
106 ELSE IF( N.LT.0 ) THEN
107 INFO = -2
108 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
109 INFO = -4
110 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
111 INFO = -7
112 END IF
113 IF( INFO.NE.0 ) THEN
114 CALL XERBLA( 'DGEQRFP', -INFO )
115 RETURN
116 ELSE IF( LQUERY ) THEN
117 RETURN
118 END IF
119 *
120 * Quick return if possible
121 *
122 K = MIN( M, N )
123 IF( K.EQ.0 ) THEN
124 WORK( 1 ) = 1
125 RETURN
126 END IF
127 *
128 NBMIN = 2
129 NX = 0
130 IWS = N
131 IF( NB.GT.1 .AND. NB.LT.K ) THEN
132 *
133 * Determine when to cross over from blocked to unblocked code.
134 *
135 NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )
136 IF( NX.LT.K ) THEN
137 *
138 * Determine if workspace is large enough for blocked code.
139 *
140 LDWORK = N
141 IWS = LDWORK*NB
142 IF( LWORK.LT.IWS ) THEN
143 *
144 * Not enough workspace to use optimal NB: reduce NB and
145 * determine the minimum value of NB.
146 *
147 NB = LWORK / LDWORK
148 NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1,
149 $ -1 ) )
150 END IF
151 END IF
152 END IF
153 *
154 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
155 *
156 * Use blocked code initially
157 *
158 DO 10 I = 1, K - NX, NB
159 IB = MIN( K-I+1, NB )
160 *
161 * Compute the QR factorization of the current block
162 * A(i:m,i:i+ib-1)
163 *
164 CALL DGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
165 $ IINFO )
166 IF( I+IB.LE.N ) THEN
167 *
168 * Form the triangular factor of the block reflector
169 * H = H(i) H(i+1) . . . H(i+ib-1)
170 *
171 CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB,
172 $ A( I, I ), LDA, TAU( I ), WORK, LDWORK )
173 *
174 * Apply H**T to A(i:m,i+ib:n) from the left
175 *
176 CALL DLARFB( 'Left', 'Transpose', 'Forward',
177 $ 'Columnwise', M-I+1, N-I-IB+1, IB,
178 $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ),
179 $ LDA, WORK( IB+1 ), LDWORK )
180 END IF
181 10 CONTINUE
182 ELSE
183 I = 1
184 END IF
185 *
186 * Use unblocked code to factor the last or only block.
187 *
188 IF( I.LE.K )
189 $ CALL DGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK,
190 $ IINFO )
191 *
192 WORK( 1 ) = IWS
193 RETURN
194 *
195 * End of DGEQRFP
196 *
197 END