192 implicit none
193
194
195
196
197
198
199 CHARACTER JOBVL, JOBVR
200 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
201
202
203 REAL A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
204 $ WI( * ), WORK( * ), WR( * )
205
206
207
208
209
210 REAL ZERO, ONE
211 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
219 $ SN
220
221
222 LOGICAL SELECT( 1 )
223 REAL DUM( 1 )
224
225
228
229
230 LOGICAL LSAME
231 INTEGER ISAMAX, ILAENV
232 REAL SLAMCH, SLANGE, SLAPY2, SNRM2, SROUNDUP_LWORK
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,
'SGEHRD',
' ', n, 1, n, 0 )
277 IF( wantvl ) THEN
278 minwrk = 4*n
279 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
280 $ 'SORGHR', ' ', n, 1, n, -1 ) )
281 CALL shseqr(
'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 strevc3(
'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 $ 'SORGHR', ' ', n, 1, n, -1 ) )
295 CALL shseqr(
'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 strevc3(
'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 shseqr(
'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
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(
'SGEEV ', -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 =
slange(
'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 slascl(
'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
354
355
356
357
358 ibal = 1
359 CALL sgebal(
'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
360
361
362
363
364 itau = ibal + n
365 iwrk = itau + n
366 CALL sgehrd( 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 slacpy(
'L', n, n, a, lda, vl, ldvl )
376
377
378
379
380 CALL sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
381 $ lwork-iwrk+1, ierr )
382
383
384
385
386 iwrk = itau
387 CALL shseqr(
'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 slacpy(
'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 slacpy(
'L', n, n, a, lda, vr, ldvr )
406
407
408
409
410 CALL sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
411 $ lwork-iwrk+1, ierr )
412
413
414
415
416 iwrk = itau
417 CALL shseqr(
'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 shseqr(
'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 strevc3( 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 sgebak(
'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 /
snrm2( n, vl( 1, i ), 1 )
457 CALL sscal( n, scl, vl( 1, i ), 1 )
458 ELSE IF( wi( i ).GT.zero ) THEN
460 $
snrm2( n, vl( 1, i+1 ), 1 ) )
461 CALL sscal( n, scl, vl( 1, i ), 1 )
462 CALL sscal( 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 =
isamax( n, work( iwrk ), 1 )
467 CALL slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
468 CALL srot( 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 sgebak(
'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 /
snrm2( n, vr( 1, i ), 1 )
487 CALL sscal( n, scl, vr( 1, i ), 1 )
488 ELSE IF( wi( i ).GT.zero ) THEN
490 $
snrm2( n, vr( 1, i+1 ), 1 ) )
491 CALL sscal( n, scl, vr( 1, i ), 1 )
492 CALL sscal( 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 =
isamax( n, work( iwrk ), 1 )
497 CALL slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
498 CALL srot( 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 slascl(
'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
509 $ max( n-info, 1 ), ierr )
510 CALL slascl(
'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 slascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
514 $ ierr )
515 CALL slascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
516 $ ierr )
517 END IF
518 END IF
519
521 RETURN
522
523
524
subroutine xerbla(srname, info)
subroutine sgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
SGEBAK
subroutine sgebal(job, n, a, lda, ilo, ihi, scale, info)
SGEBAL
subroutine sgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
SGEHRD
subroutine shseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
SHSEQR
integer function isamax(n, sx, incx)
ISAMAX
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
logical function lsame(ca, cb)
LSAME
real(wp) function snrm2(n, x, incx)
SNRM2
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sscal(n, sa, sx, incx)
SSCAL
subroutine strevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
STREVC3
subroutine sorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
SORGHR