LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine dtrevc3 ( character  SIDE,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DTREVC3

Download DTREVC3 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 DTREVC3 computes some or all of the right and/or left eigenvectors of
 a real upper quasi-triangular matrix T.
 Matrices of this type are produced by the Schur factorization of
 a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.

 The right eigenvector x and the left eigenvector y of T corresponding
 to an eigenvalue w are defined by:

    T*x = w*x,     (y**H)*T = w*(y**H)

 where y**H denotes the conjugate transpose of y.
 The eigenvalues are not input to this routine, but are read directly
 from the diagonal blocks of T.

 This routine returns the matrices X and/or Y of right and left
 eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
 input matrix. If Q is the orthogonal factor that reduces a matrix
 A to Schur form T, then Q*X and Q*Y are the matrices of right and
 left eigenvectors of A.

 This uses a Level 3 BLAS version of the back transformation.
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'R':  compute right eigenvectors only;
          = 'L':  compute left eigenvectors only;
          = 'B':  compute both right and left eigenvectors.
[in]HOWMNY
          HOWMNY is CHARACTER*1
          = 'A':  compute all right and/or left eigenvectors;
          = 'B':  compute all right and/or left eigenvectors,
                  backtransformed by the matrices in VR and/or VL;
          = 'S':  compute selected right and/or left eigenvectors,
                  as indicated by the logical array SELECT.
[in,out]SELECT
          SELECT is LOGICAL array, dimension (N)
          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
          computed.
          If w(j) is a real eigenvalue, the corresponding real
          eigenvector is computed if SELECT(j) is .TRUE..
          If w(j) and w(j+1) are the real and imaginary parts of a
          complex eigenvalue, the corresponding complex eigenvector is
          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
          .FALSE..
          Not referenced if HOWMNY = 'A' or 'B'.
[in]N
          N is INTEGER
          The order of the matrix T. N >= 0.
[in]T
          T is DOUBLE PRECISION array, dimension (LDT,N)
          The upper quasi-triangular matrix T in Schur canonical form.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T. LDT >= max(1,N).
[in,out]VL
          VL is DOUBLE PRECISION array, dimension (LDVL,MM)
          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
          contain an N-by-N matrix Q (usually the orthogonal matrix Q
          of Schur vectors returned by DHSEQR).
          On exit, if SIDE = 'L' or 'B', VL contains:
          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
          if HOWMNY = 'B', the matrix Q*Y;
          if HOWMNY = 'S', the left eigenvectors of T specified by
                           SELECT, stored consecutively in the columns
                           of VL, in the same order as their
                           eigenvalues.
          A complex eigenvector corresponding to a complex eigenvalue
          is stored in two consecutive columns, the first holding the
          real part, and the second the imaginary part.
          Not referenced if SIDE = 'R'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.
          LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
[in,out]VR
          VR is DOUBLE PRECISION array, dimension (LDVR,MM)
          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
          contain an N-by-N matrix Q (usually the orthogonal matrix Q
          of Schur vectors returned by DHSEQR).
          On exit, if SIDE = 'R' or 'B', VR contains:
          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
          if HOWMNY = 'B', the matrix Q*X;
          if HOWMNY = 'S', the right eigenvectors of T specified by
                           SELECT, stored consecutively in the columns
                           of VR, in the same order as their
                           eigenvalues.
          A complex eigenvector corresponding to a complex eigenvalue
          is stored in two consecutive columns, the first holding the
          real part and the second the imaginary part.
          Not referenced if SIDE = 'L'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.
          LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
[in]MM
          MM is INTEGER
          The number of columns in the arrays VL and/or VR. MM >= M.
[out]M
          M is INTEGER
          The number of columns in the arrays VL and/or VR actually
          used to store the eigenvectors.
          If HOWMNY = 'A' or 'B', M is set to N.
          Each selected real eigenvector occupies one column and each
          selected complex eigenvector occupies two columns.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
[in]LWORK
          LWORK is INTEGER
          The dimension of array WORK. LWORK >= max(1,3*N).
          For optimum performance, LWORK >= N + 2*N*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Further Details: 
  The algorithm used in this program is basically backward (forward)
  substitution, with scaling to make the the code robust against
  possible overflow.

  Each eigenvector is normalized so that the element of largest
  magnitude has magnitude 1; here the magnitude of a complex number
  (x,y) is taken to be |x| + |y|.

Definition at line 241 of file dtrevc3.f.

241  IMPLICIT NONE
242 *
243 * -- LAPACK computational routine (version 3.7.0) --
244 * -- LAPACK is a software package provided by Univ. of Tennessee, --
245 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
246 * December 2016
247 *
248 * .. Scalar Arguments ..
249  CHARACTER howmny, side
250  INTEGER info, ldt, ldvl, ldvr, lwork, m, mm, n
251 * ..
252 * .. Array Arguments ..
253  LOGICAL select( * )
254  DOUBLE PRECISION t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
255  $ work( * )
256 * ..
257 *
258 * =====================================================================
259 *
260 * .. Parameters ..
261  DOUBLE PRECISION zero, one
262  parameter ( zero = 0.0d+0, one = 1.0d+0 )
263  INTEGER nbmin, nbmax
264  parameter ( nbmin = 8, nbmax = 128 )
265 * ..
266 * .. Local Scalars ..
267  LOGICAL allv, bothv, leftv, lquery, over, pair,
268  $ rightv, somev
269  INTEGER i, ierr, ii, ip, is, j, j1, j2, jnxt, k, ki,
270  $ iv, maxwrk, nb, ki2
271  DOUBLE PRECISION beta, bignum, emax, ovfl, rec, remax, scale,
272  $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
273  $ xnorm
274 * ..
275 * .. External Functions ..
276  LOGICAL lsame
277  INTEGER idamax, ilaenv
278  DOUBLE PRECISION ddot, dlamch
279  EXTERNAL lsame, idamax, ilaenv, ddot, dlamch
280 * ..
281 * .. External Subroutines ..
282  EXTERNAL daxpy, dcopy, dgemv, dlaln2, dscal, xerbla,
283  $ dgemm, dlaset, dlabad
284 * ..
285 * .. Intrinsic Functions ..
286  INTRINSIC abs, max, sqrt
287 * ..
288 * .. Local Arrays ..
289  DOUBLE PRECISION x( 2, 2 )
290  INTEGER iscomplex( nbmax )
291 * ..
292 * .. Executable Statements ..
293 *
294 * Decode and test the input parameters
295 *
296  bothv = lsame( side, 'B' )
297  rightv = lsame( side, 'R' ) .OR. bothv
298  leftv = lsame( side, 'L' ) .OR. bothv
299 *
300  allv = lsame( howmny, 'A' )
301  over = lsame( howmny, 'B' )
302  somev = lsame( howmny, 'S' )
303 *
304  info = 0
305  nb = ilaenv( 1, 'DTREVC', side // howmny, n, -1, -1, -1 )
306  maxwrk = n + 2*n*nb
307  work(1) = maxwrk
308  lquery = ( lwork.EQ.-1 )
309  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
310  info = -1
311  ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
312  info = -2
313  ELSE IF( n.LT.0 ) THEN
314  info = -4
315  ELSE IF( ldt.LT.max( 1, n ) ) THEN
316  info = -6
317  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
318  info = -8
319  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
320  info = -10
321  ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery ) THEN
322  info = -14
323  ELSE
324 *
325 * Set M to the number of columns required to store the selected
326 * eigenvectors, standardize the array SELECT if necessary, and
327 * test MM.
328 *
329  IF( somev ) THEN
330  m = 0
331  pair = .false.
332  DO 10 j = 1, n
333  IF( pair ) THEN
334  pair = .false.
335  SELECT( j ) = .false.
336  ELSE
337  IF( j.LT.n ) THEN
338  IF( t( j+1, j ).EQ.zero ) THEN
339  IF( SELECT( j ) )
340  $ m = m + 1
341  ELSE
342  pair = .true.
343  IF( SELECT( j ) .OR. SELECT( j+1 ) ) THEN
344  SELECT( j ) = .true.
345  m = m + 2
346  END IF
347  END IF
348  ELSE
349  IF( SELECT( n ) )
350  $ m = m + 1
351  END IF
352  END IF
353  10 CONTINUE
354  ELSE
355  m = n
356  END IF
357 *
358  IF( mm.LT.m ) THEN
359  info = -11
360  END IF
361  END IF
362  IF( info.NE.0 ) THEN
363  CALL xerbla( 'DTREVC3', -info )
364  RETURN
365  ELSE IF( lquery ) THEN
366  RETURN
367  END IF
368 *
369 * Quick return if possible.
370 *
371  IF( n.EQ.0 )
372  $ RETURN
373 *
374 * Use blocked version of back-transformation if sufficient workspace.
375 * Zero-out the workspace to avoid potential NaN propagation.
376 *
377  IF( over .AND. lwork .GE. n + 2*n*nbmin ) THEN
378  nb = (lwork - n) / (2*n)
379  nb = min( nb, nbmax )
380  CALL dlaset( 'F', n, 1+2*nb, zero, zero, work, n )
381  ELSE
382  nb = 1
383  END IF
384 *
385 * Set the constants to control overflow.
386 *
387  unfl = dlamch( 'Safe minimum' )
388  ovfl = one / unfl
389  CALL dlabad( unfl, ovfl )
390  ulp = dlamch( 'Precision' )
391  smlnum = unfl*( n / ulp )
392  bignum = ( one-ulp ) / smlnum
393 *
394 * Compute 1-norm of each column of strictly upper triangular
395 * part of T to control overflow in triangular solver.
396 *
397  work( 1 ) = zero
398  DO 30 j = 2, n
399  work( j ) = zero
400  DO 20 i = 1, j - 1
401  work( j ) = work( j ) + abs( t( i, j ) )
402  20 CONTINUE
403  30 CONTINUE
404 *
405 * Index IP is used to specify the real or complex eigenvalue:
406 * IP = 0, real eigenvalue,
407 * 1, first of conjugate complex pair: (wr,wi)
408 * -1, second of conjugate complex pair: (wr,wi)
409 * ISCOMPLEX array stores IP for each column in current block.
410 *
411  IF( rightv ) THEN
412 *
413 * ============================================================
414 * Compute right eigenvectors.
415 *
416 * IV is index of column in current block.
417 * For complex right vector, uses IV-1 for real part and IV for complex part.
418 * Non-blocked version always uses IV=2;
419 * blocked version starts with IV=NB, goes down to 1 or 2.
420 * (Note the "0-th" column is used for 1-norms computed above.)
421  iv = 2
422  IF( nb.GT.2 ) THEN
423  iv = nb
424  END IF
425 
426  ip = 0
427  is = m
428  DO 140 ki = n, 1, -1
429  IF( ip.EQ.-1 ) THEN
430 * previous iteration (ki+1) was second of conjugate pair,
431 * so this ki is first of conjugate pair; skip to end of loop
432  ip = 1
433  GO TO 140
434  ELSE IF( ki.EQ.1 ) THEN
435 * last column, so this ki must be real eigenvalue
436  ip = 0
437  ELSE IF( t( ki, ki-1 ).EQ.zero ) THEN
438 * zero on sub-diagonal, so this ki is real eigenvalue
439  ip = 0
440  ELSE
441 * non-zero on sub-diagonal, so this ki is second of conjugate pair
442  ip = -1
443  END IF
444 
445  IF( somev ) THEN
446  IF( ip.EQ.0 ) THEN
447  IF( .NOT.SELECT( ki ) )
448  $ GO TO 140
449  ELSE
450  IF( .NOT.SELECT( ki-1 ) )
451  $ GO TO 140
452  END IF
453  END IF
454 *
455 * Compute the KI-th eigenvalue (WR,WI).
456 *
457  wr = t( ki, ki )
458  wi = zero
459  IF( ip.NE.0 )
460  $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
461  $ sqrt( abs( t( ki-1, ki ) ) )
462  smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
463 *
464  IF( ip.EQ.0 ) THEN
465 *
466 * --------------------------------------------------------
467 * Real right eigenvector
468 *
469  work( ki + iv*n ) = one
470 *
471 * Form right-hand side.
472 *
473  DO 50 k = 1, ki - 1
474  work( k + iv*n ) = -t( k, ki )
475  50 CONTINUE
476 *
477 * Solve upper quasi-triangular system:
478 * [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
479 *
480  jnxt = ki - 1
481  DO 60 j = ki - 1, 1, -1
482  IF( j.GT.jnxt )
483  $ GO TO 60
484  j1 = j
485  j2 = j
486  jnxt = j - 1
487  IF( j.GT.1 ) THEN
488  IF( t( j, j-1 ).NE.zero ) THEN
489  j1 = j - 1
490  jnxt = j - 2
491  END IF
492  END IF
493 *
494  IF( j1.EQ.j2 ) THEN
495 *
496 * 1-by-1 diagonal block
497 *
498  CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
499  $ ldt, one, one, work( j+iv*n ), n, wr,
500  $ zero, x, 2, scale, xnorm, ierr )
501 *
502 * Scale X(1,1) to avoid overflow when updating
503 * the right-hand side.
504 *
505  IF( xnorm.GT.one ) THEN
506  IF( work( j ).GT.bignum / xnorm ) THEN
507  x( 1, 1 ) = x( 1, 1 ) / xnorm
508  scale = scale / xnorm
509  END IF
510  END IF
511 *
512 * Scale if necessary
513 *
514  IF( scale.NE.one )
515  $ CALL dscal( ki, scale, work( 1+iv*n ), 1 )
516  work( j+iv*n ) = x( 1, 1 )
517 *
518 * Update right-hand side
519 *
520  CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
521  $ work( 1+iv*n ), 1 )
522 *
523  ELSE
524 *
525 * 2-by-2 diagonal block
526 *
527  CALL dlaln2( .false., 2, 1, smin, one,
528  $ t( j-1, j-1 ), ldt, one, one,
529  $ work( j-1+iv*n ), n, wr, zero, x, 2,
530  $ scale, xnorm, ierr )
531 *
532 * Scale X(1,1) and X(2,1) to avoid overflow when
533 * updating the right-hand side.
534 *
535  IF( xnorm.GT.one ) THEN
536  beta = max( work( j-1 ), work( j ) )
537  IF( beta.GT.bignum / xnorm ) THEN
538  x( 1, 1 ) = x( 1, 1 ) / xnorm
539  x( 2, 1 ) = x( 2, 1 ) / xnorm
540  scale = scale / xnorm
541  END IF
542  END IF
543 *
544 * Scale if necessary
545 *
546  IF( scale.NE.one )
547  $ CALL dscal( ki, scale, work( 1+iv*n ), 1 )
548  work( j-1+iv*n ) = x( 1, 1 )
549  work( j +iv*n ) = x( 2, 1 )
550 *
551 * Update right-hand side
552 *
553  CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
554  $ work( 1+iv*n ), 1 )
555  CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
556  $ work( 1+iv*n ), 1 )
557  END IF
558  60 CONTINUE
559 *
560 * Copy the vector x or Q*x to VR and normalize.
561 *
562  IF( .NOT.over ) THEN
563 * ------------------------------
564 * no back-transform: copy x to VR and normalize.
565  CALL dcopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
566 *
567  ii = idamax( ki, vr( 1, is ), 1 )
568  remax = one / abs( vr( ii, is ) )
569  CALL dscal( ki, remax, vr( 1, is ), 1 )
570 *
571  DO 70 k = ki + 1, n
572  vr( k, is ) = zero
573  70 CONTINUE
574 *
575  ELSE IF( nb.EQ.1 ) THEN
576 * ------------------------------
577 * version 1: back-transform each vector with GEMV, Q*x.
578  IF( ki.GT.1 )
579  $ CALL dgemv( 'N', n, ki-1, one, vr, ldvr,
580  $ work( 1 + iv*n ), 1, work( ki + iv*n ),
581  $ vr( 1, ki ), 1 )
582 *
583  ii = idamax( n, vr( 1, ki ), 1 )
584  remax = one / abs( vr( ii, ki ) )
585  CALL dscal( n, remax, vr( 1, ki ), 1 )
586 *
587  ELSE
588 * ------------------------------
589 * version 2: back-transform block of vectors with GEMM
590 * zero out below vector
591  DO k = ki + 1, n
592  work( k + iv*n ) = zero
593  END DO
594  iscomplex( iv ) = ip
595 * back-transform and normalization is done below
596  END IF
597  ELSE
598 *
599 * --------------------------------------------------------
600 * Complex right eigenvector.
601 *
602 * Initial solve
603 * [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
604 * [ ( T(KI, KI-1) T(KI, KI) ) ]
605 *
606  IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) ) THEN
607  work( ki-1 + (iv-1)*n ) = one
608  work( ki + (iv )*n ) = wi / t( ki-1, ki )
609  ELSE
610  work( ki-1 + (iv-1)*n ) = -wi / t( ki, ki-1 )
611  work( ki + (iv )*n ) = one
612  END IF
613  work( ki + (iv-1)*n ) = zero
614  work( ki-1 + (iv )*n ) = zero
615 *
616 * Form right-hand side.
617 *
618  DO 80 k = 1, ki - 2
619  work( k+(iv-1)*n ) = -work( ki-1+(iv-1)*n )*t(k,ki-1)
620  work( k+(iv )*n ) = -work( ki +(iv )*n )*t(k,ki )
621  80 CONTINUE
622 *
623 * Solve upper quasi-triangular system:
624 * [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
625 *
626  jnxt = ki - 2
627  DO 90 j = ki - 2, 1, -1
628  IF( j.GT.jnxt )
629  $ GO TO 90
630  j1 = j
631  j2 = j
632  jnxt = j - 1
633  IF( j.GT.1 ) THEN
634  IF( t( j, j-1 ).NE.zero ) THEN
635  j1 = j - 1
636  jnxt = j - 2
637  END IF
638  END IF
639 *
640  IF( j1.EQ.j2 ) THEN
641 *
642 * 1-by-1 diagonal block
643 *
644  CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
645  $ ldt, one, one, work( j+(iv-1)*n ), n,
646  $ wr, wi, x, 2, scale, xnorm, ierr )
647 *
648 * Scale X(1,1) and X(1,2) to avoid overflow when
649 * updating the right-hand side.
650 *
651  IF( xnorm.GT.one ) THEN
652  IF( work( j ).GT.bignum / xnorm ) THEN
653  x( 1, 1 ) = x( 1, 1 ) / xnorm
654  x( 1, 2 ) = x( 1, 2 ) / xnorm
655  scale = scale / xnorm
656  END IF
657  END IF
658 *
659 * Scale if necessary
660 *
661  IF( scale.NE.one ) THEN
662  CALL dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
663  CALL dscal( ki, scale, work( 1+(iv )*n ), 1 )
664  END IF
665  work( j+(iv-1)*n ) = x( 1, 1 )
666  work( j+(iv )*n ) = x( 1, 2 )
667 *
668 * Update the right-hand side
669 *
670  CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
671  $ work( 1+(iv-1)*n ), 1 )
672  CALL daxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
673  $ work( 1+(iv )*n ), 1 )
674 *
675  ELSE
676 *
677 * 2-by-2 diagonal block
678 *
679  CALL dlaln2( .false., 2, 2, smin, one,
680  $ t( j-1, j-1 ), ldt, one, one,
681  $ work( j-1+(iv-1)*n ), n, wr, wi, x, 2,
682  $ scale, xnorm, ierr )
683 *
684 * Scale X to avoid overflow when updating
685 * the right-hand side.
686 *
687  IF( xnorm.GT.one ) THEN
688  beta = max( work( j-1 ), work( j ) )
689  IF( beta.GT.bignum / xnorm ) THEN
690  rec = one / xnorm
691  x( 1, 1 ) = x( 1, 1 )*rec
692  x( 1, 2 ) = x( 1, 2 )*rec
693  x( 2, 1 ) = x( 2, 1 )*rec
694  x( 2, 2 ) = x( 2, 2 )*rec
695  scale = scale*rec
696  END IF
697  END IF
698 *
699 * Scale if necessary
700 *
701  IF( scale.NE.one ) THEN
702  CALL dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
703  CALL dscal( ki, scale, work( 1+(iv )*n ), 1 )
704  END IF
705  work( j-1+(iv-1)*n ) = x( 1, 1 )
706  work( j +(iv-1)*n ) = x( 2, 1 )
707  work( j-1+(iv )*n ) = x( 1, 2 )
708  work( j +(iv )*n ) = x( 2, 2 )
709 *
710 * Update the right-hand side
711 *
712  CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
713  $ work( 1+(iv-1)*n ), 1 )
714  CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
715  $ work( 1+(iv-1)*n ), 1 )
716  CALL daxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
717  $ work( 1+(iv )*n ), 1 )
718  CALL daxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
719  $ work( 1+(iv )*n ), 1 )
720  END IF
721  90 CONTINUE
722 *
723 * Copy the vector x or Q*x to VR and normalize.
724 *
725  IF( .NOT.over ) THEN
726 * ------------------------------
727 * no back-transform: copy x to VR and normalize.
728  CALL dcopy( ki, work( 1+(iv-1)*n ), 1, vr(1,is-1), 1 )
729  CALL dcopy( ki, work( 1+(iv )*n ), 1, vr(1,is ), 1 )
730 *
731  emax = zero
732  DO 100 k = 1, ki
733  emax = max( emax, abs( vr( k, is-1 ) )+
734  $ abs( vr( k, is ) ) )
735  100 CONTINUE
736  remax = one / emax
737  CALL dscal( ki, remax, vr( 1, is-1 ), 1 )
738  CALL dscal( ki, remax, vr( 1, is ), 1 )
739 *
740  DO 110 k = ki + 1, n
741  vr( k, is-1 ) = zero
742  vr( k, is ) = zero
743  110 CONTINUE
744 *
745  ELSE IF( nb.EQ.1 ) THEN
746 * ------------------------------
747 * version 1: back-transform each vector with GEMV, Q*x.
748  IF( ki.GT.2 ) THEN
749  CALL dgemv( 'N', n, ki-2, one, vr, ldvr,
750  $ work( 1 + (iv-1)*n ), 1,
751  $ work( ki-1 + (iv-1)*n ), vr(1,ki-1), 1)
752  CALL dgemv( 'N', n, ki-2, one, vr, ldvr,
753  $ work( 1 + (iv)*n ), 1,
754  $ work( ki + (iv)*n ), vr( 1, ki ), 1 )
755  ELSE
756  CALL dscal( n, work(ki-1+(iv-1)*n), vr(1,ki-1), 1)
757  CALL dscal( n, work(ki +(iv )*n), vr(1,ki ), 1)
758  END IF
759 *
760  emax = zero
761  DO 120 k = 1, n
762  emax = max( emax, abs( vr( k, ki-1 ) )+
763  $ abs( vr( k, ki ) ) )
764  120 CONTINUE
765  remax = one / emax
766  CALL dscal( n, remax, vr( 1, ki-1 ), 1 )
767  CALL dscal( n, remax, vr( 1, ki ), 1 )
768 *
769  ELSE
770 * ------------------------------
771 * version 2: back-transform block of vectors with GEMM
772 * zero out below vector
773  DO k = ki + 1, n
774  work( k + (iv-1)*n ) = zero
775  work( k + (iv )*n ) = zero
776  END DO
777  iscomplex( iv-1 ) = -ip
778  iscomplex( iv ) = ip
779  iv = iv - 1
780 * back-transform and normalization is done below
781  END IF
782  END IF
783 
784  IF( nb.GT.1 ) THEN
785 * --------------------------------------------------------
786 * Blocked version of back-transform
787 * For complex case, KI2 includes both vectors (KI-1 and KI)
788  IF( ip.EQ.0 ) THEN
789  ki2 = ki
790  ELSE
791  ki2 = ki - 1
792  END IF
793 
794 * Columns IV:NB of work are valid vectors.
795 * When the number of vectors stored reaches NB-1 or NB,
796 * or if this was last vector, do the GEMM
797  IF( (iv.LE.2) .OR. (ki2.EQ.1) ) THEN
798  CALL dgemm( 'N', 'N', n, nb-iv+1, ki2+nb-iv, one,
799  $ vr, ldvr,
800  $ work( 1 + (iv)*n ), n,
801  $ zero,
802  $ work( 1 + (nb+iv)*n ), n )
803 * normalize vectors
804  DO k = iv, nb
805  IF( iscomplex(k).EQ.0 ) THEN
806 * real eigenvector
807  ii = idamax( n, work( 1 + (nb+k)*n ), 1 )
808  remax = one / abs( work( ii + (nb+k)*n ) )
809  ELSE IF( iscomplex(k).EQ.1 ) THEN
810 * first eigenvector of conjugate pair
811  emax = zero
812  DO ii = 1, n
813  emax = max( emax,
814  $ abs( work( ii + (nb+k )*n ) )+
815  $ abs( work( ii + (nb+k+1)*n ) ) )
816  END DO
817  remax = one / emax
818 * else if ISCOMPLEX(K).EQ.-1
819 * second eigenvector of conjugate pair
820 * reuse same REMAX as previous K
821  END IF
822  CALL dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
823  END DO
824  CALL dlacpy( 'F', n, nb-iv+1,
825  $ work( 1 + (nb+iv)*n ), n,
826  $ vr( 1, ki2 ), ldvr )
827  iv = nb
828  ELSE
829  iv = iv - 1
830  END IF
831  END IF ! blocked back-transform
832 *
833  is = is - 1
834  IF( ip.NE.0 )
835  $ is = is - 1
836  140 CONTINUE
837  END IF
838 
839  IF( leftv ) THEN
840 *
841 * ============================================================
842 * Compute left eigenvectors.
843 *
844 * IV is index of column in current block.
845 * For complex left vector, uses IV for real part and IV+1 for complex part.
846 * Non-blocked version always uses IV=1;
847 * blocked version starts with IV=1, goes up to NB-1 or NB.
848 * (Note the "0-th" column is used for 1-norms computed above.)
849  iv = 1
850  ip = 0
851  is = 1
852  DO 260 ki = 1, n
853  IF( ip.EQ.1 ) THEN
854 * previous iteration (ki-1) was first of conjugate pair,
855 * so this ki is second of conjugate pair; skip to end of loop
856  ip = -1
857  GO TO 260
858  ELSE IF( ki.EQ.n ) THEN
859 * last column, so this ki must be real eigenvalue
860  ip = 0
861  ELSE IF( t( ki+1, ki ).EQ.zero ) THEN
862 * zero on sub-diagonal, so this ki is real eigenvalue
863  ip = 0
864  ELSE
865 * non-zero on sub-diagonal, so this ki is first of conjugate pair
866  ip = 1
867  END IF
868 *
869  IF( somev ) THEN
870  IF( .NOT.SELECT( ki ) )
871  $ GO TO 260
872  END IF
873 *
874 * Compute the KI-th eigenvalue (WR,WI).
875 *
876  wr = t( ki, ki )
877  wi = zero
878  IF( ip.NE.0 )
879  $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
880  $ sqrt( abs( t( ki+1, ki ) ) )
881  smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
882 *
883  IF( ip.EQ.0 ) THEN
884 *
885 * --------------------------------------------------------
886 * Real left eigenvector
887 *
888  work( ki + iv*n ) = one
889 *
890 * Form right-hand side.
891 *
892  DO 160 k = ki + 1, n
893  work( k + iv*n ) = -t( ki, k )
894  160 CONTINUE
895 *
896 * Solve transposed quasi-triangular system:
897 * [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
898 *
899  vmax = one
900  vcrit = bignum
901 *
902  jnxt = ki + 1
903  DO 170 j = ki + 1, n
904  IF( j.LT.jnxt )
905  $ GO TO 170
906  j1 = j
907  j2 = j
908  jnxt = j + 1
909  IF( j.LT.n ) THEN
910  IF( t( j+1, j ).NE.zero ) THEN
911  j2 = j + 1
912  jnxt = j + 2
913  END IF
914  END IF
915 *
916  IF( j1.EQ.j2 ) THEN
917 *
918 * 1-by-1 diagonal block
919 *
920 * Scale if necessary to avoid overflow when forming
921 * the right-hand side.
922 *
923  IF( work( j ).GT.vcrit ) THEN
924  rec = one / vmax
925  CALL dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
926  vmax = one
927  vcrit = bignum
928  END IF
929 *
930  work( j+iv*n ) = work( j+iv*n ) -
931  $ ddot( j-ki-1, t( ki+1, j ), 1,
932  $ work( ki+1+iv*n ), 1 )
933 *
934 * Solve [ T(J,J) - WR ]**T * X = WORK
935 *
936  CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
937  $ ldt, one, one, work( j+iv*n ), n, wr,
938  $ zero, x, 2, scale, xnorm, ierr )
939 *
940 * Scale if necessary
941 *
942  IF( scale.NE.one )
943  $ CALL dscal( n-ki+1, scale, work( ki+iv*n ), 1 )
944  work( j+iv*n ) = x( 1, 1 )
945  vmax = max( abs( work( j+iv*n ) ), vmax )
946  vcrit = bignum / vmax
947 *
948  ELSE
949 *
950 * 2-by-2 diagonal block
951 *
952 * Scale if necessary to avoid overflow when forming
953 * the right-hand side.
954 *
955  beta = max( work( j ), work( j+1 ) )
956  IF( beta.GT.vcrit ) THEN
957  rec = one / vmax
958  CALL dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
959  vmax = one
960  vcrit = bignum
961  END IF
962 *
963  work( j+iv*n ) = work( j+iv*n ) -
964  $ ddot( j-ki-1, t( ki+1, j ), 1,
965  $ work( ki+1+iv*n ), 1 )
966 *
967  work( j+1+iv*n ) = work( j+1+iv*n ) -
968  $ ddot( j-ki-1, t( ki+1, j+1 ), 1,
969  $ work( ki+1+iv*n ), 1 )
970 *
971 * Solve
972 * [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
973 * [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
974 *
975  CALL dlaln2( .true., 2, 1, smin, one, t( j, j ),
976  $ ldt, one, one, work( j+iv*n ), n, wr,
977  $ zero, x, 2, scale, xnorm, ierr )
978 *
979 * Scale if necessary
980 *
981  IF( scale.NE.one )
982  $ CALL dscal( n-ki+1, scale, work( ki+iv*n ), 1 )
983  work( j +iv*n ) = x( 1, 1 )
984  work( j+1+iv*n ) = x( 2, 1 )
985 *
986  vmax = max( abs( work( j +iv*n ) ),
987  $ abs( work( j+1+iv*n ) ), vmax )
988  vcrit = bignum / vmax
989 *
990  END IF
991  170 CONTINUE
992 *
993 * Copy the vector x or Q*x to VL and normalize.
994 *
995  IF( .NOT.over ) THEN
996 * ------------------------------
997 * no back-transform: copy x to VL and normalize.
998  CALL dcopy( n-ki+1, work( ki + iv*n ), 1,
999  $ vl( ki, is ), 1 )
1000 *
1001  ii = idamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
1002  remax = one / abs( vl( ii, is ) )
1003  CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1004 *
1005  DO 180 k = 1, ki - 1
1006  vl( k, is ) = zero
1007  180 CONTINUE
1008 *
1009  ELSE IF( nb.EQ.1 ) THEN
1010 * ------------------------------
1011 * version 1: back-transform each vector with GEMV, Q*x.
1012  IF( ki.LT.n )
1013  $ CALL dgemv( 'N', n, n-ki, one,
1014  $ vl( 1, ki+1 ), ldvl,
1015  $ work( ki+1 + iv*n ), 1,
1016  $ work( ki + iv*n ), vl( 1, ki ), 1 )
1017 *
1018  ii = idamax( n, vl( 1, ki ), 1 )
1019  remax = one / abs( vl( ii, ki ) )
1020  CALL dscal( n, remax, vl( 1, ki ), 1 )
1021 *
1022  ELSE
1023 * ------------------------------
1024 * version 2: back-transform block of vectors with GEMM
1025 * zero out above vector
1026 * could go from KI-NV+1 to KI-1
1027  DO k = 1, ki - 1
1028  work( k + iv*n ) = zero
1029  END DO
1030  iscomplex( iv ) = ip
1031 * back-transform and normalization is done below
1032  END IF
1033  ELSE
1034 *
1035 * --------------------------------------------------------
1036 * Complex left eigenvector.
1037 *
1038 * Initial solve:
1039 * [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
1040 * [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
1041 *
1042  IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) ) THEN
1043  work( ki + (iv )*n ) = wi / t( ki, ki+1 )
1044  work( ki+1 + (iv+1)*n ) = one
1045  ELSE
1046  work( ki + (iv )*n ) = one
1047  work( ki+1 + (iv+1)*n ) = -wi / t( ki+1, ki )
1048  END IF
1049  work( ki+1 + (iv )*n ) = zero
1050  work( ki + (iv+1)*n ) = zero
1051 *
1052 * Form right-hand side.
1053 *
1054  DO 190 k = ki + 2, n
1055  work( k+(iv )*n ) = -work( ki +(iv )*n )*t(ki, k)
1056  work( k+(iv+1)*n ) = -work( ki+1+(iv+1)*n )*t(ki+1,k)
1057  190 CONTINUE
1058 *
1059 * Solve transposed quasi-triangular system:
1060 * [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
1061 *
1062  vmax = one
1063  vcrit = bignum
1064 *
1065  jnxt = ki + 2
1066  DO 200 j = ki + 2, n
1067  IF( j.LT.jnxt )
1068  $ GO TO 200
1069  j1 = j
1070  j2 = j
1071  jnxt = j + 1
1072  IF( j.LT.n ) THEN
1073  IF( t( j+1, j ).NE.zero ) THEN
1074  j2 = j + 1
1075  jnxt = j + 2
1076  END IF
1077  END IF
1078 *
1079  IF( j1.EQ.j2 ) THEN
1080 *
1081 * 1-by-1 diagonal block
1082 *
1083 * Scale if necessary to avoid overflow when
1084 * forming the right-hand side elements.
1085 *
1086  IF( work( j ).GT.vcrit ) THEN
1087  rec = one / vmax
1088  CALL dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1089  CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1090  vmax = one
1091  vcrit = bignum
1092  END IF
1093 *
1094  work( j+(iv )*n ) = work( j+(iv)*n ) -
1095  $ ddot( j-ki-2, t( ki+2, j ), 1,
1096  $ work( ki+2+(iv)*n ), 1 )
1097  work( j+(iv+1)*n ) = work( j+(iv+1)*n ) -
1098  $ ddot( j-ki-2, t( ki+2, j ), 1,
1099  $ work( ki+2+(iv+1)*n ), 1 )
1100 *
1101 * Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
1102 *
1103  CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
1104  $ ldt, one, one, work( j+iv*n ), n, wr,
1105  $ -wi, x, 2, scale, xnorm, ierr )
1106 *
1107 * Scale if necessary
1108 *
1109  IF( scale.NE.one ) THEN
1110  CALL dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1111  CALL dscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1112  END IF
1113  work( j+(iv )*n ) = x( 1, 1 )
1114  work( j+(iv+1)*n ) = x( 1, 2 )
1115  vmax = max( abs( work( j+(iv )*n ) ),
1116  $ abs( work( j+(iv+1)*n ) ), vmax )
1117  vcrit = bignum / vmax
1118 *
1119  ELSE
1120 *
1121 * 2-by-2 diagonal block
1122 *
1123 * Scale if necessary to avoid overflow when forming
1124 * the right-hand side elements.
1125 *
1126  beta = max( work( j ), work( j+1 ) )
1127  IF( beta.GT.vcrit ) THEN
1128  rec = one / vmax
1129  CALL dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1130  CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1131  vmax = one
1132  vcrit = bignum
1133  END IF
1134 *
1135  work( j +(iv )*n ) = work( j+(iv)*n ) -
1136  $ ddot( j-ki-2, t( ki+2, j ), 1,
1137  $ work( ki+2+(iv)*n ), 1 )
1138 *
1139  work( j +(iv+1)*n ) = work( j+(iv+1)*n ) -
1140  $ ddot( j-ki-2, t( ki+2, j ), 1,
1141  $ work( ki+2+(iv+1)*n ), 1 )
1142 *
1143  work( j+1+(iv )*n ) = work( j+1+(iv)*n ) -
1144  $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1145  $ work( ki+2+(iv)*n ), 1 )
1146 *
1147  work( j+1+(iv+1)*n ) = work( j+1+(iv+1)*n ) -
1148  $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1149  $ work( ki+2+(iv+1)*n ), 1 )
1150 *
1151 * Solve 2-by-2 complex linear equation
1152 * [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
1153 * [ (T(j+1,j) T(j+1,j+1)) ]
1154 *
1155  CALL dlaln2( .true., 2, 2, smin, one, t( j, j ),
1156  $ ldt, one, one, work( j+iv*n ), n, wr,
1157  $ -wi, x, 2, scale, xnorm, ierr )
1158 *
1159 * Scale if necessary
1160 *
1161  IF( scale.NE.one ) THEN
1162  CALL dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1163  CALL dscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1164  END IF
1165  work( j +(iv )*n ) = x( 1, 1 )
1166  work( j +(iv+1)*n ) = x( 1, 2 )
1167  work( j+1+(iv )*n ) = x( 2, 1 )
1168  work( j+1+(iv+1)*n ) = x( 2, 2 )
1169  vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1170  $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ),
1171  $ vmax )
1172  vcrit = bignum / vmax
1173 *
1174  END IF
1175  200 CONTINUE
1176 *
1177 * Copy the vector x or Q*x to VL and normalize.
1178 *
1179  IF( .NOT.over ) THEN
1180 * ------------------------------
1181 * no back-transform: copy x to VL and normalize.
1182  CALL dcopy( n-ki+1, work( ki + (iv )*n ), 1,
1183  $ vl( ki, is ), 1 )
1184  CALL dcopy( n-ki+1, work( ki + (iv+1)*n ), 1,
1185  $ vl( ki, is+1 ), 1 )
1186 *
1187  emax = zero
1188  DO 220 k = ki, n
1189  emax = max( emax, abs( vl( k, is ) )+
1190  $ abs( vl( k, is+1 ) ) )
1191  220 CONTINUE
1192  remax = one / emax
1193  CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1194  CALL dscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1195 *
1196  DO 230 k = 1, ki - 1
1197  vl( k, is ) = zero
1198  vl( k, is+1 ) = zero
1199  230 CONTINUE
1200 *
1201  ELSE IF( nb.EQ.1 ) THEN
1202 * ------------------------------
1203 * version 1: back-transform each vector with GEMV, Q*x.
1204  IF( ki.LT.n-1 ) THEN
1205  CALL dgemv( 'N', n, n-ki-1, one,
1206  $ vl( 1, ki+2 ), ldvl,
1207  $ work( ki+2 + (iv)*n ), 1,
1208  $ work( ki + (iv)*n ),
1209  $ vl( 1, ki ), 1 )
1210  CALL dgemv( 'N', n, n-ki-1, one,
1211  $ vl( 1, ki+2 ), ldvl,
1212  $ work( ki+2 + (iv+1)*n ), 1,
1213  $ work( ki+1 + (iv+1)*n ),
1214  $ vl( 1, ki+1 ), 1 )
1215  ELSE
1216  CALL dscal( n, work(ki+ (iv )*n), vl(1, ki ), 1)
1217  CALL dscal( n, work(ki+1+(iv+1)*n), vl(1, ki+1), 1)
1218  END IF
1219 *
1220  emax = zero
1221  DO 240 k = 1, n
1222  emax = max( emax, abs( vl( k, ki ) )+
1223  $ abs( vl( k, ki+1 ) ) )
1224  240 CONTINUE
1225  remax = one / emax
1226  CALL dscal( n, remax, vl( 1, ki ), 1 )
1227  CALL dscal( n, remax, vl( 1, ki+1 ), 1 )
1228 *
1229  ELSE
1230 * ------------------------------
1231 * version 2: back-transform block of vectors with GEMM
1232 * zero out above vector
1233 * could go from KI-NV+1 to KI-1
1234  DO k = 1, ki - 1
1235  work( k + (iv )*n ) = zero
1236  work( k + (iv+1)*n ) = zero
1237  END DO
1238  iscomplex( iv ) = ip
1239  iscomplex( iv+1 ) = -ip
1240  iv = iv + 1
1241 * back-transform and normalization is done below
1242  END IF
1243  END IF
1244 
1245  IF( nb.GT.1 ) THEN
1246 * --------------------------------------------------------
1247 * Blocked version of back-transform
1248 * For complex case, KI2 includes both vectors (KI and KI+1)
1249  IF( ip.EQ.0 ) THEN
1250  ki2 = ki
1251  ELSE
1252  ki2 = ki + 1
1253  END IF
1254 
1255 * Columns 1:IV of work are valid vectors.
1256 * When the number of vectors stored reaches NB-1 or NB,
1257 * or if this was last vector, do the GEMM
1258  IF( (iv.GE.nb-1) .OR. (ki2.EQ.n) ) THEN
1259  CALL dgemm( 'N', 'N', n, iv, n-ki2+iv, one,
1260  $ vl( 1, ki2-iv+1 ), ldvl,
1261  $ work( ki2-iv+1 + (1)*n ), n,
1262  $ zero,
1263  $ work( 1 + (nb+1)*n ), n )
1264 * normalize vectors
1265  DO k = 1, iv
1266  IF( iscomplex(k).EQ.0) THEN
1267 * real eigenvector
1268  ii = idamax( n, work( 1 + (nb+k)*n ), 1 )
1269  remax = one / abs( work( ii + (nb+k)*n ) )
1270  ELSE IF( iscomplex(k).EQ.1) THEN
1271 * first eigenvector of conjugate pair
1272  emax = zero
1273  DO ii = 1, n
1274  emax = max( emax,
1275  $ abs( work( ii + (nb+k )*n ) )+
1276  $ abs( work( ii + (nb+k+1)*n ) ) )
1277  END DO
1278  remax = one / emax
1279 * else if ISCOMPLEX(K).EQ.-1
1280 * second eigenvector of conjugate pair
1281 * reuse same REMAX as previous K
1282  END IF
1283  CALL dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
1284  END DO
1285  CALL dlacpy( 'F', n, iv,
1286  $ work( 1 + (nb+1)*n ), n,
1287  $ vl( 1, ki2-iv+1 ), ldvl )
1288  iv = 1
1289  ELSE
1290  iv = iv + 1
1291  END IF
1292  END IF ! blocked back-transform
1293 *
1294  is = is + 1
1295  IF( ip.NE.0 )
1296  $ is = is + 1
1297  260 CONTINUE
1298  END IF
1299 *
1300  RETURN
1301 *
1302 * End of DTREVC3
1303 *
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
idamax
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
dlaln2
subroutine dlaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition: dlaln2.f:220
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76

Here is the call graph for this function:

Here is the caller graph for this function: