192 implicit none
193
194
195
196
197
198
199 CHARACTER JOBVL, JOBVR
200 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
201
202
203 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
204 $ WI( * ), WORK( * ), WR( * )
205
206
207
208
209
210 DOUBLE PRECISION ZERO, ONE
211 parameter( zero = 0.0d0, one = 1.0d0 )
212
213
214 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
215 CHARACTER SIDE
216 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
217 $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
218 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
219 $ SN
220
221
222 LOGICAL SELECT( 1 )
223 DOUBLE PRECISION DUM( 1 )
224
225
229
230
231 LOGICAL LSAME
232 INTEGER IDAMAX, ILAENV
233 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
236
237
238 INTRINSIC max, sqrt
239
240
241
242
243
244 info = 0
245 lquery = ( lwork.EQ.-1 )
246 wantvl =
lsame( jobvl,
'V' )
247 wantvr =
lsame( jobvr,
'V' )
248 IF( ( .NOT.wantvl ) .AND. ( .NOT.
lsame( jobvl,
'N' ) ) )
THEN
249 info = -1
250 ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.
lsame( jobvr,
'N' ) ) )
THEN
251 info = -2
252 ELSE IF( n.LT.0 ) THEN
253 info = -3
254 ELSE IF( lda.LT.max( 1, n ) ) THEN
255 info = -5
256 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
257 info = -9
258 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
259 info = -11
260 END IF
261
262
263
264
265
266
267
268
269
270
271
272 IF( info.EQ.0 ) THEN
273 IF( n.EQ.0 ) THEN
274 minwrk = 1
275 maxwrk = 1
276 ELSE
277 maxwrk = 2*n + n*
ilaenv( 1,
'DGEHRD',
' ', n, 1, n, 0 )
278 IF( wantvl ) THEN
279 minwrk = 4*n
280 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
281 $ 'DORGHR', ' ', n, 1, n, -1 ) )
282 CALL dhseqr(
'S',
'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
283 $ work, -1, info )
284 hswork = int( work(1) )
285 maxwrk = max( maxwrk, n + 1, n + hswork )
286 CALL dtrevc3(
'L',
'B',
SELECT, n, a, lda,
287 $ vl, ldvl, vr, ldvr, n, nout,
288 $ work, -1, ierr )
289 lwork_trevc = int( work(1) )
290 maxwrk = max( maxwrk, n + lwork_trevc )
291 maxwrk = max( maxwrk, 4*n )
292 ELSE IF( wantvr ) THEN
293 minwrk = 4*n
294 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
295 $ 'DORGHR', ' ', n, 1, n, -1 ) )
296 CALL dhseqr(
'S',
'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
297 $ work, -1, info )
298 hswork = int( work(1) )
299 maxwrk = max( maxwrk, n + 1, n + hswork )
300 CALL dtrevc3(
'R',
'B',
SELECT, n, a, lda,
301 $ vl, ldvl, vr, ldvr, n, nout,
302 $ work, -1, ierr )
303 lwork_trevc = int( work(1) )
304 maxwrk = max( maxwrk, n + lwork_trevc )
305 maxwrk = max( maxwrk, 4*n )
306 ELSE
307 minwrk = 3*n
308 CALL dhseqr(
'E',
'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
309 $ work, -1, info )
310 hswork = int( work(1) )
311 maxwrk = max( maxwrk, n + 1, n + hswork )
312 END IF
313 maxwrk = max( maxwrk, minwrk )
314 END IF
315 work( 1 ) = maxwrk
316
317 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
318 info = -13
319 END IF
320 END IF
321
322 IF( info.NE.0 ) THEN
323 CALL xerbla(
'DGEEV ', -info )
324 RETURN
325 ELSE IF( lquery ) THEN
326 RETURN
327 END IF
328
329
330
331 IF( n.EQ.0 )
332 $ RETURN
333
334
335
338 bignum = one / smlnum
339 CALL dlabad( smlnum, bignum )
340 smlnum = sqrt( smlnum ) / eps
341 bignum = one / smlnum
342
343
344
345 anrm =
dlange(
'M', n, n, a, lda, dum )
346 scalea = .false.
347 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
348 scalea = .true.
349 cscale = smlnum
350 ELSE IF( anrm.GT.bignum ) THEN
351 scalea = .true.
352 cscale = bignum
353 END IF
354 IF( scalea )
355 $
CALL dlascl(
'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
356
357
358
359
360 ibal = 1
361 CALL dgebal(
'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
362
363
364
365
366 itau = ibal + n
367 iwrk = itau + n
368 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
369 $ lwork-iwrk+1, ierr )
370
371 IF( wantvl ) THEN
372
373
374
375
376 side = 'L'
377 CALL dlacpy(
'L', n, n, a, lda, vl, ldvl )
378
379
380
381
382 CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
383 $ lwork-iwrk+1, ierr )
384
385
386
387
388 iwrk = itau
389 CALL dhseqr(
'S',
'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
390 $ work( iwrk ), lwork-iwrk+1, info )
391
392 IF( wantvr ) THEN
393
394
395
396
397 side = 'B'
398 CALL dlacpy(
'F', n, n, vl, ldvl, vr, ldvr )
399 END IF
400
401 ELSE IF( wantvr ) THEN
402
403
404
405
406 side = 'R'
407 CALL dlacpy(
'L', n, n, a, lda, vr, ldvr )
408
409
410
411
412 CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
413 $ lwork-iwrk+1, ierr )
414
415
416
417
418 iwrk = itau
419 CALL dhseqr(
'S',
'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
420 $ work( iwrk ), lwork-iwrk+1, info )
421
422 ELSE
423
424
425
426
427 iwrk = itau
428 CALL dhseqr(
'E',
'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
429 $ work( iwrk ), lwork-iwrk+1, info )
430 END IF
431
432
433
434 IF( info.NE.0 )
435 $ GO TO 50
436
437 IF( wantvl .OR. wantvr ) THEN
438
439
440
441
442 CALL dtrevc3( side,
'B',
SELECT, n, a, lda, vl, ldvl, vr, ldvr,
443 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
444 END IF
445
446 IF( wantvl ) THEN
447
448
449
450
451 CALL dgebak(
'B',
'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
452 $ ierr )
453
454
455
456 DO 20 i = 1, n
457 IF( wi( i ).EQ.zero ) THEN
458 scl = one /
dnrm2( n, vl( 1, i ), 1 )
459 CALL dscal( n, scl, vl( 1, i ), 1 )
460 ELSE IF( wi( i ).GT.zero ) THEN
462 $
dnrm2( n, vl( 1, i+1 ), 1 ) )
463 CALL dscal( n, scl, vl( 1, i ), 1 )
464 CALL dscal( n, scl, vl( 1, i+1 ), 1 )
465 DO 10 k = 1, n
466 work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
467 10 CONTINUE
468 k =
idamax( n, work( iwrk ), 1 )
469 CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
470 CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
471 vl( k, i+1 ) = zero
472 END IF
473 20 CONTINUE
474 END IF
475
476 IF( wantvr ) THEN
477
478
479
480
481 CALL dgebak(
'B',
'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
482 $ ierr )
483
484
485
486 DO 40 i = 1, n
487 IF( wi( i ).EQ.zero ) THEN
488 scl = one /
dnrm2( n, vr( 1, i ), 1 )
489 CALL dscal( n, scl, vr( 1, i ), 1 )
490 ELSE IF( wi( i ).GT.zero ) THEN
492 $
dnrm2( n, vr( 1, i+1 ), 1 ) )
493 CALL dscal( n, scl, vr( 1, i ), 1 )
494 CALL dscal( n, scl, vr( 1, i+1 ), 1 )
495 DO 30 k = 1, n
496 work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
497 30 CONTINUE
498 k =
idamax( n, work( iwrk ), 1 )
499 CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
500 CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
501 vr( k, i+1 ) = zero
502 END IF
503 40 CONTINUE
504 END IF
505
506
507
508 50 CONTINUE
509 IF( scalea ) THEN
510 CALL dlascl(
'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
511 $ max( n-info, 1 ), ierr )
512 CALL dlascl(
'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
513 $ max( n-info, 1 ), ierr )
514 IF( info.GT.0 ) THEN
515 CALL dlascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
516 $ ierr )
517 CALL dlascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
518 $ ierr )
519 END IF
520 END IF
521
522 work( 1 ) = maxwrk
523 RETURN
524
525
526
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
integer function idamax(N, DX, INCX)
IDAMAX
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dscal(N, DA, DX, INCX)
DSCAL
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
subroutine dtrevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
DTREVC3
real(wp) function dnrm2(n, x, incx)
DNRM2