LAPACK  3.6.1
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**T)*T = w*(y**T)

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

Here is the call graph for this function:

Here is the caller graph for this function: