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