LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ strevc3()

subroutine strevc3 ( character  SIDE,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
real, dimension( ldt, * )  T,
integer  LDT,
real, dimension( ldvl, * )  VL,
integer  LDVL,
real, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

STREVC3

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

Purpose:
 STREVC3 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 SHSEQR.

 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 REAL 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 REAL 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 SHSEQR).
          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 REAL 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 SHSEQR).
          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 REAL 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.
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 235 of file strevc3.f.

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