Mercurial > octave

comparison libcruft/lapack/sormrz.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE SORMRZ( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
2 $ WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.1.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * January 2007
7 *
8 * .. Scalar Arguments ..
9 CHARACTER SIDE, TRANS
10 INTEGER INFO, K, L, LDA, LDC, LWORK, M, N
11 * ..
12 * .. Array Arguments ..
13 REAL A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SORMRZ overwrites the general real M-by-N matrix C with
20 *
21 * SIDE = 'L' SIDE = 'R'
22 * TRANS = 'N': Q * C C * Q
23 * TRANS = 'T': Q**T * C C * Q**T
24 *
25 * where Q is a real orthogonal matrix defined as the product of k
26 * elementary reflectors
27 *
28 * Q = H(1) H(2) . . . H(k)
29 *
30 * as returned by STZRZF. Q is of order M if SIDE = 'L' and of order N
31 * if SIDE = 'R'.
32 *
33 * Arguments
34 * =========
35 *
36 * SIDE (input) CHARACTER*1
37 * = 'L': apply Q or Q**T from the Left;
38 * = 'R': apply Q or Q**T from the Right.
39 *
40 * TRANS (input) CHARACTER*1
41 * = 'N': No transpose, apply Q;
42 * = 'T': Transpose, apply Q**T.
43 *
44 * M (input) INTEGER
45 * The number of rows of the matrix C. M >= 0.
46 *
47 * N (input) INTEGER
48 * The number of columns of the matrix C. N >= 0.
49 *
50 * K (input) INTEGER
51 * The number of elementary reflectors whose product defines
52 * the matrix Q.
53 * If SIDE = 'L', M >= K >= 0;
54 * if SIDE = 'R', N >= K >= 0.
55 *
56 * L (input) INTEGER
57 * The number of columns of the matrix A containing
58 * the meaningful part of the Householder reflectors.
59 * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
60 *
61 * A (input) REAL array, dimension
62 * (LDA,M) if SIDE = 'L',
63 * (LDA,N) if SIDE = 'R'
64 * The i-th row must contain the vector which defines the
65 * elementary reflector H(i), for i = 1,2,...,k, as returned by
66 * STZRZF in the last k rows of its array argument A.
67 * A is modified by the routine but restored on exit.
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,K).
71 *
72 * TAU (input) REAL array, dimension (K)
73 * TAU(i) must contain the scalar factor of the elementary
74 * reflector H(i), as returned by STZRZF.
75 *
76 * C (input/output) REAL array, dimension (LDC,N)
77 * On entry, the M-by-N matrix C.
78 * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
79 *
80 * LDC (input) INTEGER
81 * The leading dimension of the array C. LDC >= max(1,M).
82 *
83 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
84 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
85 *
86 * LWORK (input) INTEGER
87 * The dimension of the array WORK.
88 * If SIDE = 'L', LWORK >= max(1,N);
89 * if SIDE = 'R', LWORK >= max(1,M).
90 * For optimum performance LWORK >= N*NB if SIDE = 'L', and
91 * LWORK >= M*NB if SIDE = 'R', where NB is the optimal
92 * blocksize.
93 *
94 * If LWORK = -1, then a workspace query is assumed; the routine
95 * only calculates the optimal size of the WORK array, returns
96 * this value as the first entry of the WORK array, and no error
97 * message related to LWORK is issued by XERBLA.
98 *
99 * INFO (output) INTEGER
100 * = 0: successful exit
101 * < 0: if INFO = -i, the i-th argument had an illegal value
102 *
103 * Further Details
104 * ===============
105 *
106 * Based on contributions by
107 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
108 *
109 * =====================================================================
110 *
111 * .. Parameters ..
112 INTEGER NBMAX, LDT
113 PARAMETER ( NBMAX = 64, LDT = NBMAX+1 )
114 * ..
115 * .. Local Scalars ..
116 LOGICAL LEFT, LQUERY, NOTRAN
117 CHARACTER TRANST
118 INTEGER I, I1, I2, I3, IB, IC, IINFO, IWS, JA, JC,
119 $ LDWORK, LWKOPT, MI, NB, NBMIN, NI, NQ, NW
120 * ..
121 * .. Local Arrays ..
122 REAL T( LDT, NBMAX )
123 * ..
124 * .. External Functions ..
125 LOGICAL LSAME
126 INTEGER ILAENV
127 EXTERNAL LSAME, ILAENV
128 * ..
129 * .. External Subroutines ..
130 EXTERNAL SLARZB, SLARZT, SORMR3, XERBLA
131 * ..
132 * .. Intrinsic Functions ..
133 INTRINSIC MAX, MIN
134 * ..
135 * .. Executable Statements ..
136 *
137 * Test the input arguments
138 *
139 INFO = 0
140 LEFT = LSAME( SIDE, 'L' )
141 NOTRAN = LSAME( TRANS, 'N' )
142 LQUERY = ( LWORK.EQ.-1 )
143 *
144 * NQ is the order of Q and NW is the minimum dimension of WORK
145 *
146 IF( LEFT ) THEN
147 NQ = M
148 NW = MAX( 1, N )
149 ELSE
150 NQ = N
151 NW = MAX( 1, M )
152 END IF
153 IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
154 INFO = -1
155 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
156 INFO = -2
157 ELSE IF( M.LT.0 ) THEN
158 INFO = -3
159 ELSE IF( N.LT.0 ) THEN
160 INFO = -4
161 ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
162 INFO = -5
163 ELSE IF( L.LT.0 .OR. ( LEFT .AND. ( L.GT.M ) ) .OR.
164 $ ( .NOT.LEFT .AND. ( L.GT.N ) ) ) THEN
165 INFO = -6
166 ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
167 INFO = -8
168 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
169 INFO = -11
170 END IF
171 *
172 IF( INFO.EQ.0 ) THEN
173 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
174 LWKOPT = 1
175 ELSE
176 *
177 * Determine the block size. NB may be at most NBMAX, where
178 * NBMAX is used to define the local array T.
179 *
180 NB = MIN( NBMAX, ILAENV( 1, 'SORMRQ', SIDE // TRANS, M, N,
181 $ K, -1 ) )
182 LWKOPT = NW*NB
183 END IF
184 WORK( 1 ) = LWKOPT
185 *
186 IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
187 INFO = -13
188 END IF
189 END IF
190 *
191 IF( INFO.NE.0 ) THEN
192 CALL XERBLA( 'SORMRZ', -INFO )
193 RETURN
194 ELSE IF( LQUERY ) THEN
195 RETURN
196 END IF
197 *
198 * Quick return if possible
199 *
200 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
201 RETURN
202 END IF
203 *
204 NBMIN = 2
205 LDWORK = NW
206 IF( NB.GT.1 .AND. NB.LT.K ) THEN
207 IWS = NW*NB
208 IF( LWORK.LT.IWS ) THEN
209 NB = LWORK / LDWORK
210 NBMIN = MAX( 2, ILAENV( 2, 'SORMRQ', SIDE // TRANS, M, N, K,
211 $ -1 ) )
212 END IF
213 ELSE
214 IWS = NW
215 END IF
216 *
217 IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
218 *
219 * Use unblocked code
220 *
221 CALL SORMR3( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
222 $ WORK, IINFO )
223 ELSE
224 *
225 * Use blocked code
226 *
227 IF( ( LEFT .AND. .NOT.NOTRAN ) .OR.
228 $ ( .NOT.LEFT .AND. NOTRAN ) ) THEN
229 I1 = 1
230 I2 = K
231 I3 = NB
232 ELSE
233 I1 = ( ( K-1 ) / NB )*NB + 1
234 I2 = 1
235 I3 = -NB
236 END IF
237 *
238 IF( LEFT ) THEN
239 NI = N
240 JC = 1
241 JA = M - L + 1
242 ELSE
243 MI = M
244 IC = 1
245 JA = N - L + 1
246 END IF
247 *
248 IF( NOTRAN ) THEN
249 TRANST = 'T'
250 ELSE
251 TRANST = 'N'
252 END IF
253 *
254 DO 10 I = I1, I2, I3
255 IB = MIN( NB, K-I+1 )
256 *
257 * Form the triangular factor of the block reflector
258 * H = H(i+ib-1) . . . H(i+1) H(i)
259 *
260 CALL SLARZT( 'Backward', 'Rowwise', L, IB, A( I, JA ), LDA,
261 $ TAU( I ), T, LDT )
262 *
263 IF( LEFT ) THEN
264 *
265 * H or H' is applied to C(i:m,1:n)
266 *
267 MI = M - I + 1
268 IC = I
269 ELSE
270 *
271 * H or H' is applied to C(1:m,i:n)
272 *
273 NI = N - I + 1
274 JC = I
275 END IF
276 *
277 * Apply H or H'
278 *
279 CALL SLARZB( SIDE, TRANST, 'Backward', 'Rowwise', MI, NI,
280 $ IB, L, A( I, JA ), LDA, T, LDT, C( IC, JC ),
281 $ LDC, WORK, LDWORK )
282 10 CONTINUE
283 *
284 END IF
285 *
286 WORK( 1 ) = LWKOPT
287 *
288 RETURN
289 *
290 * End of SORMRZ
291 *
292 END