Mercurial > octave

comparison libcruft/lapack/sormbr.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 SORMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C,
2 $ LDC, WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER SIDE, TRANS, VECT
10 INTEGER INFO, K, LDA, LDC, LWORK, M, N
11 * ..
12 * .. Array Arguments ..
13 REAL A( LDA, * ), C( LDC, * ), TAU( * ),
14 $ WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * If VECT = 'Q', SORMBR overwrites the general real M-by-N matrix C
21 * with
22 * SIDE = 'L' SIDE = 'R'
23 * TRANS = 'N': Q * C C * Q
24 * TRANS = 'T': Q**T * C C * Q**T
25 *
26 * If VECT = 'P', SORMBR overwrites the general real M-by-N matrix C
27 * with
28 * SIDE = 'L' SIDE = 'R'
29 * TRANS = 'N': P * C C * P
30 * TRANS = 'T': P**T * C C * P**T
31 *
32 * Here Q and P**T are the orthogonal matrices determined by SGEBRD when
33 * reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and
34 * P**T are defined as products of elementary reflectors H(i) and G(i)
35 * respectively.
36 *
37 * Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the
38 * order of the orthogonal matrix Q or P**T that is applied.
39 *
40 * If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:
41 * if nq >= k, Q = H(1) H(2) . . . H(k);
42 * if nq < k, Q = H(1) H(2) . . . H(nq-1).
43 *
44 * If VECT = 'P', A is assumed to have been a K-by-NQ matrix:
45 * if k < nq, P = G(1) G(2) . . . G(k);
46 * if k >= nq, P = G(1) G(2) . . . G(nq-1).
47 *
48 * Arguments
49 * =========
50 *
51 * VECT (input) CHARACTER*1
52 * = 'Q': apply Q or Q**T;
53 * = 'P': apply P or P**T.
54 *
55 * SIDE (input) CHARACTER*1
56 * = 'L': apply Q, Q**T, P or P**T from the Left;
57 * = 'R': apply Q, Q**T, P or P**T from the Right.
58 *
59 * TRANS (input) CHARACTER*1
60 * = 'N': No transpose, apply Q or P;
61 * = 'T': Transpose, apply Q**T or P**T.
62 *
63 * M (input) INTEGER
64 * The number of rows of the matrix C. M >= 0.
65 *
66 * N (input) INTEGER
67 * The number of columns of the matrix C. N >= 0.
68 *
69 * K (input) INTEGER
70 * If VECT = 'Q', the number of columns in the original
71 * matrix reduced by SGEBRD.
72 * If VECT = 'P', the number of rows in the original
73 * matrix reduced by SGEBRD.
74 * K >= 0.
75 *
76 * A (input) REAL array, dimension
77 * (LDA,min(nq,K)) if VECT = 'Q'
78 * (LDA,nq) if VECT = 'P'
79 * The vectors which define the elementary reflectors H(i) and
80 * G(i), whose products determine the matrices Q and P, as
81 * returned by SGEBRD.
82 *
83 * LDA (input) INTEGER
84 * The leading dimension of the array A.
85 * If VECT = 'Q', LDA >= max(1,nq);
86 * if VECT = 'P', LDA >= max(1,min(nq,K)).
87 *
88 * TAU (input) REAL array, dimension (min(nq,K))
89 * TAU(i) must contain the scalar factor of the elementary
90 * reflector H(i) or G(i) which determines Q or P, as returned
91 * by SGEBRD in the array argument TAUQ or TAUP.
92 *
93 * C (input/output) REAL array, dimension (LDC,N)
94 * On entry, the M-by-N matrix C.
95 * On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q
96 * or P*C or P**T*C or C*P or C*P**T.
97 *
98 * LDC (input) INTEGER
99 * The leading dimension of the array C. LDC >= max(1,M).
100 *
101 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
102 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
103 *
104 * LWORK (input) INTEGER
105 * The dimension of the array WORK.
106 * If SIDE = 'L', LWORK >= max(1,N);
107 * if SIDE = 'R', LWORK >= max(1,M).
108 * For optimum performance LWORK >= N*NB if SIDE = 'L', and
109 * LWORK >= M*NB if SIDE = 'R', where NB is the optimal
110 * blocksize.
111 *
112 * If LWORK = -1, then a workspace query is assumed; the routine
113 * only calculates the optimal size of the WORK array, returns
114 * this value as the first entry of the WORK array, and no error
115 * message related to LWORK is issued by XERBLA.
116 *
117 * INFO (output) INTEGER
118 * = 0: successful exit
119 * < 0: if INFO = -i, the i-th argument had an illegal value
120 *
121 * =====================================================================
122 *
123 * .. Local Scalars ..
124 LOGICAL APPLYQ, LEFT, LQUERY, NOTRAN
125 CHARACTER TRANST
126 INTEGER I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW
127 * ..
128 * .. External Functions ..
129 LOGICAL LSAME
130 INTEGER ILAENV
131 EXTERNAL ILAENV, LSAME
132 * ..
133 * .. External Subroutines ..
134 EXTERNAL SORMLQ, SORMQR, XERBLA
135 * ..
136 * .. Intrinsic Functions ..
137 INTRINSIC MAX, MIN
138 * ..
139 * .. Executable Statements ..
140 *
141 * Test the input arguments
142 *
143 INFO = 0
144 APPLYQ = LSAME( VECT, 'Q' )
145 LEFT = LSAME( SIDE, 'L' )
146 NOTRAN = LSAME( TRANS, 'N' )
147 LQUERY = ( LWORK.EQ.-1 )
148 *
149 * NQ is the order of Q or P and NW is the minimum dimension of WORK
150 *
151 IF( LEFT ) THEN
152 NQ = M
153 NW = N
154 ELSE
155 NQ = N
156 NW = M
157 END IF
158 IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
159 INFO = -1
160 ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
161 INFO = -2
162 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
163 INFO = -3
164 ELSE IF( M.LT.0 ) THEN
165 INFO = -4
166 ELSE IF( N.LT.0 ) THEN
167 INFO = -5
168 ELSE IF( K.LT.0 ) THEN
169 INFO = -6
170 ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR.
171 $ ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) )
172 $ THEN
173 INFO = -8
174 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
175 INFO = -11
176 ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
177 INFO = -13
178 END IF
179 *
180 IF( INFO.EQ.0 ) THEN
181 IF( APPLYQ ) THEN
182 IF( LEFT ) THEN
183 NB = ILAENV( 1, 'SORMQR', SIDE // TRANS, M-1, N, M-1,
184 $ -1 )
185 ELSE
186 NB = ILAENV( 1, 'SORMQR', SIDE // TRANS, M, N-1, N-1,
187 $ -1 )
188 END IF
189 ELSE
190 IF( LEFT ) THEN
191 NB = ILAENV( 1, 'SORMLQ', SIDE // TRANS, M-1, N, M-1,
192 $ -1 )
193 ELSE
194 NB = ILAENV( 1, 'SORMLQ', SIDE // TRANS, M, N-1, N-1,
195 $ -1 )
196 END IF
197 END IF
198 LWKOPT = MAX( 1, NW )*NB
199 WORK( 1 ) = LWKOPT
200 END IF
201 *
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'SORMBR', -INFO )
204 RETURN
205 ELSE IF( LQUERY ) THEN
206 RETURN
207 END IF
208 *
209 * Quick return if possible
210 *
211 WORK( 1 ) = 1
212 IF( M.EQ.0 .OR. N.EQ.0 )
213 $ RETURN
214 *
215 IF( APPLYQ ) THEN
216 *
217 * Apply Q
218 *
219 IF( NQ.GE.K ) THEN
220 *
221 * Q was determined by a call to SGEBRD with nq >= k
222 *
223 CALL SORMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
224 $ WORK, LWORK, IINFO )
225 ELSE IF( NQ.GT.1 ) THEN
226 *
227 * Q was determined by a call to SGEBRD with nq < k
228 *
229 IF( LEFT ) THEN
230 MI = M - 1
231 NI = N
232 I1 = 2
233 I2 = 1
234 ELSE
235 MI = M
236 NI = N - 1
237 I1 = 1
238 I2 = 2
239 END IF
240 CALL SORMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU,
241 $ C( I1, I2 ), LDC, WORK, LWORK, IINFO )
242 END IF
243 ELSE
244 *
245 * Apply P
246 *
247 IF( NOTRAN ) THEN
248 TRANST = 'T'
249 ELSE
250 TRANST = 'N'
251 END IF
252 IF( NQ.GT.K ) THEN
253 *
254 * P was determined by a call to SGEBRD with nq > k
255 *
256 CALL SORMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC,
257 $ WORK, LWORK, IINFO )
258 ELSE IF( NQ.GT.1 ) THEN
259 *
260 * P was determined by a call to SGEBRD with nq <= k
261 *
262 IF( LEFT ) THEN
263 MI = M - 1
264 NI = N
265 I1 = 2
266 I2 = 1
267 ELSE
268 MI = M
269 NI = N - 1
270 I1 = 1
271 I2 = 2
272 END IF
273 CALL SORMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA,
274 $ TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO )
275 END IF
276 END IF
277 WORK( 1 ) = LWKOPT
278 RETURN
279 *
280 * End of SORMBR
281 *
282 END