216
217
218
219
220
221
222 CHARACTER JOBVL, JOBVR
223 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
224
225
226 DOUBLE PRECISION RWORK( * )
227 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
228 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
229 $ WORK( * )
230
231
232
233
234
235 DOUBLE PRECISION ZERO, ONE
236 parameter( zero = 0.0d0, one = 1.0d0 )
237 COMPLEX*16 CZERO, CONE
238 parameter( czero = ( 0.0d0, 0.0d0 ),
239 $ cone = ( 1.0d0, 0.0d0 ) )
240
241
242 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
243 CHARACTER CHTEMP
244 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
245 $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
246 $ LWKOPT
247 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
248 $ SMLNUM, TEMP
249 COMPLEX*16 X
250
251
252 LOGICAL LDUMMA( 1 )
253
254
258
259
260 LOGICAL LSAME
261 DOUBLE PRECISION DLAMCH, ZLANGE
263
264
265 INTRINSIC abs, dble, dimag, max, sqrt
266
267
268 DOUBLE PRECISION ABS1
269
270
271 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
272
273
274
275
276
277 IF(
lsame( jobvl,
'N' ) )
THEN
278 ijobvl = 1
279 ilvl = .false.
280 ELSE IF(
lsame( jobvl,
'V' ) )
THEN
281 ijobvl = 2
282 ilvl = .true.
283 ELSE
284 ijobvl = -1
285 ilvl = .false.
286 END IF
287
288 IF(
lsame( jobvr,
'N' ) )
THEN
289 ijobvr = 1
290 ilvr = .false.
291 ELSE IF(
lsame( jobvr,
'V' ) )
THEN
292 ijobvr = 2
293 ilvr = .true.
294 ELSE
295 ijobvr = -1
296 ilvr = .false.
297 END IF
298 ilv = ilvl .OR. ilvr
299
300
301
302 info = 0
303 lquery = ( lwork.EQ.-1 )
304 IF( ijobvl.LE.0 ) THEN
305 info = -1
306 ELSE IF( ijobvr.LE.0 ) THEN
307 info = -2
308 ELSE IF( n.LT.0 ) THEN
309 info = -3
310 ELSE IF( lda.LT.max( 1, n ) ) THEN
311 info = -5
312 ELSE IF( ldb.LT.max( 1, n ) ) THEN
313 info = -7
314 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
315 info = -11
316 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
317 info = -13
318 ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
319 info = -15
320 END IF
321
322
323
324 IF( info.EQ.0 ) THEN
325 CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
326 lwkopt = max( 1, n+int( work( 1 ) ) )
327 CALL zunmqr(
'L',
'C', n, n, n, b, ldb, work, a, lda, work,
328 $ -1, ierr )
329 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
330 IF( ilvl ) THEN
331 CALL zungqr( n, n, n, vl, ldvl, work, work, -1, ierr )
332 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
333 END IF
334 IF( ilv ) THEN
335 CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
336 $ ldvl, vr, ldvr, work, -1, ierr )
337 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
338 CALL zlaqz0(
'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
339 $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
340 $ rwork, 0, ierr )
341 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
342 ELSE
343 CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
344 $ ldvl, vr, ldvr, work, -1, ierr )
345 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
346 CALL zlaqz0(
'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
347 $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
348 $ rwork, 0, ierr )
349 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
350 END IF
351 work( 1 ) = dcmplx( lwkopt )
352 END IF
353
354 IF( info.NE.0 ) THEN
355 CALL xerbla(
'ZGGEV3 ', -info )
356 RETURN
357 ELSE IF( lquery ) THEN
358 RETURN
359 END IF
360
361
362
363 IF( n.EQ.0 )
364 $ RETURN
365
366
367
370 bignum = one / smlnum
371 CALL dlabad( smlnum, bignum )
372 smlnum = sqrt( smlnum ) / eps
373 bignum = one / smlnum
374
375
376
377 anrm =
zlange(
'M', n, n, a, lda, rwork )
378 ilascl = .false.
379 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
380 anrmto = smlnum
381 ilascl = .true.
382 ELSE IF( anrm.GT.bignum ) THEN
383 anrmto = bignum
384 ilascl = .true.
385 END IF
386 IF( ilascl )
387 $
CALL zlascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
388
389
390
391 bnrm =
zlange(
'M', n, n, b, ldb, rwork )
392 ilbscl = .false.
393 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
394 bnrmto = smlnum
395 ilbscl = .true.
396 ELSE IF( bnrm.GT.bignum ) THEN
397 bnrmto = bignum
398 ilbscl = .true.
399 END IF
400 IF( ilbscl )
401 $
CALL zlascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
402
403
404
405 ileft = 1
406 iright = n + 1
407 irwrk = iright + n
408 CALL zggbal(
'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
409 $ rwork( iright ), rwork( irwrk ), ierr )
410
411
412
413 irows = ihi + 1 - ilo
414 IF( ilv ) THEN
415 icols = n + 1 - ilo
416 ELSE
417 icols = irows
418 END IF
419 itau = 1
420 iwrk = itau + irows
421 CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
422 $ work( iwrk ), lwork+1-iwrk, ierr )
423
424
425
426 CALL zunmqr(
'L',
'C', irows, icols, irows, b( ilo, ilo ), ldb,
427 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
428 $ lwork+1-iwrk, ierr )
429
430
431
432 IF( ilvl ) THEN
433 CALL zlaset(
'Full', n, n, czero, cone, vl, ldvl )
434 IF( irows.GT.1 ) THEN
435 CALL zlacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
436 $ vl( ilo+1, ilo ), ldvl )
437 END IF
438 CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
439 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
440 END IF
441
442
443
444 IF( ilvr )
445 $
CALL zlaset(
'Full', n, n, czero, cone, vr, ldvr )
446
447
448
449 IF( ilv ) THEN
450
451
452
453 CALL zgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
454 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
455 ELSE
456 CALL zgghd3(
'N',
'N', irows, 1, irows, a( ilo, ilo ), lda,
457 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
458 $ work( iwrk ), lwork+1-iwrk, ierr )
459 END IF
460
461
462
463
464 iwrk = itau
465 IF( ilv ) THEN
466 chtemp = 'S'
467 ELSE
468 chtemp = 'E'
469 END IF
470 CALL zlaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
471 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
472 $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
473 IF( ierr.NE.0 ) THEN
474 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
475 info = ierr
476 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
477 info = ierr - n
478 ELSE
479 info = n + 1
480 END IF
481 GO TO 70
482 END IF
483
484
485
486 IF( ilv ) THEN
487 IF( ilvl ) THEN
488 IF( ilvr ) THEN
489 chtemp = 'B'
490 ELSE
491 chtemp = 'L'
492 END IF
493 ELSE
494 chtemp = 'R'
495 END IF
496
497 CALL ztgevc( chtemp,
'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
498 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
499 $ ierr )
500 IF( ierr.NE.0 ) THEN
501 info = n + 2
502 GO TO 70
503 END IF
504
505
506
507 IF( ilvl ) THEN
508 CALL zggbak(
'P',
'L', n, ilo, ihi, rwork( ileft ),
509 $ rwork( iright ), n, vl, ldvl, ierr )
510 DO 30 jc = 1, n
511 temp = zero
512 DO 10 jr = 1, n
513 temp = max( temp, abs1( vl( jr, jc ) ) )
514 10 CONTINUE
515 IF( temp.LT.smlnum )
516 $ GO TO 30
517 temp = one / temp
518 DO 20 jr = 1, n
519 vl( jr, jc ) = vl( jr, jc )*temp
520 20 CONTINUE
521 30 CONTINUE
522 END IF
523 IF( ilvr ) THEN
524 CALL zggbak(
'P',
'R', n, ilo, ihi, rwork( ileft ),
525 $ rwork( iright ), n, vr, ldvr, ierr )
526 DO 60 jc = 1, n
527 temp = zero
528 DO 40 jr = 1, n
529 temp = max( temp, abs1( vr( jr, jc ) ) )
530 40 CONTINUE
531 IF( temp.LT.smlnum )
532 $ GO TO 60
533 temp = one / temp
534 DO 50 jr = 1, n
535 vr( jr, jc ) = vr( jr, jc )*temp
536 50 CONTINUE
537 60 CONTINUE
538 END IF
539 END IF
540
541
542
543 70 CONTINUE
544
545 IF( ilascl )
546 $
CALL zlascl(
'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
547
548 IF( ilbscl )
549 $
CALL zlascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
550
551 work( 1 ) = dcmplx( lwkopt )
552 RETURN
553
554
555
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
recursive subroutine zlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, INFO)
ZLAQZ0
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
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 zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
ZGGHD3
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.