LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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,
278 $ slacpy, sgemm, slabad, slaset
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*
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
slaset
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
slacpy
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
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
slaln2
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
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
sdot
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
slamch
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: