KTH framework for Nek5000 toolboxes; testing version  0.0.1
dorgql.f
Go to the documentation of this file.
1  SUBROUTINE dorgql( M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
2 *
3 * -- LAPACK routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * June 30, 1999
7 *
8 * .. Scalar Arguments ..
9  INTEGER INFO, K, LDA, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12  DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DORGQL generates an M-by-N real matrix Q with orthonormal columns,
19 * which is defined as the last N columns of a product of K elementary
20 * reflectors of order M
21 *
22 * Q = H(k) . . . H(2) H(1)
23 *
24 * as returned by DGEQLF.
25 *
26 * Arguments
27 * =========
28 *
29 * M (input) INTEGER
30 * The number of rows of the matrix Q. M >= 0.
31 *
32 * N (input) INTEGER
33 * The number of columns of the matrix Q. M >= N >= 0.
34 *
35 * K (input) INTEGER
36 * The number of elementary reflectors whose product defines the
37 * matrix Q. N >= K >= 0.
38 *
39 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
40 * On entry, the (n-k+i)-th column must contain the vector which
41 * defines the elementary reflector H(i), for i = 1,2,...,k, as
42 * returned by DGEQLF in the last k columns of its array
43 * argument A.
44 * On exit, the M-by-N matrix Q.
45 *
46 * LDA (input) INTEGER
47 * The first dimension of the array A. LDA >= max(1,M).
48 *
49 * TAU (input) DOUBLE PRECISION array, dimension (K)
50 * TAU(i) must contain the scalar factor of the elementary
51 * reflector H(i), as returned by DGEQLF.
52 *
53 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
54 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
55 *
56 * LWORK (input) INTEGER
57 * The dimension of the array WORK. LWORK >= max(1,N).
58 * For optimum performance LWORK >= N*NB, where NB is the
59 * optimal blocksize.
60 *
61 * If LWORK = -1, then a workspace query is assumed; the routine
62 * only calculates the optimal size of the WORK array, returns
63 * this value as the first entry of the WORK array, and no error
64 * message related to LWORK is issued by XERBLA.
65 *
66 * INFO (output) INTEGER
67 * = 0: successful exit
68 * < 0: if INFO = -i, the i-th argument has an illegal value
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73  DOUBLE PRECISION ZERO
74  parameter( zero = 0.0d+0 )
75 * ..
76 * .. Local Scalars ..
77  LOGICAL LQUERY
78  INTEGER I, IB, IINFO, IWS, J, KK, L, LDWORK, LWKOPT,
79  $ NB, NBMIN, NX
80 * ..
81 * .. External Subroutines ..
82  EXTERNAL dlarfb, dlarft, dorg2l, xerbla
83 * ..
84 * .. Intrinsic Functions ..
85  INTRINSIC max, min
86 * ..
87 * .. External Functions ..
88  INTEGER ILAENV
89  EXTERNAL ilaenv
90 * ..
91 * .. Executable Statements ..
92 *
93 * Test the input arguments
94 *
95  info = 0
96  nb = ilaenv( 1, 'DORGQL', ' ', m, n, k, -1 )
97  lwkopt = max( 1, n )*nb
98  work( 1 ) = lwkopt
99  lquery = ( lwork.EQ.-1 )
100  IF( m.LT.0 ) THEN
101  info = -1
102  ELSE IF( n.LT.0 .OR. n.GT.m ) THEN
103  info = -2
104  ELSE IF( k.LT.0 .OR. k.GT.n ) THEN
105  info = -3
106  ELSE IF( lda.LT.max( 1, m ) ) THEN
107  info = -5
108  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
109  info = -8
110  END IF
111  IF( info.NE.0 ) THEN
112  CALL xerbla( 'DORGQL', -info )
113  RETURN
114  ELSE IF( lquery ) THEN
115  RETURN
116  END IF
117 *
118 * Quick return if possible
119 *
120  IF( n.LE.0 ) THEN
121  work( 1 ) = 1
122  RETURN
123  END IF
124 *
125  nbmin = 2
126  nx = 0
127  iws = n
128  IF( nb.GT.1 .AND. nb.LT.k ) THEN
129 *
130 * Determine when to cross over from blocked to unblocked code.
131 *
132  nx = max( 0, ilaenv( 3, 'DORGQL', ' ', m, n, k, -1 ) )
133  IF( nx.LT.k ) THEN
134 *
135 * Determine if workspace is large enough for blocked code.
136 *
137  ldwork = n
138  iws = ldwork*nb
139  IF( lwork.LT.iws ) THEN
140 *
141 * Not enough workspace to use optimal NB: reduce NB and
142 * determine the minimum value of NB.
143 *
144  nb = lwork / ldwork
145  nbmin = max( 2, ilaenv( 2, 'DORGQL', ' ', m, n, k, -1 ) )
146  END IF
147  END IF
148  END IF
149 *
150  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
151 *
152 * Use blocked code after the first block.
153 * The last kk columns are handled by the block method.
154 *
155  kk = min( k, ( ( k-nx+nb-1 ) / nb )*nb )
156 *
157 * Set A(m-kk+1:m,1:n-kk) to zero.
158 *
159  DO 20 j = 1, n - kk
160  DO 10 i = m - kk + 1, m
161  a( i, j ) = zero
162  10 CONTINUE
163  20 CONTINUE
164  ELSE
165  kk = 0
166  END IF
167 *
168 * Use unblocked code for the first or only block.
169 *
170  CALL dorg2l( m-kk, n-kk, k-kk, a, lda, tau, work, iinfo )
171 *
172  IF( kk.GT.0 ) THEN
173 *
174 * Use blocked code
175 *
176  DO 50 i = k - kk + 1, k, nb
177  ib = min( nb, k-i+1 )
178  IF( n-k+i.GT.1 ) THEN
179 *
180 * Form the triangular factor of the block reflector
181 * H = H(i+ib-1) . . . H(i+1) H(i)
182 *
183  CALL dlarft( 'Backward', 'Columnwise', m-k+i+ib-1, ib,
184  $ a( 1, n-k+i ), lda, tau( i ), work, ldwork )
185 *
186 * Apply H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
187 *
188  CALL dlarfb( 'Left', 'No transpose', 'Backward',
189  $ 'Columnwise', m-k+i+ib-1, n-k+i-1, ib,
190  $ a( 1, n-k+i ), lda, work, ldwork, a, lda,
191  $ work( ib+1 ), ldwork )
192  END IF
193 *
194 * Apply H to rows 1:m-k+i+ib-1 of current block
195 *
196  CALL dorg2l( m-k+i+ib-1, ib, ib, a( 1, n-k+i ), lda,
197  $ tau( i ), work, iinfo )
198 *
199 * Set rows m-k+i+ib:m of current block to zero
200 *
201  DO 40 j = n - k + i, n - k + i + ib - 1
202  DO 30 l = m - k + i + ib, m
203  a( l, j ) = zero
204  30 CONTINUE
205  40 CONTINUE
206  50 CONTINUE
207  END IF
208 *
209  work( 1 ) = iws
210  RETURN
211 *
212 * End of DORGQL
213 *
214  END
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
Definition: dlarfb.f:3
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
Definition: dlarft.f:2
subroutine dorg2l(M, N, K, A, LDA, TAU, WORK, INFO)
Definition: dorg2l.f:2
subroutine dorgql(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dorgql.f:2
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2