160
161
162
163
164
165 IMPLICIT NONE
166
167
168 CHARACTER SIDE, TRANS
169 INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
170
171
172 COMPLEX*16 Q( LDQ, * ), C( LDC, * ), WORK( * )
173
174
175
176
177
178 COMPLEX*16 ONE
179 parameter( one = ( 1.0d+0, 0.0d+0 ) )
180
181
182 LOGICAL LEFT, LQUERY, NOTRAN
183 INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
184
185
186 LOGICAL LSAME
188
189
191
192
193 INTRINSIC dcmplx, max, min
194
195
196
197
198
199 info = 0
200 left =
lsame( side,
'L' )
201 notran =
lsame( trans,
'N' )
202 lquery = ( lwork.EQ.-1 )
203
204
205
206
207 IF( left ) THEN
208 nq = m
209 ELSE
210 nq = n
211 END IF
212 nw = nq
213 IF( n1.EQ.0 .OR. n2.EQ.0 ) nw = 1
214 IF( .NOT.left .AND. .NOT.
lsame( side,
'R' ) )
THEN
215 info = -1
216 ELSE IF( .NOT.
lsame( trans,
'N' ) .AND.
217 $ .NOT.
lsame( trans,
'C' ) )
218 $ THEN
219 info = -2
220 ELSE IF( m.LT.0 ) THEN
221 info = -3
222 ELSE IF( n.LT.0 ) THEN
223 info = -4
224 ELSE IF( n1.LT.0 .OR. n1+n2.NE.nq ) THEN
225 info = -5
226 ELSE IF( n2.LT.0 ) THEN
227 info = -6
228 ELSE IF( ldq.LT.max( 1, nq ) ) THEN
229 info = -8
230 ELSE IF( ldc.LT.max( 1, m ) ) THEN
231 info = -10
232 ELSE IF( lwork.LT.nw .AND. .NOT.lquery ) THEN
233 info = -12
234 END IF
235
236 IF( info.EQ.0 ) THEN
237 lwkopt = m*n
238 work( 1 ) = dcmplx( lwkopt )
239 END IF
240
241 IF( info.NE.0 ) THEN
242 CALL xerbla(
'ZUNM22', -info )
243 RETURN
244 ELSE IF( lquery ) THEN
245 RETURN
246 END IF
247
248
249
250 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
251 work( 1 ) = 1
252 RETURN
253 END IF
254
255
256
257 IF( n1.EQ.0 ) THEN
258 CALL ztrmm( side,
'Upper', trans,
'Non-Unit', m, n, one,
259 $ q, ldq, c, ldc )
260 work( 1 ) = one
261 RETURN
262 ELSE IF( n2.EQ.0 ) THEN
263 CALL ztrmm( side,
'Lower', trans,
'Non-Unit', m, n, one,
264 $ q, ldq, c, ldc )
265 work( 1 ) = one
266 RETURN
267 END IF
268
269
270
271 nb = max( 1, min( lwork, lwkopt ) / nq )
272
273 IF( left ) THEN
274 IF( notran ) THEN
275 DO i = 1, n, nb
276 len = min( nb, n-i+1 )
277 ldwork = m
278
279
280
281 CALL zlacpy(
'All', n1, len, c( n2+1, i ), ldc, work,
282 $ ldwork )
283 CALL ztrmm(
'Left',
'Lower',
'No Transpose',
284 $ 'Non-Unit',
285 $ n1, len, one, q( 1, n2+1 ), ldq, work,
286 $ ldwork )
287
288
289
290 CALL zgemm(
'No Transpose',
'No Transpose', n1, len,
291 $ n2,
292 $ one, q, ldq, c( 1, i ), ldc, one, work,
293 $ ldwork )
294
295
296
297 CALL zlacpy(
'All', n2, len, c( 1, i ), ldc,
298 $ work( n1+1 ), ldwork )
299 CALL ztrmm(
'Left',
'Upper',
'No Transpose',
300 $ 'Non-Unit',
301 $ n2, len, one, q( n1+1, 1 ), ldq,
302 $ work( n1+1 ), ldwork )
303
304
305
306 CALL zgemm(
'No Transpose',
'No Transpose', n2, len,
307 $ n1,
308 $ one, q( n1+1, n2+1 ), ldq, c( n2+1, i ), ldc,
309 $ one, work( n1+1 ), ldwork )
310
311
312
313 CALL zlacpy(
'All', m, len, work, ldwork, c( 1, i ),
314 $ ldc )
315 END DO
316 ELSE
317 DO i = 1, n, nb
318 len = min( nb, n-i+1 )
319 ldwork = m
320
321
322
323 CALL zlacpy(
'All', n2, len, c( n1+1, i ), ldc, work,
324 $ ldwork )
325 CALL ztrmm(
'Left',
'Upper',
'Conjugate',
'Non-Unit',
326 $ n2, len, one, q( n1+1, 1 ), ldq, work,
327 $ ldwork )
328
329
330
331 CALL zgemm(
'Conjugate',
'No Transpose', n2, len, n1,
332 $ one, q, ldq, c( 1, i ), ldc, one, work,
333 $ ldwork )
334
335
336
337 CALL zlacpy(
'All', n1, len, c( 1, i ), ldc,
338 $ work( n2+1 ), ldwork )
339 CALL ztrmm(
'Left',
'Lower',
'Conjugate',
'Non-Unit',
340 $ n1, len, one, q( 1, n2+1 ), ldq,
341 $ work( n2+1 ), ldwork )
342
343
344
345 CALL zgemm(
'Conjugate',
'No Transpose', n1, len, n2,
346 $ one, q( n1+1, n2+1 ), ldq, c( n1+1, i ), ldc,
347 $ one, work( n2+1 ), ldwork )
348
349
350
351 CALL zlacpy(
'All', m, len, work, ldwork, c( 1, i ),
352 $ ldc )
353 END DO
354 END IF
355 ELSE
356 IF( notran ) THEN
357 DO i = 1, m, nb
358 len = min( nb, m-i+1 )
359 ldwork = len
360
361
362
363 CALL zlacpy(
'All', len, n2, c( i, n1+1 ), ldc, work,
364 $ ldwork )
365 CALL ztrmm(
'Right',
'Upper',
'No Transpose',
366 $ 'Non-Unit',
367 $ len, n2, one, q( n1+1, 1 ), ldq, work,
368 $ ldwork )
369
370
371
372 CALL zgemm(
'No Transpose',
'No Transpose', len, n2,
373 $ n1,
374 $ one, c( i, 1 ), ldc, q, ldq, one, work,
375 $ ldwork )
376
377
378
379 CALL zlacpy(
'All', len, n1, c( i, 1 ), ldc,
380 $ work( 1 + n2*ldwork ), ldwork )
381 CALL ztrmm(
'Right',
'Lower',
'No Transpose',
382 $ 'Non-Unit',
383 $ len, n1, one, q( 1, n2+1 ), ldq,
384 $ work( 1 + n2*ldwork ), ldwork )
385
386
387
388 CALL zgemm(
'No Transpose',
'No Transpose', len, n1,
389 $ n2,
390 $ one, c( i, n1+1 ), ldc, q( n1+1, n2+1 ), ldq,
391 $ one, work( 1 + n2*ldwork ), ldwork )
392
393
394
395 CALL zlacpy(
'All', len, n, work, ldwork, c( i, 1 ),
396 $ ldc )
397 END DO
398 ELSE
399 DO i = 1, m, nb
400 len = min( nb, m-i+1 )
401 ldwork = len
402
403
404
405 CALL zlacpy(
'All', len, n1, c( i, n2+1 ), ldc, work,
406 $ ldwork )
407 CALL ztrmm(
'Right',
'Lower',
'Conjugate',
'Non-Unit',
408 $ len, n1, one, q( 1, n2+1 ), ldq, work,
409 $ ldwork )
410
411
412
413 CALL zgemm(
'No Transpose',
'Conjugate', len, n1, n2,
414 $ one, c( i, 1 ), ldc, q, ldq, one, work,
415 $ ldwork )
416
417
418
419 CALL zlacpy(
'All', len, n2, c( i, 1 ), ldc,
420 $ work( 1 + n1*ldwork ), ldwork )
421 CALL ztrmm(
'Right',
'Upper',
'Conjugate',
'Non-Unit',
422 $ len, n2, one, q( n1+1, 1 ), ldq,
423 $ work( 1 + n1*ldwork ), ldwork )
424
425
426
427 CALL zgemm(
'No Transpose',
'Conjugate', len, n2, n1,
428 $ one, c( i, n2+1 ), ldc, q( n1+1, n2+1 ), ldq,
429 $ one, work( 1 + n1*ldwork ), ldwork )
430
431
432
433 CALL zlacpy(
'All', len, n, work, ldwork, c( i, 1 ),
434 $ ldc )
435 END DO
436 END IF
437 END IF
438
439 work( 1 ) = dcmplx( lwkopt )
440 RETURN
441
442
443
subroutine xerbla(srname, info)
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
logical function lsame(ca, cb)
LSAME
subroutine ztrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRMM