179 implicit none
180
181
182
183
184
185
186 CHARACTER JOBVL, JOBVR
187 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
188
189
190 REAL RWORK( * )
191 COMPLEX A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
192 $ W( * ), WORK( * )
193
194
195
196
197
198 REAL ZERO, ONE
199 parameter( zero = 0.0e0, one = 1.0e0 )
200
201
202 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
203 CHARACTER SIDE
204 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
205 $ IWRK, K, LWORK_TREVC, MAXWRK, MINWRK, NOUT
206 REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
207 COMPLEX TMP
208
209
210 LOGICAL SELECT( 1 )
211 REAL DUM( 1 )
212
213
218
219
220 LOGICAL LSAME
221 INTEGER ISAMAX, ILAENV
222 REAL SLAMCH, SCNRM2, CLANGE,
223 $ SROUNDUP_LWORK
227
228
229 INTRINSIC real, cmplx, conjg, aimag, max, sqrt
230
231
232
233
234
235 info = 0
236 lquery = ( lwork.EQ.-1 )
237 wantvl =
lsame( jobvl,
'V' )
238 wantvr =
lsame( jobvr,
'V' )
239 IF( ( .NOT.wantvl ) .AND. ( .NOT.
lsame( jobvl,
'N' ) ) )
THEN
240 info = -1
241 ELSE IF( ( .NOT.wantvr ) .AND.
242 $ ( .NOT.
lsame( jobvr,
'N' ) ) )
THEN
243 info = -2
244 ELSE IF( n.LT.0 ) THEN
245 info = -3
246 ELSE IF( lda.LT.max( 1, n ) ) THEN
247 info = -5
248 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
249 info = -8
250 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
251 info = -10
252 END IF
253
254
255
256
257
258
259
260
261
262
263
264
265 IF( info.EQ.0 ) THEN
266 IF( n.EQ.0 ) THEN
267 minwrk = 1
268 maxwrk = 1
269 ELSE
270 maxwrk = n + n*
ilaenv( 1,
'CGEHRD',
' ', n, 1, n, 0 )
271 minwrk = 2*n
272 IF( wantvl ) THEN
273 maxwrk = max( maxwrk, n + ( n - 1 )*
ilaenv( 1,
274 $ 'CUNGHR',
275 $ ' ', n, 1, n, -1 ) )
276 CALL ctrevc3(
'L',
'B',
SELECT, n, a, lda,
277 $ vl, ldvl, vr, ldvr,
278 $ n, nout, work, -1, rwork, -1, ierr )
279 lwork_trevc = int( work(1) )
280 maxwrk = max( maxwrk, n + lwork_trevc )
281 CALL chseqr(
'S',
'V', n, 1, n, a, lda, w, vl, ldvl,
282 $ work, -1, info )
283 ELSE IF( wantvr ) THEN
284 maxwrk = max( maxwrk, n + ( n - 1 )*
ilaenv( 1,
285 $ 'CUNGHR',
286 $ ' ', n, 1, n, -1 ) )
287 CALL ctrevc3(
'R',
'B',
SELECT, n, a, lda,
288 $ vl, ldvl, vr, ldvr,
289 $ n, nout, work, -1, rwork, -1, ierr )
290 lwork_trevc = int( work(1) )
291 maxwrk = max( maxwrk, n + lwork_trevc )
292 CALL chseqr(
'S',
'V', n, 1, n, a, lda, w, vr, ldvr,
293 $ work, -1, info )
294 ELSE
295 CALL chseqr(
'E',
'N', n, 1, n, a, lda, w, vr, ldvr,
296 $ work, -1, info )
297 END IF
298 hswork = int( work(1) )
299 maxwrk = max( maxwrk, hswork, minwrk )
300 END IF
302
303 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
304 info = -12
305 END IF
306 END IF
307
308 IF( info.NE.0 ) THEN
309 CALL xerbla(
'CGEEV ', -info )
310 RETURN
311 ELSE IF( lquery ) THEN
312 RETURN
313 END IF
314
315
316
317 IF( n.EQ.0 )
318 $ RETURN
319
320
321
324 bignum = one / smlnum
325 smlnum = sqrt( smlnum ) / eps
326 bignum = one / smlnum
327
328
329
330 anrm =
clange(
'M', n, n, a, lda, dum )
331 scalea = .false.
332 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
333 scalea = .true.
334 cscale = smlnum
335 ELSE IF( anrm.GT.bignum ) THEN
336 scalea = .true.
337 cscale = bignum
338 END IF
339 IF( scalea )
340 $
CALL clascl(
'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
341
342
343
344
345
346 ibal = 1
347 CALL cgebal(
'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
348
349
350
351
352
353 itau = 1
354 iwrk = itau + n
355 CALL cgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
356 $ lwork-iwrk+1, ierr )
357
358 IF( wantvl ) THEN
359
360
361
362
363 side = 'L'
364 CALL clacpy(
'L', n, n, a, lda, vl, ldvl )
365
366
367
368
369
370 CALL cunghr( n, ilo, ihi, vl, ldvl, work( itau ),
371 $ work( iwrk ),
372 $ lwork-iwrk+1, ierr )
373
374
375
376
377
378 iwrk = itau
379 CALL chseqr(
'S',
'V', n, ilo, ihi, a, lda, w, vl, ldvl,
380 $ work( iwrk ), lwork-iwrk+1, info )
381
382 IF( wantvr ) THEN
383
384
385
386
387 side = 'B'
388 CALL clacpy(
'F', n, n, vl, ldvl, vr, ldvr )
389 END IF
390
391 ELSE IF( wantvr ) THEN
392
393
394
395
396 side = 'R'
397 CALL clacpy(
'L', n, n, a, lda, vr, ldvr )
398
399
400
401
402
403 CALL cunghr( n, ilo, ihi, vr, ldvr, work( itau ),
404 $ work( iwrk ),
405 $ lwork-iwrk+1, ierr )
406
407
408
409
410
411 iwrk = itau
412 CALL chseqr(
'S',
'V', n, ilo, ihi, a, lda, w, vr, ldvr,
413 $ work( iwrk ), lwork-iwrk+1, info )
414
415 ELSE
416
417
418
419
420
421 iwrk = itau
422 CALL chseqr(
'E',
'N', n, ilo, ihi, a, lda, w, vr, ldvr,
423 $ work( iwrk ), lwork-iwrk+1, info )
424 END IF
425
426
427
428 IF( info.NE.0 )
429 $ GO TO 50
430
431 IF( wantvl .OR. wantvr ) THEN
432
433
434
435
436
437 irwork = ibal + n
438 CALL ctrevc3( side,
'B',
SELECT, n, a, lda, vl, ldvl, vr,
439 $ ldvr,
440 $ n, nout, work( iwrk ), lwork-iwrk+1,
441 $ rwork( irwork ), n, ierr )
442 END IF
443
444 IF( wantvl ) THEN
445
446
447
448
449
450 CALL cgebak(
'B',
'L', n, ilo, ihi, rwork( ibal ), n, vl,
451 $ ldvl,
452 $ ierr )
453
454
455
456 DO 20 i = 1, n
457 scl = one /
scnrm2( n, vl( 1, i ), 1 )
458 CALL csscal( n, scl, vl( 1, i ), 1 )
459 DO 10 k = 1, n
460 rwork( irwork+k-1 ) = real( vl( k, i ) )**2 +
461 $ aimag( vl( k, i ) )**2
462 10 CONTINUE
463 k =
isamax( n, rwork( irwork ), 1 )
464 tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
465 CALL cscal( n, tmp, vl( 1, i ), 1 )
466 vl( k, i ) = cmplx( real( vl( k, i ) ), zero )
467 20 CONTINUE
468 END IF
469
470 IF( wantvr ) THEN
471
472
473
474
475
476 CALL cgebak(
'B',
'R', n, ilo, ihi, rwork( ibal ), n, vr,
477 $ ldvr,
478 $ ierr )
479
480
481
482 DO 40 i = 1, n
483 scl = one /
scnrm2( n, vr( 1, i ), 1 )
484 CALL csscal( n, scl, vr( 1, i ), 1 )
485 DO 30 k = 1, n
486 rwork( irwork+k-1 ) = real( vr( k, i ) )**2 +
487 $ aimag( vr( k, i ) )**2
488 30 CONTINUE
489 k =
isamax( n, rwork( irwork ), 1 )
490 tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
491 CALL cscal( n, tmp, vr( 1, i ), 1 )
492 vr( k, i ) = cmplx( real( vr( k, i ) ), zero )
493 40 CONTINUE
494 END IF
495
496
497
498 50 CONTINUE
499 IF( scalea ) THEN
500 CALL clascl(
'G', 0, 0, cscale, anrm, n-info, 1,
501 $ w( info+1 ),
502 $ max( n-info, 1 ), ierr )
503 IF( info.GT.0 ) THEN
504 CALL clascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, w, n,
505 $ ierr )
506 END IF
507 END IF
508
510 RETURN
511
512
513
subroutine xerbla(srname, info)
subroutine cgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
CGEBAK
subroutine cgebal(job, n, a, lda, ilo, ihi, scale, info)
CGEBAL
subroutine cgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
CGEHRD
subroutine chseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
CHSEQR
integer function isamax(n, sx, incx)
ISAMAX
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
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 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.
logical function lsame(ca, cb)
LSAME
real(wp) function scnrm2(n, x, incx)
SCNRM2
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine csscal(n, sa, cx, incx)
CSSCAL
subroutine cscal(n, ca, cx, incx)
CSCAL
subroutine ctrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, rwork, lrwork, info)
CTREVC3
subroutine cunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
CUNGHR