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
228
229
230 LOGICAL LSAME
231 INTEGER IDAMAX, ILAENV
232 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
235
236
237 INTRINSIC max, sqrt
238
239
240
241
242
243 info = 0
244 lquery = ( lwork.EQ.-1 )
245 wantvl =
lsame( jobvl,
'V' )
246 wantvr =
lsame( jobvr,
'V' )
247 IF( ( .NOT.wantvl ) .AND. ( .NOT.
lsame( jobvl,
'N' ) ) )
THEN
248 info = -1
249 ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.
lsame( jobvr,
'N' ) ) )
THEN
250 info = -2
251 ELSE IF( n.LT.0 ) THEN
252 info = -3
253 ELSE IF( lda.LT.max( 1, n ) ) THEN
254 info = -5
255 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
256 info = -9
257 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
258 info = -11
259 END IF
260
261
262
263
264
265
266
267
268
269
270
271 IF( info.EQ.0 ) THEN
272 IF( n.EQ.0 ) THEN
273 minwrk = 1
274 maxwrk = 1
275 ELSE
276 maxwrk = 2*n + n*
ilaenv( 1,
'DGEHRD',
' ', n, 1, n, 0 )
277 IF( wantvl ) THEN
278 minwrk = 4*n
279 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
280 $ 'DORGHR', ' ', n, 1, n, -1 ) )
281 CALL dhseqr(
'S',
'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
282 $ work, -1, info )
283 hswork = int( work(1) )
284 maxwrk = max( maxwrk, n + 1, n + hswork )
285 CALL dtrevc3(
'L',
'B',
SELECT, n, a, lda,
286 $ vl, ldvl, vr, ldvr, n, nout,
287 $ work, -1, ierr )
288 lwork_trevc = int( work(1) )
289 maxwrk = max( maxwrk, n + lwork_trevc )
290 maxwrk = max( maxwrk, 4*n )
291 ELSE IF( wantvr ) THEN
292 minwrk = 4*n
293 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
294 $ 'DORGHR', ' ', n, 1, n, -1 ) )
295 CALL dhseqr(
'S',
'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
296 $ work, -1, info )
297 hswork = int( work(1) )
298 maxwrk = max( maxwrk, n + 1, n + hswork )
299 CALL dtrevc3(
'R',
'B',
SELECT, n, a, lda,
300 $ vl, ldvl, vr, ldvr, n, nout,
301 $ work, -1, ierr )
302 lwork_trevc = int( work(1) )
303 maxwrk = max( maxwrk, n + lwork_trevc )
304 maxwrk = max( maxwrk, 4*n )
305 ELSE
306 minwrk = 3*n
307 CALL dhseqr(
'E',
'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
308 $ work, -1, info )
309 hswork = int( work(1) )
310 maxwrk = max( maxwrk, n + 1, n + hswork )
311 END IF
312 maxwrk = max( maxwrk, minwrk )
313 END IF
314 work( 1 ) = maxwrk
315
316 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
317 info = -13
318 END IF
319 END IF
320
321 IF( info.NE.0 ) THEN
322 CALL xerbla(
'DGEEV ', -info )
323 RETURN
324 ELSE IF( lquery ) THEN
325 RETURN
326 END IF
327
328
329
330 IF( n.EQ.0 )
331 $ RETURN
332
333
334
337 bignum = one / smlnum
338 smlnum = sqrt( smlnum ) / eps
339 bignum = one / smlnum
340
341
342
343 anrm =
dlange(
'M', n, n, a, lda, dum )
344 scalea = .false.
345 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
346 scalea = .true.
347 cscale = smlnum
348 ELSE IF( anrm.GT.bignum ) THEN
349 scalea = .true.
350 cscale = bignum
351 END IF
352 IF( scalea )
353 $
CALL dlascl(
'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
354
355
356
357
358 ibal = 1
359 CALL dgebal(
'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
360
361
362
363
364 itau = ibal + n
365 iwrk = itau + n
366 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
367 $ lwork-iwrk+1, ierr )
368
369 IF( wantvl ) THEN
370
371
372
373
374 side = 'L'
375 CALL dlacpy(
'L', n, n, a, lda, vl, ldvl )
376
377
378
379
380 CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
381 $ lwork-iwrk+1, ierr )
382
383
384
385
386 iwrk = itau
387 CALL dhseqr(
'S',
'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
388 $ work( iwrk ), lwork-iwrk+1, info )
389
390 IF( wantvr ) THEN
391
392
393
394
395 side = 'B'
396 CALL dlacpy(
'F', n, n, vl, ldvl, vr, ldvr )
397 END IF
398
399 ELSE IF( wantvr ) THEN
400
401
402
403
404 side = 'R'
405 CALL dlacpy(
'L', n, n, a, lda, vr, ldvr )
406
407
408
409
410 CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
411 $ lwork-iwrk+1, ierr )
412
413
414
415
416 iwrk = itau
417 CALL dhseqr(
'S',
'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
418 $ work( iwrk ), lwork-iwrk+1, info )
419
420 ELSE
421
422
423
424
425 iwrk = itau
426 CALL dhseqr(
'E',
'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
427 $ work( iwrk ), lwork-iwrk+1, info )
428 END IF
429
430
431
432 IF( info.NE.0 )
433 $ GO TO 50
434
435 IF( wantvl .OR. wantvr ) THEN
436
437
438
439
440 CALL dtrevc3( side,
'B',
SELECT, n, a, lda, vl, ldvl, vr, ldvr,
441 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
442 END IF
443
444 IF( wantvl ) THEN
445
446
447
448
449 CALL dgebak(
'B',
'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
450 $ ierr )
451
452
453
454 DO 20 i = 1, n
455 IF( wi( i ).EQ.zero ) THEN
456 scl = one /
dnrm2( n, vl( 1, i ), 1 )
457 CALL dscal( n, scl, vl( 1, i ), 1 )
458 ELSE IF( wi( i ).GT.zero ) THEN
460 $
dnrm2( n, vl( 1, i+1 ), 1 ) )
461 CALL dscal( n, scl, vl( 1, i ), 1 )
462 CALL dscal( n, scl, vl( 1, i+1 ), 1 )
463 DO 10 k = 1, n
464 work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
465 10 CONTINUE
466 k =
idamax( n, work( iwrk ), 1 )
467 CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
468 CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
469 vl( k, i+1 ) = zero
470 END IF
471 20 CONTINUE
472 END IF
473
474 IF( wantvr ) THEN
475
476
477
478
479 CALL dgebak(
'B',
'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
480 $ ierr )
481
482
483
484 DO 40 i = 1, n
485 IF( wi( i ).EQ.zero ) THEN
486 scl = one /
dnrm2( n, vr( 1, i ), 1 )
487 CALL dscal( n, scl, vr( 1, i ), 1 )
488 ELSE IF( wi( i ).GT.zero ) THEN
490 $
dnrm2( n, vr( 1, i+1 ), 1 ) )
491 CALL dscal( n, scl, vr( 1, i ), 1 )
492 CALL dscal( n, scl, vr( 1, i+1 ), 1 )
493 DO 30 k = 1, n
494 work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
495 30 CONTINUE
496 k =
idamax( n, work( iwrk ), 1 )
497 CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
498 CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
499 vr( k, i+1 ) = zero
500 END IF
501 40 CONTINUE
502 END IF
503
504
505
506 50 CONTINUE
507 IF( scalea ) THEN
508 CALL dlascl(
'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
509 $ max( n-info, 1 ), ierr )
510 CALL dlascl(
'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
511 $ max( n-info, 1 ), ierr )
512 IF( info.GT.0 ) THEN
513 CALL dlascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
514 $ ierr )
515 CALL dlascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
516 $ ierr )
517 END IF
518 END IF
519
520 work( 1 ) = maxwrk
521 RETURN
522
523
524
subroutine xerbla(srname, info)
subroutine dgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
DGEBAK
subroutine dgebal(job, n, a, lda, ilo, ihi, scale, info)
DGEBAL
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
integer function idamax(n, dx, incx)
IDAMAX
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
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 ...
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
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.
logical function lsame(ca, cb)
LSAME
real(wp) function dnrm2(n, x, incx)
DNRM2
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine dtrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
DTREVC3
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR