216
217
218
219
220
221
222 CHARACTER JOBVL, JOBVR
223 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
224
225
226 REAL RWORK( * )
227 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
228 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
229 $ WORK( * )
230
231
232
233
234
235 REAL ZERO, ONE
236 parameter( zero = 0.0e0, one = 1.0e0 )
237 COMPLEX CZERO, CONE
238 parameter( czero = ( 0.0e0, 0.0e0 ),
239 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
248 $ SMLNUM, TEMP
249 COMPLEX X
250
251
252 LOGICAL LDUMMA( 1 )
253
254
258
259
260 LOGICAL LSAME
261 REAL CLANGE, SLAMCH
263
264
265 INTRINSIC abs, aimag, max, real, sqrt
266
267
268 REAL ABS1
269
270
271 abs1( x ) = abs( real( x ) ) + abs( aimag( 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 cgeqrf( n, n, b, ldb, work, work, -1, ierr )
326 lwkopt = max( n, n+int( work( 1 ) ) )
327 CALL cunmqr(
'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 cungqr( 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 cgghd3( 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 claqz0(
'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 cgghd3(
'N',
'N', n, 1, n, a, lda, b, ldb, vl, ldvl,
344 $ vr, ldvr, work, -1, ierr )
345 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
346 CALL claqz0(
'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 ) = cmplx( lwkopt )
352 END IF
353
354 IF( info.NE.0 ) THEN
355 CALL xerbla(
'CGGEV3 ', -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 slabad( smlnum, bignum )
372 smlnum = sqrt( smlnum ) / eps
373 bignum = one / smlnum
374
375
376
377 anrm =
clange(
'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 clascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
388
389
390
391 bnrm =
clange(
'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 clascl(
'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 cggbal(
'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 cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
422 $ work( iwrk ), lwork+1-iwrk, ierr )
423
424
425
426 CALL cunmqr(
'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 claset(
'Full', n, n, czero, cone, vl, ldvl )
434 IF( irows.GT.1 ) THEN
435 CALL clacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
436 $ vl( ilo+1, ilo ), ldvl )
437 END IF
438 CALL cungqr( 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 claset(
'Full', n, n, czero, cone, vr, ldvr )
446
447
448
449 IF( ilv ) THEN
450
451
452
453 CALL cgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
454 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk,
455 $ ierr )
456 ELSE
457 CALL cgghd3(
'N',
'N', irows, 1, irows, a( ilo, ilo ), lda,
458 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
459 $ work( iwrk ), lwork+1-iwrk, ierr )
460 END IF
461
462
463
464
465 iwrk = itau
466 IF( ilv ) THEN
467 chtemp = 'S'
468 ELSE
469 chtemp = 'E'
470 END IF
471 CALL claqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
472 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
473 $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
474 IF( ierr.NE.0 ) THEN
475 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
476 info = ierr
477 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
478 info = ierr - n
479 ELSE
480 info = n + 1
481 END IF
482 GO TO 70
483 END IF
484
485
486
487 IF( ilv ) THEN
488 IF( ilvl ) THEN
489 IF( ilvr ) THEN
490 chtemp = 'B'
491 ELSE
492 chtemp = 'L'
493 END IF
494 ELSE
495 chtemp = 'R'
496 END IF
497
498 CALL ctgevc( chtemp,
'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
499 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
500 $ ierr )
501 IF( ierr.NE.0 ) THEN
502 info = n + 2
503 GO TO 70
504 END IF
505
506
507
508 IF( ilvl ) THEN
509 CALL cggbak(
'P',
'L', n, ilo, ihi, rwork( ileft ),
510 $ rwork( iright ), n, vl, ldvl, ierr )
511 DO 30 jc = 1, n
512 temp = zero
513 DO 10 jr = 1, n
514 temp = max( temp, abs1( vl( jr, jc ) ) )
515 10 CONTINUE
516 IF( temp.LT.smlnum )
517 $ GO TO 30
518 temp = one / temp
519 DO 20 jr = 1, n
520 vl( jr, jc ) = vl( jr, jc )*temp
521 20 CONTINUE
522 30 CONTINUE
523 END IF
524 IF( ilvr ) THEN
525 CALL cggbak(
'P',
'R', n, ilo, ihi, rwork( ileft ),
526 $ rwork( iright ), n, vr, ldvr, ierr )
527 DO 60 jc = 1, n
528 temp = zero
529 DO 40 jr = 1, n
530 temp = max( temp, abs1( vr( jr, jc ) ) )
531 40 CONTINUE
532 IF( temp.LT.smlnum )
533 $ GO TO 60
534 temp = one / temp
535 DO 50 jr = 1, n
536 vr( jr, jc ) = vr( jr, jc )*temp
537 50 CONTINUE
538 60 CONTINUE
539 END IF
540 END IF
541
542
543
544 70 CONTINUE
545
546 IF( ilascl )
547 $
CALL clascl(
'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
548
549 IF( ilbscl )
550 $
CALL clascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
551
552 work( 1 ) = cmplx( lwkopt )
553 RETURN
554
555
556
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine cggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
CGGBAL
subroutine cggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
CGGBAK
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC
recursive subroutine claqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, INFO)
CLAQZ0
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
CGGHD3
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
real function slamch(CMACH)
SLAMCH