208 IMPLICIT NONE
209
210
211 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
212 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
213 $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
214
215 COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
216 $ * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
217 $ ALPHA( * ), BETA( * )
218
219 INTEGER, INTENT( OUT ) :: INFO
220
221
222 COMPLEX*16 CZERO, CONE
223 parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
224 $ 0.0d+0 ) )
225 DOUBLE PRECISION :: ZERO, ONE, HALF
226 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
227
228
229 INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
230 $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
231 DOUBLE PRECISION :: SAFMIN, SAFMAX, C, SCALE
232 COMPLEX*16 :: TEMP, TEMP2, TEMP3, S
233
234
236 DOUBLE PRECISION, EXTERNAL :: DLAMCH
237
238 info = 0
239 IF ( nblock_desired .LT. nshifts+1 ) THEN
240 info = -8
241 END IF
242 IF ( lwork .EQ.-1 ) THEN
243
244 work( 1 ) = n*nblock_desired
245 RETURN
246 ELSE IF ( lwork .LT. n*nblock_desired ) THEN
247 info = -25
248 END IF
249
250 IF( info.NE.0 ) THEN
251 CALL xerbla(
'ZLAQZ3', -info )
252 RETURN
253 END IF
254
255
256
257
258
259
260 safmin =
dlamch(
'SAFE MINIMUM' )
261 safmax = one/safmin
262
263 IF ( ilo .GE. ihi ) THEN
264 RETURN
265 END IF
266
267 IF ( ilschur ) THEN
268 istartm = 1
269 istopm = n
270 ELSE
271 istartm = ilo
272 istopm = ihi
273 END IF
274
275 ns = nshifts
276 npos = max( nblock_desired-ns, 1 )
277
278
279
280
281
282
283
284 CALL zlaset(
'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
285 CALL zlaset(
'FULL', ns, ns, czero, cone, zc, ldzc )
286
287 DO i = 1, ns
288
289 scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
290 IF( scale .GE. safmin .AND. scale .LE. safmax ) THEN
291 alpha( i ) = alpha( i )/scale
292 beta( i ) = beta( i )/scale
293 END IF
294
295 temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
296 temp3 = beta( i )*a( ilo+1, ilo )
297
298 IF ( abs( temp2 ) .GT. safmax .OR.
299 $ abs( temp3 ) .GT. safmax ) THEN
300 temp2 = cone
301 temp3 = czero
302 END IF
303
304 CALL zlartg( temp2, temp3, c, s, temp )
305 CALL zrot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c,
306 $ s )
307 CALL zrot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c,
308 $ s )
309 CALL zrot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
310 $ dconjg( s ) )
311
312
313 DO j = 1, ns-i
314
315 CALL zlaqz1( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
316 $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
317 $ ldqc, ns, 1, zc, ldzc )
318
319 END DO
320
321 END DO
322
323
324
325
326
327 sheight = ns+1
328 swidth = istopm-( ilo+ns )+1
329 IF ( swidth > 0 ) THEN
330 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
331 $ a( ilo, ilo+ns ), lda, czero, work, sheight )
332 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( ilo,
333 $ ilo+ns ), lda )
334 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
335 $ b( ilo, ilo+ns ), ldb, czero, work, sheight )
336 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( ilo,
337 $ ilo+ns ), ldb )
338 END IF
339 IF ( ilq ) THEN
340 CALL zgemm(
'N',
'N', n, sheight, sheight, cone, q( 1, ilo ),
341 $ ldq, qc, ldqc, czero, work, n )
342 CALL zlacpy(
'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
343 END IF
344
345
346
347 sheight = ilo-1-istartm+1
348 swidth = ns
349 IF ( sheight > 0 ) THEN
350 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
351 $ a( istartm, ilo ), lda, zc, ldzc, czero, work,
352 $ sheight )
353 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( istartm,
354 $ ilo ), lda )
355 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
356 $ b( istartm, ilo ), ldb, zc, ldzc, czero, work,
357 $ sheight )
358 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( istartm,
359 $ ilo ), ldb )
360 END IF
361 IF ( ilz ) THEN
362 CALL zgemm(
'N',
'N', n, swidth, swidth, cone, z( 1, ilo ),
363 $ ldz, zc, ldzc, czero, work, n )
364 CALL zlacpy(
'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
365 END IF
366
367
368
369
370
371 k = ilo
372 DO WHILE ( k < ihi-ns )
373 np = min( ihi-ns-k, npos )
374
375 nblock = ns+np
376
377 istartb = k+1
378
379 istopb = k+nblock-1
380
381 CALL zlaset(
'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
382 CALL zlaset(
'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
383
384
385 DO i = ns-1, 0, -1
386 DO j = 0, np-1
387
388
389
390 CALL zlaqz1( .true., .true., k+i+j, istartb, istopb, ihi,
391 $ a, lda, b, ldb, nblock, k+1, qc, ldqc,
392 $ nblock, k, zc, ldzc )
393 END DO
394 END DO
395
396
397
398
399
400
401 sheight = ns+np
402 swidth = istopm-( k+ns+np )+1
403 IF ( swidth > 0 ) THEN
404 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc,
405 $ ldqc, a( k+1, k+ns+np ), lda, czero, work,
406 $ sheight )
407 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( k+1,
408 $ k+ns+np ), lda )
409 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc,
410 $ ldqc, b( k+1, k+ns+np ), ldb, czero, work,
411 $ sheight )
412 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( k+1,
413 $ k+ns+np ), ldb )
414 END IF
415 IF ( ilq ) THEN
416 CALL zgemm(
'N',
'N', n, nblock, nblock, cone, q( 1, k+1 ),
417 $ ldq, qc, ldqc, czero, work, n )
418 CALL zlacpy(
'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
419 END IF
420
421
422
423 sheight = k-istartm+1
424 swidth = nblock
425 IF ( sheight > 0 ) THEN
426 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
427 $ a( istartm, k ), lda, zc, ldzc, czero, work,
428 $ sheight )
429 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
430 $ a( istartm, k ), lda )
431 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
432 $ b( istartm, k ), ldb, zc, ldzc, czero, work,
433 $ sheight )
434 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
435 $ b( istartm, k ), ldb )
436 END IF
437 IF ( ilz ) THEN
438 CALL zgemm(
'N',
'N', n, nblock, nblock, cone, z( 1, k ),
439 $ ldz, zc, ldzc, czero, work, n )
440 CALL zlacpy(
'ALL', n, nblock, work, n, z( 1, k ), ldz )
441 END IF
442
443 k = k+np
444
445 END DO
446
447
448
449
450 CALL zlaset(
'FULL', ns, ns, czero, cone, qc, ldqc )
451 CALL zlaset(
'FULL', ns+1, ns+1, czero, cone, zc, ldzc )
452
453
454 istartb = ihi-ns+1
455
456 istopb = ihi
457
458 DO i = 1, ns
459
460 DO ishift = ihi-i, ihi-1
461 CALL zlaqz1( .true., .true., ishift, istartb, istopb, ihi,
462 $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
463 $ ihi-ns, zc, ldzc )
464 END DO
465
466 END DO
467
468
469
470
471
472 sheight = ns
473 swidth = istopm-( ihi+1 )+1
474 IF ( swidth > 0 ) THEN
475 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
476 $ a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
477 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
478 $ a( ihi-ns+1, ihi+1 ), lda )
479 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
480 $ b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
481 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
482 $ b( ihi-ns+1, ihi+1 ), ldb )
483 END IF
484 IF ( ilq ) THEN
485 CALL zgemm(
'N',
'N', n, ns, ns, cone, q( 1, ihi-ns+1 ), ldq,
486 $ qc, ldqc, czero, work, n )
487 CALL zlacpy(
'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
488 END IF
489
490
491
492 sheight = ihi-ns-istartm+1
493 swidth = ns+1
494 IF ( sheight > 0 ) THEN
495 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
496 $ a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
497 $ sheight )
498 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( istartm,
499 $ ihi-ns ), lda )
500 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
501 $ b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
502 $ sheight )
503 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( istartm,
504 $ ihi-ns ), ldb )
505 END IF
506 IF ( ilz ) THEN
507 CALL zgemm(
'N',
'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ), ldz,
508 $ zc, ldzc, czero, work, n )
509 CALL zlacpy(
'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
510 END IF
511
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.
double precision function dlamch(cmach)
DLAMCH
subroutine zlaqz1(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
ZLAQZ1
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.