+C> \brief \b DLAMCH
+C>\details
+C> \b Purpose:
+C>\verbatim
+C>
+C> DLAMCH determines double precision machine parameters.
+C>
+C>\endverbatim
+C> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
+C> \date November 2011
+C> \ingroup auxOTHERauxiliary
+C>
+C> \param[in] CMACH
+C> \verbatim
+C> Specifies the value to be returned by DLAMCH:
+C> = 'E' or 'e', DLAMCH := eps
+C> = 'S' or 's , DLAMCH := sfmin
+C> = 'B' or 'b', DLAMCH := base
+C> = 'P' or 'p', DLAMCH := eps*base
+C> = 'N' or 'n', DLAMCH := t
+C> = 'R' or 'r', DLAMCH := rnd
+C> = 'M' or 'm', DLAMCH := emin
+C> = 'U' or 'u', DLAMCH := rmin
+C> = 'L' or 'l', DLAMCH := emax
+C> = 'O' or 'o', DLAMCH := rmax
+C> where
+C> eps = relative machine precision
+C> sfmin = safe minimum, such that 1/sfmin does not overflow
+C> base = base of the machine
+C> prec = eps*base
+C> t = number of (base) digits in the mantissa
+C> rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
+C> emin = minimum exponent before (gradual) underflow
+C> rmin = underflow threshold - base**(emin-1)
+C> emax = largest exponent before overflow
+C> rmax = overflow threshold - (base**emax)*(1-eps)
+C> \endverbatim
+C>
DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
-*
-* -- LAPACK auxiliary routine (version 3.3.0) --
-* -- LAPACK is a software package provided by Univ. of Tennessee, --
-* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* Based on LAPACK DLAMCH but with Fortran 95 query functions
-* See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
-* and http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
-* July 2010
-*
-* .. Scalar Arguments ..
+C
+C -- LAPACK auxiliary routine (version 3.3.0) --
+C -- LAPACK is a software package provided by Univ. of Tennessee, --
+C -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+C November 2011
+C
+C .. Scalar Arguments ..
CHARACTER CMACH
-* ..
-*
-* Purpose
-* =======
-*
-* DLAMCH determines double precision machine parameters.
-*
-* Arguments
-* =========
-*
-* CMACH (input) CHARACTER*1
-* Specifies the value to be returned by DLAMCH:
-* = 'E' or 'e', DLAMCH := eps
-* = 'S' or 's , DLAMCH := sfmin
-* = 'B' or 'b', DLAMCH := base
-* = 'P' or 'p', DLAMCH := eps*base
-* = 'N' or 'n', DLAMCH := t
-* = 'R' or 'r', DLAMCH := rnd
-* = 'M' or 'm', DLAMCH := emin
-* = 'U' or 'u', DLAMCH := rmin
-* = 'L' or 'l', DLAMCH := emax
-* = 'O' or 'o', DLAMCH := rmax
-*
-* where
-*
-* eps = relative machine precision
-* sfmin = safe minimum, such that 1/sfmin does not overflow
-* base = base of the machine
-* prec = eps*base
-* t = number of (base) digits in the mantissa
-* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
-* emin = minimum exponent before (gradual) underflow
-* rmin = underflow threshold - base**(emin-1)
-* emax = largest exponent before overflow
-* rmax = overflow threshold - (base**emax)*(1-eps)
-*
-* =====================================================================
-*
-* .. Parameters ..
+C ..
+C
+C .. Scalar Arguments ..
+ DOUBLE PRECISION A, B
+C ..
+C
+C =====================================================================
+C
+C .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
-* ..
-* .. Local Scalars ..
+C ..
+C .. Local Scalars ..
DOUBLE PRECISION RND, EPS, SFMIN, SMALL, RMACH
-* ..
-* .. External Functions ..
+C ..
+C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
-* ..
-* .. Intrinsic Functions ..
+C ..
+C .. Intrinsic Functions ..
INTRINSIC DIGITS, EPSILON, HUGE, MAXEXPONENT,
$ MINEXPONENT, RADIX, TINY
-* ..
-* .. Executable Statements ..
-*
-*
-* Assume rounding, not chopping. Always.
-*
+C ..
+C .. Executable Statements ..
+C
+C
+C Assume rounding, not chopping. Always.
+C
RND = ONE
-*
+C
IF( ONE.EQ.RND ) THEN
EPS = EPSILON(ZERO) * 0.5
ELSE
EPS = EPSILON(ZERO)
END IF
-*
+C
IF( LSAME( CMACH, 'E' ) ) THEN
RMACH = EPS
ELSE IF( LSAME( CMACH, 'S' ) ) THEN
SFMIN = TINY(ZERO)
SMALL = ONE / HUGE(ZERO)
IF( SMALL.GE.SFMIN ) THEN
-*
-* Use SMALL plus a bit, to avoid the possibility of rounding
-* causing overflow when computing 1/sfmin.
-*
+C
+C Use SMALL plus a bit, to avoid the possibility of rounding
+C causing overflow when computing 1/sfmin.
+C
SFMIN = SMALL*( ONE+EPS )
END IF
RMACH = SFMIN
ELSE
RMACH = ZERO
END IF
-*
+C
DLAMCH = RMACH
RETURN
-*
-* End of DLAMCH
-*
+C
+C End of DLAMCH
+C
END
-************************************************************************
-*
+C***********************************************************************
+C> \brief \b DLAMC3
+C>\details
+C> \b Purpose:
+C>\verbatim
+C>
+C> DLAMC3 is intended to force A and B to be stored prior to doing
+C> the addition of A and B , for use in situations where optimizers
+C> might hold one of these in a register.
+C>
+C>\endverbatim
+C> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
+C> \date November 2011
+C> \ingroup auxOTHERauxiliary
+C>
+C> \param[in] A
+C> \verbatim
+C> A is a DOUBLE PRECISION
+C> \endverbatim
+C>
+C> \param[in] B
+C> \verbatim
+C> B is a DOUBLE PRECISION
+C> The values A and B.
+C> \endverbatim
+C>
DOUBLE PRECISION FUNCTION DLAMC3( A, B )
-*
-* -- LAPACK auxiliary routine (version 3.3.0) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2010
-*
-* .. Scalar Arguments ..
+C
+C -- LAPACK auxiliary routine (version 3.3.0) --
+C Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+C November 2010
+C
+C .. Scalar Arguments ..
DOUBLE PRECISION A, B
-* ..
-*
-* Purpose
-* =======
-*
-* DLAMC3 is intended to force A and B to be stored prior to doing
-* the addition of A and B , for use in situations where optimizers
-* might hold one of these in a register.
-*
-* Arguments
-* =========
-*
-* A (input) DOUBLE PRECISION
-* B (input) DOUBLE PRECISION
-* The values A and B.
-*
-* =====================================================================
-*
-* .. Executable Statements ..
-*
+C ..
+C =====================================================================
+C
+C .. Executable Statements ..
+C
DLAMC3 = A + B
-*
+C
RETURN
-*
-* End of DLAMC3
-*
+C
+C End of DLAMC3
+C
END
-*
-************************************************************************
+C
+C***********************************************************************
- PROGRAM TEST3
+ PROGRAM DLAMCHTST
*
* -- LAPACK test routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
+* =====================================================================
+*
* .. Local Scalars ..
DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMAX, RMIN, RND,
$ SFMIN, T
- PROGRAM TEST5
+ PROGRAM DSECNDTST
*
* -- LAPACK test routine (version 3.3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* -- April 2011 --
*
+* =====================================================================
+*
* .. Parameters ..
INTEGER NMAX, ITS
PARAMETER ( NMAX = 1000, ITS = 50000 )
* =========
*
* CA (input) CHARACTER*1
+*
* CB (input) CHARACTER*1
* CA and CB specify the single characters to be compared.
*
- PROGRAM TEST1
+ PROGRAM LSAMETST
*
* -- LAPACK test routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
- REAL FUNCTION SECOND( )
+ REAL FUNCTION SECOND( )
*
* -- LAPACK auxiliary routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
- REAL FUNCTION SECOND( )
+ REAL FUNCTION SECOND( )
*
* -- LAPACK auxiliary routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
- REAL FUNCTION SECOND( )
+ REAL FUNCTION SECOND( )
*
* -- LAPACK auxiliary routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
- REAL FUNCTION SECOND( )
+ REAL FUNCTION SECOND( )
*
* -- LAPACK auxiliary routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
- REAL FUNCTION SECOND( )
+ REAL FUNCTION SECOND( )
*
* -- LAPACK auxiliary routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
- PROGRAM TEST4
+ PROGRAM SECONDTST
*
* -- LAPACK test routine (version 3.3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* -- April 2011 --
*
+* =====================================================================
+*
* .. Parameters ..
INTEGER NMAX, ITS
PARAMETER ( NMAX = 1000, ITS = 50000 )
+C> \brief \b SLAMCH
+C>\details
+C> \b Purpose:
+C>\verbatim
+C>
+C> SLAMCH determines single precision machine parameters.
+C>
+C>\endverbatim
+C> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
+C> \date November 2011
+C> \ingroup auxOTHERauxiliary
+C>
+C> \param[in] CMACH
+C> \verbatim
+C> Specifies the value to be returned by SLAMCH:
+C> = 'E' or 'e', SLAMCH := eps
+C> = 'S' or 's , SLAMCH := sfmin
+C> = 'B' or 'b', SLAMCH := base
+C> = 'P' or 'p', SLAMCH := eps*base
+C> = 'N' or 'n', SLAMCH := t
+C> = 'R' or 'r', SLAMCH := rnd
+C> = 'M' or 'm', SLAMCH := emin
+C> = 'U' or 'u', SLAMCH := rmin
+C> = 'L' or 'l', SLAMCH := emax
+C> = 'O' or 'o', SLAMCH := rmax
+C> where
+C> eps = relative machine precision
+C> sfmin = safe minimum, such that 1/sfmin does not overflow
+C> base = base of the machine
+C> prec = eps*base
+C> t = number of (base) digits in the mantissa
+C> rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
+C> emin = minimum exponent before (gradual) underflow
+C> rmin = underflow threshold - base**(emin-1)
+C> emax = largest exponent before overflow
+C> rmax = overflow threshold - (base**emax)*(1-eps)
+C> \endverbatim
+C>
REAL FUNCTION SLAMCH( CMACH )
-*
-* -- LAPACK auxiliary routine (version 3.3.0) --
-* -- LAPACK is a software package provided by Univ. of Tennessee, --
-* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* Based on LAPACK DLAMCH but with Fortran 95 query functions
-* See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
-* and http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
-* July 2010
-*
-* .. Scalar Arguments ..
+C
+C -- LAPACK auxiliary routine (version 3.3.0) --
+C -- LAPACK is a software package provided by Univ. of Tennessee, --
+C -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+C November 2011
+C
+C .. Scalar Arguments ..
CHARACTER CMACH
-* ..
-*
-* Purpose
-* =======
-*
-* SLAMCH determines single precision machine parameters.
-*
-* Arguments
-* =========
-*
-* CMACH (input) CHARACTER*1
-* Specifies the value to be returned by SLAMCH:
-* = 'E' or 'e', SLAMCH := eps
-* = 'S' or 's , SLAMCH := sfmin
-* = 'B' or 'b', SLAMCH := base
-* = 'P' or 'p', SLAMCH := eps*base
-* = 'N' or 'n', SLAMCH := t
-* = 'R' or 'r', SLAMCH := rnd
-* = 'M' or 'm', SLAMCH := emin
-* = 'U' or 'u', SLAMCH := rmin
-* = 'L' or 'l', SLAMCH := emax
-* = 'O' or 'o', SLAMCH := rmax
-*
-* where
-*
-* eps = relative machine precision
-* sfmin = safe minimum, such that 1/sfmin does not overflow
-* base = base of the machine
-* prec = eps*base
-* t = number of (base) digits in the mantissa
-* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
-* emin = minimum exponent before (gradual) underflow
-* rmin = underflow threshold - base**(emin-1)
-* emax = largest exponent before overflow
-* rmax = overflow threshold - (base**emax)*(1-eps)
-*
-* =====================================================================
-*
-* .. Parameters ..
+C ..
+C
+C .. Scalar Arguments ..
+ REAL A, B
+C ..
+C
+C =====================================================================
+C
+C .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
-* ..
-* .. Local Scalars ..
+C ..
+C .. Local Scalars ..
REAL RND, EPS, SFMIN, SMALL, RMACH
-* ..
-* .. External Functions ..
+C ..
+C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
-* ..
-* .. Intrinsic Functions ..
+C ..
+C .. Intrinsic Functions ..
INTRINSIC DIGITS, EPSILON, HUGE, MAXEXPONENT,
$ MINEXPONENT, RADIX, TINY
-* ..
-* .. Executable Statements ..
-*
-*
-* Assume rounding, not chopping. Always.
-*
+C ..
+C .. Executable Statements ..
+C
+C
+C Assume rounding, not chopping. Always.
+C
RND = ONE
-*
+C
IF( ONE.EQ.RND ) THEN
EPS = EPSILON(ZERO) * 0.5
ELSE
EPS = EPSILON(ZERO)
END IF
-*
+C
IF( LSAME( CMACH, 'E' ) ) THEN
RMACH = EPS
ELSE IF( LSAME( CMACH, 'S' ) ) THEN
SFMIN = TINY(ZERO)
SMALL = ONE / HUGE(ZERO)
IF( SMALL.GE.SFMIN ) THEN
-*
-* Use SMALL plus a bit, to avoid the possibility of rounding
-* causing overflow when computing 1/sfmin.
-*
+C
+C Use SMALL plus a bit, to avoid the possibility of rounding
+C causing overflow when computing 1/sfmin.
+C
SFMIN = SMALL*( ONE+EPS )
END IF
RMACH = SFMIN
ELSE
RMACH = ZERO
END IF
-*
+C
SLAMCH = RMACH
RETURN
-*
-* End of SLAMCH
-*
+C
+C End of SLAMCH
+C
END
-************************************************************************
-*
+C***********************************************************************
+C> \brief \b SLAMC3
+C>\details
+C> \b Purpose:
+C>\verbatim
+C>
+C> SLAMC3 is intended to force A and B to be stored prior to doing
+C> the addition of A and B , for use in situations where optimizers
+C> might hold one of these in a register.
+C>
+C>\endverbatim
+C> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
+C> \date November 2011
+C> \ingroup auxOTHERauxiliary
+C>
+C> \param[in] A
+C> \verbatim
+C> \endverbatim
+C>
+C> \param[in] B
+C> \verbatim
+C> The values A and B.
+C> \endverbatim
+C>
+C
REAL FUNCTION SLAMC3( A, B )
-*
-* -- LAPACK auxiliary routine (version 3.3.0) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2010
-*
-* .. Scalar Arguments ..
+C
+C -- LAPACK auxiliary routine (version 3.3.0) --
+C Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+C November 2010
+C
+C .. Scalar Arguments ..
REAL A, B
-* ..
-*
-* Purpose
-* =======
-*
-* SLAMC3 is intended to force A and B to be stored prior to doing
-* the addition of A and B , for use in situations where optimizers
-* might hold one of these in a register.
-*
-* Arguments
-* =========
-*
-* A (input) REAL
-* B (input) REAL
-* The values A and B.
-*
-* =====================================================================
-*
-* .. Executable Statements ..
-*
+C ..
+C =====================================================================
+C
+C .. Executable Statements ..
+C
SLAMC3 = A + B
-*
+C
RETURN
-*
-* End of SLAMC3
-*
+C
+C End of SLAMC3
+C
END
-*
-************************************************************************
+C
+C***********************************************************************
- PROGRAM TEST2
+ PROGRAM SLAMCHTST
*
* -- LAPACK test routine (version 3.2) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
+* =====================================================================
+*
* .. Local Scalars ..
REAL BASE, EMAX, EMIN, EPS, RMAX, RMIN, RND, SFMIN,
$ T, PREC
- PROGRAM MAIN
-*
-* -- LAPACK test routine (version 3.2) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2006
-*
-* .. External Functions ..
+C> \brief \b TSTIEE
+C>\details
+C> \b Purpose:
+C>\verbatim
+C>
+C> TEST IEEE
+C>
+C>\endverbatim
+C> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
+C> \date November 2011
+C> \ingroup auxOTHERauxiliary
+
+ PROGRAM TSTIEE
+C
+C -- LAPACK test routine (version 3.2) --
+C Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+C November 2006
+C
+C .. External Functions ..
INTEGER ILAENV
EXTERNAL ILAENV
-* ..
-* .. Local Scalars ..
+C ..
+C .. Local Scalars ..
INTEGER IEEEOK
-* ..
-* .. Executable Statements ..
-*
+C ..
+C .. Executable Statements ..
+C
WRITE( 6, FMT = * )
$ 'We are about to check whether infinity arithmetic'
WRITE( 6, FMT = * )'can be trusted. If this test hangs, set'
WRITE( 6, FMT = * )
$ 'ILAENV = 0 for ISPEC = 10 in LAPACK/SRC/ilaenv.f'
-*
+C
IEEEOK = ILAENV( 10, 'ILAENV', 'N', 1, 2, 3, 4 )
WRITE( 6, FMT = * )
-*
+C
IF( IEEEOK.EQ.0 ) THEN
WRITE( 6, FMT = * )
$ 'Infinity arithmetic did not perform per the ieee spec'
$ 'guarantee that infinity arithmetic meets the',
$ ' ieee spec.'
END IF
-*
+C
WRITE( 6, FMT = * )
WRITE( 6, FMT = * )
$ 'We are about to check whether NaN arithmetic'
WRITE( 6, FMT = * )
$ 'ILAENV = 0 for ISPEC = 11 in LAPACK/SRC/ilaenv.f'
IEEEOK = ILAENV( 11, 'ILAENV', 'N', 1, 2, 3, 4 )
-*
+C
WRITE( 6, FMT = * )
IF( IEEEOK.EQ.0 ) THEN
WRITE( 6, FMT = * )
$ ' ieee spec.'
END IF
WRITE( 6, FMT = * )
-*
+C
END
INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
$ N4 )
-*
-* -- LAPACK auxiliary routine (version 3.2) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2006
-*
-* .. Scalar Arguments ..
+C
+C -- LAPACK auxiliary routine (version 3.2) --
+C Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+C November 2006
+C
+C .. Scalar Arguments ..
CHARACTER*( * ) NAME, OPTS
INTEGER ISPEC, N1, N2, N3, N4
-* ..
-*
-* Purpose
-* =======
-*
-* ILAENV is called from the LAPACK routines to choose problem-dependent
-* parameters for the local environment. See ISPEC for a description of
-* the parameters.
-*
-* This version provides a set of parameters which should give good,
-* but not optimal, performance on many of the currently available
-* computers. Users are encouraged to modify this subroutine to set
-* the tuning parameters for their particular machine using the option
-* and problem size information in the arguments.
-*
-* This routine will not function correctly if it is converted to all
-* lower case. Converting it to all upper case is allowed.
-*
-* Arguments
-* =========
-*
-* ISPEC (input) INTEGER
-* Specifies the parameter to be returned as the value of
-* ILAENV.
-* = 1: the optimal blocksize; if this value is 1, an unblocked
-* algorithm will give the best performance.
-* = 2: the minimum block size for which the block routine
-* should be used; if the usable block size is less than
-* this value, an unblocked routine should be used.
-* = 3: the crossover point (in a block routine, for N less
-* than this value, an unblocked routine should be used)
-* = 4: the number of shifts, used in the nonsymmetric
-* eigenvalue routines
-* = 5: the minimum column dimension for blocking to be used;
-* rectangular blocks must have dimension at least k by m,
-* where k is given by ILAENV(2,...) and m by ILAENV(5,...)
-* = 6: the crossover point for the SVD (when reducing an m by n
-* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
-* this value, a QR factorization is used first to reduce
-* the matrix to a triangular form.)
-* = 7: the number of processors
-* = 8: the crossover point for the multishift QR and QZ methods
-* for nonsymmetric eigenvalue problems.
-* = 9: maximum size of the subproblems at the bottom of the
-* computation tree in the divide-and-conquer algorithm
-* (used by xGELSD and xGESDD)
-* =10: ieee NaN arithmetic can be trusted not to trap
-* =11: infinity arithmetic can be trusted not to trap
-*
-* NAME (input) CHARACTER*(*)
-* The name of the calling subroutine, in either upper case or
-* lower case.
-*
-* OPTS (input) CHARACTER*(*)
-* The character options to the subroutine NAME, concatenated
-* into a single character string. For example, UPLO = 'U',
-* TRANS = 'T', and DIAG = 'N' for a triangular routine would
-* be specified as OPTS = 'UTN'.
-*
-* N1 (input) INTEGER
-* N2 (input) INTEGER
-* N3 (input) INTEGER
-* N4 (input) INTEGER
-* Problem dimensions for the subroutine NAME; these may not all
-* be required.
-*
-* (ILAENV) (output) INTEGER
-* >= 0: the value of the parameter specified by ISPEC
-* < 0: if ILAENV = -k, the k-th argument had an illegal value.
-*
-* Further Details
-* ===============
-*
-* The following conventions have been used when calling ILAENV from the
-* LAPACK routines:
-* 1) OPTS is a concatenation of all of the character options to
-* subroutine NAME, in the same order that they appear in the
-* argument list for NAME, even if they are not used in determining
-* the value of the parameter specified by ISPEC.
-* 2) The problem dimensions N1, N2, N3, N4 are specified in the order
-* that they appear in the argument list for NAME. N1 is used
-* first, N2 second, and so on, and unused problem dimensions are
-* passed a value of -1.
-* 3) The parameter value returned by ILAENV is checked for validity in
-* the calling subroutine. For example, ILAENV is used to retrieve
-* the optimal blocksize for STRTRI as follows:
-*
-* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
-* IF( NB.LE.1 ) NB = MAX( 1, N )
-*
-* =====================================================================
-*
-* .. Local Scalars ..
+C ..
+C
+C Purpose
+C =======
+C
+C ILAENV is called from the LAPACK routines to choose problem-dependent
+C parameters for the local environment. See ISPEC for a description of
+C the parameters.
+C
+C This version provides a set of parameters which should give good,
+C but not optimal, performance on many of the currently available
+C computers. Users are encouraged to modify this subroutine to set
+C the tuning parameters for their particular machine using the option
+C and problem size information in the arguments.
+C
+C This routine will not function correctly if it is converted to all
+C lower case. Converting it to all upper case is allowed.
+C
+C Arguments
+C =========
+C
+C ISPEC (input) INTEGER
+C Specifies the parameter to be returned as the value of
+C ILAENV.
+C = 1: the optimal blocksize; if this value is 1, an unblocked
+C algorithm will give the best performance.
+C = 2: the minimum block size for which the block routine
+C should be used; if the usable block size is less than
+C this value, an unblocked routine should be used.
+C = 3: the crossover point (in a block routine, for N less
+C than this value, an unblocked routine should be used)
+C = 4: the number of shifts, used in the nonsymmetric
+C eigenvalue routines
+C = 5: the minimum column dimension for blocking to be used;
+C rectangular blocks must have dimension at least k by m,
+C where k is given by ILAENV(2,...) and m by ILAENV(5,...)
+C = 6: the crossover point for the SVD (when reducing an m by n
+C matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
+C this value, a QR factorization is used first to reduce
+C the matrix to a triangular form.)
+C = 7: the number of processors
+C = 8: the crossover point for the multishift QR and QZ methods
+C for nonsymmetric eigenvalue problems.
+C = 9: maximum size of the subproblems at the bottom of the
+C computation tree in the divide-and-conquer algorithm
+C (used by xGELSD and xGESDD)
+C =10: ieee NaN arithmetic can be trusted not to trap
+C =11: infinity arithmetic can be trusted not to trap
+C
+C NAME (input) CHARACTER*(*)
+C The name of the calling subroutine, in either upper case or
+C lower case.
+C
+C OPTS (input) CHARACTER*(*)
+C The character options to the subroutine NAME, concatenated
+C into a single character string. For example, UPLO = 'U',
+C TRANS = 'T', and DIAG = 'N' for a triangular routine would
+C be specified as OPTS = 'UTN'.
+C
+C N1 (input) INTEGER
+C N2 (input) INTEGER
+C N3 (input) INTEGER
+C N4 (input) INTEGER
+C Problem dimensions for the subroutine NAME; these may not all
+C be required.
+C
+C (ILAENV) (output) INTEGER
+C >= 0: the value of the parameter specified by ISPEC
+C < 0: if ILAENV = -k, the k-th argument had an illegal value.
+C
+C Further Details
+C ===============
+C
+C The following conventions have been used when calling ILAENV from the
+C LAPACK routines:
+C 1) OPTS is a concatenation of all of the character options to
+C subroutine NAME, in the same order that they appear in the
+C argument list for NAME, even if they are not used in determining
+C the value of the parameter specified by ISPEC.
+C 2) The problem dimensions N1, N2, N3, N4 are specified in the order
+C that they appear in the argument list for NAME. N1 is used
+C first, N2 second, and so on, and unused problem dimensions are
+C passed a value of -1.
+C 3) The parameter value returned by ILAENV is checked for validity in
+C the calling subroutine. For example, ILAENV is used to retrieve
+C the optimal blocksize for STRTRI as follows:
+C
+C NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
+C IF( NB.LE.1 ) NB = MAX( 1, N )
+C
+C =====================================================================
+C
+C .. Local Scalars ..
LOGICAL CNAME, SNAME
CHARACTER*1 C1
CHARACTER*2 C2, C4
CHARACTER*3 C3
CHARACTER*6 SUBNAM
INTEGER I, IC, IZ, NB, NBMIN, NX
-* ..
-* .. Intrinsic Functions ..
+C ..
+C .. Intrinsic Functions ..
INTRINSIC CHAR, ICHAR, INT, MIN, REAL
-* ..
-* .. External Functions ..
+C ..
+C .. External Functions ..
INTEGER IEEECK
EXTERNAL IEEECK
-* ..
-* .. Executable Statements ..
-*
+C ..
+C .. Executable Statements ..
+C
GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000,
$ 1100 ) ISPEC
-*
-* Invalid value for ISPEC
-*
+C
+C Invalid value for ISPEC
+C
ILAENV = -1
RETURN
-*
+C
100 CONTINUE
-*
-* Convert NAME to upper case if the first character is lower case.
-*
+C
+C Convert NAME to upper case if the first character is lower case.
+C
ILAENV = 1
SUBNAM = NAME
IC = ICHAR( SUBNAM( 1:1 ) )
IZ = ICHAR( 'Z' )
IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
-*
-* ASCII character set
-*
+C
+C ASCII character set
+C
IF( IC.GE.97 .AND. IC.LE.122 ) THEN
SUBNAM( 1:1 ) = CHAR( IC-32 )
DO 10 I = 2, 6
$ SUBNAM( I:I ) = CHAR( IC-32 )
10 CONTINUE
END IF
-*
+C
ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
-*
-* EBCDIC character set
-*
+C
+C EBCDIC character set
+C
IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
$ ( IC.GE.145 .AND. IC.LE.153 ) .OR.
$ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
$ SUBNAM( I:I ) = CHAR( IC+64 )
20 CONTINUE
END IF
-*
+C
ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
-*
-* Prime machines: ASCII+128
-*
+C
+C Prime machines: ASCII+128
+C
IF( IC.GE.225 .AND. IC.LE.250 ) THEN
SUBNAM( 1:1 ) = CHAR( IC-32 )
DO 30 I = 2, 6
30 CONTINUE
END IF
END IF
-*
+C
C1 = SUBNAM( 1:1 )
SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
C2 = SUBNAM( 2:3 )
C3 = SUBNAM( 4:6 )
C4 = C3( 2:3 )
-*
+C
GO TO ( 110, 200, 300 ) ISPEC
-*
+C
110 CONTINUE
-*
-* ISPEC = 1: block size
-*
-* In these examples, separate code is provided for setting NB for
-* real and complex. We assume that NB will take the same value in
-* single or double precision.
-*
+C
+C ISPEC = 1: block size
+C
+C In these examples, separate code is provided for setting NB for
+C real and complex. We assume that NB will take the same value in
+C single or double precision.
+C
NB = 1
-*
+C
IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
END IF
ILAENV = NB
RETURN
-*
+C
200 CONTINUE
-*
-* ISPEC = 2: minimum block size
-*
+C
+C ISPEC = 2: minimum block size
+C
NBMIN = 2
IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
END IF
ILAENV = NBMIN
RETURN
-*
+C
300 CONTINUE
-*
-* ISPEC = 3: crossover point
-*
+C
+C ISPEC = 3: crossover point
+C
NX = 0
IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
END IF
ILAENV = NX
RETURN
-*
+C
400 CONTINUE
-*
-* ISPEC = 4: number of shifts (used by xHSEQR)
-*
+C
+C ISPEC = 4: number of shifts (used by xHSEQR)
+C
ILAENV = 6
RETURN
-*
+C
500 CONTINUE
-*
-* ISPEC = 5: minimum column dimension (not used)
-*
+C
+C ISPEC = 5: minimum column dimension (not used)
+C
ILAENV = 2
RETURN
-*
+C
600 CONTINUE
-*
-* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
-*
+C
+C ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
+C
ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
RETURN
-*
+C
700 CONTINUE
-*
-* ISPEC = 7: number of processors (not used)
-*
+C
+C ISPEC = 7: number of processors (not used)
+C
ILAENV = 1
RETURN
-*
+C
800 CONTINUE
-*
-* ISPEC = 8: crossover point for multishift (used by xHSEQR)
-*
+C
+C ISPEC = 8: crossover point for multishift (used by xHSEQR)
+C
ILAENV = 50
RETURN
-*
+C
900 CONTINUE
-*
-* ISPEC = 9: maximum size of the subproblems at the bottom of the
-* computation tree in the divide-and-conquer algorithm
-* (used by xGELSD and xGESDD)
-*
+C
+C ISPEC = 9: maximum size of the subproblems at the bottom of the
+C computation tree in the divide-and-conquer algorithm
+C (used by xGELSD and xGESDD)
+C
ILAENV = 25
RETURN
-*
+C
1000 CONTINUE
-*
-* ISPEC = 10: ieee NaN arithmetic can be trusted not to trap
-*
+C
+C ISPEC = 10: ieee NaN arithmetic can be trusted not to trap
+C
ILAENV = 1
IF (ILAENV .EQ. 1) THEN
ILAENV = IEEECK( 0, 0.0, 1.0 )
ENDIF
RETURN
-*
+C
1100 CONTINUE
-*
-* ISPEC = 11: infinity arithmetic can be trusted not to trap
-*
+C
+C ISPEC = 11: infinity arithmetic can be trusted not to trap
+C
ILAENV = 1
IF (ILAENV .EQ. 1) THEN
ILAENV = IEEECK( 1, 0.0, 1.0 )
ENDIF
RETURN
-*
-* End of ILAENV
-*
+C
+C End of ILAENV
+C
END
INTEGER FUNCTION IEEECK( ISPEC, ZERO, ONE )
-*
-* -- LAPACK auxiliary routine (version 3.2) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2006
-*
-* .. Scalar Arguments ..
+C
+C -- LAPACK auxiliary routine (version 3.2) --
+C Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+C November 2006
+C
+C .. Scalar Arguments ..
INTEGER ISPEC
REAL ZERO, ONE
-* ..
-*
-* Purpose
-* =======
-*
-* IEEECK is called from the ILAENV to verify that Inifinity and
-* possibly NaN arithmetic is safe (i.e. will not trap).
-*
-* Arguments
-* =========
-*
-* ISPEC (input) INTEGER
-* Specifies whether to test just for inifinity arithmetic
-* or whether to test for infinity and NaN arithmetic.
-* = 0: Verify infinity arithmetic only.
-* = 1: Verify infinity and NaN arithmetic.
-*
-* ZERO (input) REAL
-* Must contain the value 0.0
-* This is passed to prevent the compiler from optimizing
-* away this code.
-*
-* ONE (input) REAL
-* Must contain the value 1.0
-* This is passed to prevent the compiler from optimizing
-* away this code.
-*
-* RETURN VALUE: INTEGER
-* = 0: Arithmetic failed to produce the correct answers
-* = 1: Arithmetic produced the correct answers
-*
-* .. Local Scalars ..
+C ..
+C
+C Purpose
+C =======
+C
+C IEEECK is called from the ILAENV to verify that Inifinity and
+C possibly NaN arithmetic is safe (i.e. will not trap).
+C
+C Arguments
+C =========
+C
+C ISPEC (input) INTEGER
+C Specifies whether to test just for inifinity arithmetic
+C or whether to test for infinity and NaN arithmetic.
+C = 0: Verify infinity arithmetic only.
+C = 1: Verify infinity and NaN arithmetic.
+C
+C ZERO (input) REAL
+C Must contain the value 0.0
+C This is passed to prevent the compiler from optimizing
+C away this code.
+C
+C ONE (input) REAL
+C Must contain the value 1.0
+C This is passed to prevent the compiler from optimizing
+C away this code.
+C
+C RETURN VALUE: INTEGER
+C = 0: Arithmetic failed to produce the correct answers
+C = 1: Arithmetic produced the correct answers
+C
+C .. Local Scalars ..
REAL POSINF, NEGINF, NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGZRO,
$ NEWZRO
-* ..
-* .. Executable Statements ..
+C ..
+C .. Executable Statements ..
IEEECK = 1
POSINF = ONE /ZERO
-*
-* Return if we were only asked to check infinity arithmetic
-*
+C
+C Return if we were only asked to check infinity arithmetic
+C
IF (ISPEC .EQ. 0 ) RETURN
NAN1 = POSINF + NEGINF
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CGBRFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CGBSVXX uses the LU factorization to compute the solution to a
* user-provided factorizations and equilibration factors if they
* differ from what CGBSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by CGEBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array A. LDA >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
* JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values determined when A was
* balanced. The balanced A(i,j) = 0 if I > J and
* u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
* are left eigenvectors of (A,B).
*
-* Note: this routine performs "full balancing" on A and B -- see
-* "Further Details", below.
+* Note: this routine performs "full balancing" on A and B
*
* Arguments
* =========
* one of the forms lambda = alpha/beta or mu = beta/alpha.
* Since either lambda or mu may overflow, they should not,
* in general, be computed.
-
*
* VL (output) COMPLEX array, dimension (LDVL,N)
* If JOBVL = 'V', the left eigenvectors u(j) are stored
* .. Executable Statements ..
*
* Test input arguments
-* ====================
+* ====================
*
INFO = 0
LQUERY = ( LWORK.EQ.-1 )
NFXD = NFXD - 1
*
* Factorize fixed columns
-* =======================
+* =======================
*
* Compute the QR factorization of fixed columns and update
* remaining columns.
END IF
*
* Factorize free columns
-* ======================
+* ======================
*
IF( NFXD.LT.MINMN ) THEN
*
* On exit, the elements on and above the diagonal of the array
* contain the min(M,N)-by-N upper trapezoidal matrix R (R is
* upper triangular if M >= N); the elements below the diagonal
-* are the columns of V. See below for further details.
+* are the columns of V.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CGERFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CGESVXX uses the LU factorization to compute the solution to a
* user-provided factorizations and equilibration factors if they
* differ from what CGESVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by CGGBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array B. LDB >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* for which SELCTG is true.
*
* ALPHA (output) COMPLEX array, dimension (N)
+*
* BETA (output) COMPLEX array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
* for which SELCTG is true.
*
* ALPHA (output) COMPLEX array, dimension (N)
+*
* BETA (output) COMPLEX array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX array, dimension (N)
+*
* BETA (output) COMPLEX array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues.
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX array, dimension (N)
+*
* BETA (output) COMPLEX array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
* eigenvalues.
* if JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* For further explanation of the reciprocal condition numbers RCONDE
* and RCONDV, see section 4.11 of LAPACK User's Guide.
*
+* =====================================================================
+*
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
* LDQ >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found.
* = 'I': the IL-th through IU-th eigenvalues will be found.
-********** For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
-********** CSTEIN are called
+* For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
+* CSTEIN are called
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* indicating the nonzero elements in Z. The i-th eigenvector
* is nonzero only in elements ISUPPZ( 2*i-1 ) through
* ISUPPZ( 2*i ).
-********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
+* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
*
* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
$ ERR_BNDS_NORM( NRHS, * ),
$ ERR_BNDS_COMP( NRHS, * )
*
-* Purpose
+* Purpose
* =======
*
* CHERFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* CHESVXX uses the diagonal pivoting factorization to compute the
* solution to a complex system of linear equations A * X = B, where
* user-provided factorizations and equilibration factors if they
* differ from what CHESVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(R) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* The order of the matrices H, T, Q, and Z. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* ILO and IHI mark the rows and columns of H which are in
* Hessenberg form. It is assumed that A is already upper
* corresponding elements of A.
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* .. Array Arguments ..
COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
-* Purpose
+* Purpose
* =======
*
* CHSEQR computes the eigenvalues of a Hessenberg matrix H
* of a matrix A which has been reduced to the Hessenberg form H
* by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* JOB (input) CHARACTER*1
* = 'E': compute eigenvalues only;
* If INFO .GT. 0 and COMPZ = 'N', then Z is not
* accessed.
*
-* ================================================================
+* ================================================================
* Default values supplied by
* ILAENV(ISPEC,'CHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
* It is suggested that these defaults be adjusted in order
* for ISPEC=16 is 0. Otherwise the default for
* ISPEC=16 is 2.
*
-* ================================================================
+* ================================================================
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
-* A - COMPLEX array of DIMENSION ( LDA, n ).
+* A (input) COMPLEX array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
* where abs(Z) is the componentwise absolute value of the matrix
* or vector Z.
*
+* Arguments
+* =========
+*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
-* A - COMPLEX array of DIMENSION ( LDA, n ).
+* A (input) COMPLEX array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
COMPLEX X( * ), Y( * ), W( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CLA_WWADDW adds a vector W into a doubled-single vector (X, Y).
* This works for all extant IBM's hex and binary floating point
* arithmetics, but not for decimal.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* N (input) INTEGER
* The length of vectors X, Y, and W.
* The increment between successive values of CY. INCY <> 0.
*
* C (input) COMPLEX
+*
* S (input) COMPLEX
* C and S define the matrix
* [ C S ].
* =========
*
* X (input) COMPLEX
+*
* Y (input) COMPLEX
* The complex scalars X and Y.
*
* value THRESH (set below).
*
* CS1 (output) COMPLEX
+*
* SN1 (output) COMPLEX
* If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector
* for RT1.
* The eigenvalue of smaller absolute value.
*
* CS1 (output) REAL
+*
* SN1 (output) COMPLEX
* The vector (CS1, SN1) is a unit right eigenvector for RT1.
*
* = .FALSE.: the input matrices A and B are lower triangular.
*
* A1 (input) REAL
+*
* A2 (input) COMPLEX
+*
* A3 (input) REAL
* On entry, A1, A2 and A3 are elements of the input 2-by-2
* upper (lower) triangular matrix A.
*
* B1 (input) REAL
+*
* B2 (input) COMPLEX
+*
* B3 (input) REAL
* On entry, B1, B2 and B3 are elements of the input 2-by-2
* upper (lower) triangular matrix B.
*
* CSU (output) REAL
+*
* SNU (output) COMPLEX
* The desired unitary matrix U.
*
* CSV (output) REAL
+*
* SNV (output) COMPLEX
* The desired unitary matrix V.
*
* CSQ (output) REAL
+*
* SNQ (output) COMPLEX
* The desired unitary matrix Q.
*
COMPLEX H( LDH, * ), W( * ), Z( LDZ, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CLAHQR is an auxiliary routine called by CHSEQR to update the
* dealing with the Hessenberg submatrix in rows and columns ILO to
* IHI.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows and
* columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
* of the Schur form returned in H, with W(i) = H(i,i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* where U is the orthogonal matrix in (*)
* (regardless of the value of WANTT.)
*
-* Further Details
-* ===============
+* Further Details
+* ===============
*
* 02-96 Based on modifications by
* David Day, Sandia National Laboratory, USA
* (2) adopts the more conservative Ahues & Tisseur stopping
* criterion (LAWN 122, 1997).
*
-* =========================================================
+* =========================================================
*
* .. Parameters ..
INTEGER ITMAX
* The order of the matrix A. N >= 0. When N = 0, CLANHF is
* set to zero.
*
-* A (input) COMPLEX*16 array, dimension ( N*(N+1)/2 );
+* A (input) COMPLEX*16 array, dimension ( N*(N+1)/2 );
* On entry, the matrix A in RFP Format.
* RFP Format is described by TRANSR, UPLO and N as follows:
* If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* CLAQR0 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* of a matrix A which has been reduced to the Hessenberg form H
* by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
COMPLEX H( LDH, * ), V( * )
* ..
*
+* Purpose
+* =======
+*
* Given a 2-by-2 or 3-by-3 matrix H, CLAQR1 sets v to a
* scalar multiple of the first column of the product
*
* This is useful for starting double implicit shift bulges
* in the QR algorithm.
*
+* Arguments
+* =========
*
* N (input) integer
* Order of the matrix H. N must be either 2 or 3.
* the calling procedure. LDH.GE.N
*
* S1 (input) COMPLEX
+*
* S2 S1 and S2 are the shifts defining K in (*) above.
*
* V (output) COMPLEX array of dimension N
* A scalar multiple of the first column of the
* matrix K in (*).
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
*
* .. Parameters ..
COMPLEX ZERO
$ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
* ..
*
-* This subroutine is identical to CLAQR3 except that it avoids
-* recursion by calling CLAHQR instead of CLAQR4.
+* Purpose
+* =======
*
+* CLAQR2 is identical to CLAQR3 except that it avoids
+* recursion by calling CLAHQR instead of CLAQR4.
*
-* ******************************************************************
* Aggressive early deflation:
*
* This subroutine accepts as input an upper Hessenberg matrix
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
COMPLEX ZERO, ONE
PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ),
$ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
* ..
*
-* ******************************************************************
+* Purpose
+* =======
+*
* Aggressive early deflation:
*
-* This subroutine accepts as input an upper Hessenberg matrix
+* CLAQR3 accepts as input an upper Hessenberg matrix
* H and performs an unitary similarity transformation
* designed to detect and deflate fully converged eigenvalues from
* a trailing principal submatrix. On output H has been over-
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
COMPLEX ZERO, ONE
PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ),
COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
*
-* This subroutine implements one level of recursion for CLAQR0.
+*
+* Purpose
+* =======
+*
+* CLAQR4 implements one level of recursion for CLAQR0.
* It is a complete implementation of the small bulge multi-shift
* QR algorithm. It may be called by CLAQR0 and, for large enough
* deflation window size, it may be called by CLAQR3. This
* subroutine is identical to CLAQR0 except that it calls CLAQR2
* instead of CLAQR3.
*
-* Purpose
-* =======
-*
* CLAQR4 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* H = Z T Z**H, where T is an upper triangular matrix (the
* of a matrix A which has been reduced to the Hessenberg form H
* by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
$ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
* ..
*
-* This auxiliary subroutine called by CLAQR0 performs a
+* Purpose
+* =======
+*
+* CLAQR5 called by CLAQR0 performs a
* single small-bulge multi-shift QR sweep.
*
+* Arguments
+* =========
+*
* WANTT (input) logical scalar
* WANTT = .true. if the triangular Schur factor
* is being computed. WANTT is set to .false. otherwise.
* subroutine operates.
*
* KTOP (input) integer scalar
+*
* KBOT (input) integer scalar
* These are the first and last rows and columns of an
* isolated diagonal block upon which the QR sweep is to be
* calling procedure. LDH.GE.MAX(1,N).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
* LDWV is the leading dimension of WV as declared in the
* in the calling subroutine. LDWV.GE.NV.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* Reference:
*
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Level 3 Performance, SIAM Journal of Matrix Analysis,
* volume 23, pages 929--947, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
COMPLEX ZERO, ONE
PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ),
* The increment between elements of C. INCC > 0.
*
* Further Details
-* ======= =======
+* ===============
*
* 6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
*
* The order of the matrix. N >= 0.
*
* VL (input) REAL
+*
* VU (input) REAL
* Lower and upper bounds of the interval that contains the desired
* eigenvalues. VL < VU. Needed to compute gaps on the left or right
* end of the extremal eigenvalues in the desired RANGE.
*
-* D (input/output) REAL array, dimension (N)
+* D (input/output) REAL array, dimension (N)
* On entry, the N diagonal elements of the diagonal matrix D.
* On exit, D may be overwritten.
*
-* L (input/output) REAL array, dimension (N)
+* L (input/output) REAL array, dimension (N)
* On entry, the (N-1) subdiagonal elements of the unit
* bidiagonal matrix L are in elements 1 to N-1 of L
* (if the matrix is not splitted.) At the end of each block
* The total number of input eigenvalues. 0 <= M <= N.
*
* DOL (input) INTEGER
+*
* DOU (input) INTEGER
* If the user wants to compute only selected eigenvectors from all
* the eigenvalues supplied, he can specify an index range DOL:DOU.
* MINRGP (input) REAL
*
* RTOL1 (input) REAL
+*
* RTOL2 (input) REAL
* Parameters for bisection.
* An interval [LEFT,RIGHT] has converged if
* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*
-* W (input/output) REAL array, dimension (N)
+* W (input/output) REAL array, dimension (N)
* The first M elements of W contain the APPROXIMATE eigenvalues for
* which eigenvectors are to be computed. The eigenvalues
* should be grouped by split-off block and ordered from
* for their block. On exit, W holds the eigenvalues of the
* UNshifted matrix.
*
-* WERR (input/output) REAL array, dimension (N)
+* WERR (input/output) REAL array, dimension (N)
* The first M elements contain the semiwidth of the uncertainty
* interval of the corresponding eigenvalue in W
*
-* WGAP (input/output) REAL array, dimension (N)
+* WGAP (input/output) REAL array, dimension (N)
* The separation from the right neighbor eigenvalue in W.
*
* IBLOCK (input) INTEGER array, dimension (N)
* for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
* i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
*
-* GERS (input) REAL array, dimension (2*N)
+* GERS (input) REAL array, dimension (2*N)
* The N Gerschgorin intervals (the i-th Gerschgorin interval
* is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
* be computed from the original UNshifted matrix.
*
-* Z (output) COMPLEX array, dimension (LDZ, max(1,M) )
+* Z (output)COMPLEX array, dimension (LDZ, max(1,M) )
* If INFO = 0, the first M columns of Z contain the
* orthonormal eigenvectors of the matrix T
* corresponding to the input eigenvalues, with the i-th
* is nonzero only in elements ISUPPZ( 2*I-1 ) through
* ISUPPZ( 2*I ).
*
-* WORK (workspace) REAL array, dimension (12*N)
+* WORK (workspace) REAL array, dimension (12*N)
*
* IWORK (workspace) INTEGER array, dimension (7*N)
*
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, CTBSV
* is called, otherwise, specific code is used which checks for possible
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, CTPSV
* is called, otherwise, specific code is used which checks for possible
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, CTRSV
* is called, otherwise, specific code is used which checks for possible
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CPORFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T
* user-provided factorizations and equilibration factors if they
* differ from what CPOSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(S) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
SUBROUTINE CSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
*
* -- LAPACK PROTOTYPE routine (version 3.2.2) --
-*
* -- Written by Julie Langou of the Univ. of TN --
* May 2010
*
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
-*
+*
* WAY (input) CHARACTER*1
* = 'C': Convert
* = 'R': Revert
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CSYRFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* CSYSVXX uses the diagonal pivoting factorization to compute the
* user-provided factorizations and equilibration factors if they
* differ from what CSYSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(R) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* If WANTZ = .TRUE., LDZ >= N.
*
* IFST (input) INTEGER
+*
* ILST (input/output) INTEGER
* Specify the reordering of the diagonal blocks of (A, B).
* The block with row index IFST is moved to row ILST, by a
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DGBRFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* DGBSVXX uses the LU factorization to compute the solution to a
* double precision system of linear equations A * X = B, where A is an
* user-provided factorizations and equilibration factors if they
* differ from what DGBSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by DGEBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array A. LDA >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
* eigenvalue count as 2.)
*
* WR (output) DOUBLE PRECISION array, dimension (N)
+*
* WI (output) DOUBLE PRECISION array, dimension (N)
* WR and WI contain the real and imaginary parts,
* respectively, of the computed eigenvalues in the same order
* eigenvalue count as 2.)
*
* WR (output) DOUBLE PRECISION array, dimension (N)
+*
* WI (output) DOUBLE PRECISION array, dimension (N)
* WR and WI contain the real and imaginary parts, respectively,
* of the computed eigenvalues, in the same order that they
* The leading dimension of the array A. LDA >= max(1,N).
*
* WR (output) DOUBLE PRECISION array, dimension (N)
+*
* WI (output) DOUBLE PRECISION array, dimension (N)
* WR and WI contain the real and imaginary parts,
* respectively, of the computed eigenvalues. Complex
* The leading dimension of the array A. LDA >= max(1,N).
*
* WR (output) DOUBLE PRECISION array, dimension (N)
+*
* WI (output) DOUBLE PRECISION array, dimension (N)
* WR and WI contain the real and imaginary parts,
* respectively, of the computed eigenvalues. Complex
* JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values determined when A was
* balanced. The balanced A(i,j) = 0 if I > J and
*
* are left eigenvectors of (A,B).
*
-* Note: this routine performs "full balancing" on A and B -- see
-* "Further Details", below.
+* Note: this routine performs "full balancing" on A and B
*
* Arguments
* =========
* .. Executable Statements ..
*
* Test input arguments
-* ====================
+* ====================
*
INFO = 0
LQUERY = ( LWORK.EQ.-1 )
NFXD = NFXD - 1
*
* Factorize fixed columns
-* =======================
+* =======================
*
* Compute the QR factorization of fixed columns and update
* remaining columns.
END IF
*
* Factorize free columns
-* ======================
+* ======================
*
IF( NFXD.LT.MINMN ) THEN
*
* On exit, the elements on and above the diagonal of the array
* contain the min(M,N)-by-N upper trapezoidal matrix R (R is
* upper triangular if M >= N); the elements below the diagonal
-* are the columns of V. See below for further details.
+* are the columns of V.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
* The upper triangular block reflectors stored in compact form
* as a sequence of upper triangular blocks. See below
* for further details.
-*
+*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= NB.
*
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DGERFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* of SIGMA are the singular values of A. The columns of U and V are the
* left and the right singular vectors of A, respectively.
*
-* Further Details
-* ~~~~~~~~~~~~~~~
-* The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
-* rotations. The rotations are implemented as fast scaled rotations of
-* Anda and Park [1]. In the case of underflow of the Jacobi angle, a
-* modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
-* column interchanges of de Rijk [2]. The relative accuracy of the computed
-* singular values and the accuracy of the computed singular vectors (in
-* angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
-* The condition number that determines the accuracy in the full rank case
-* is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
-* spectral condition number. The best performance of this Jacobi SVD
-* procedure is achieved if used in an accelerated version of Drmac and
-* Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
-* Some tunning parameters (marked with [TP]) are available for the
-* implementer.
-* The computational range for the nonzero singular values is the machine
-* number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
-* denormalized singular values can be computed with the corresponding
-* gradual loss of accurate digits.
-*
-* Contributors
-* ~~~~~~~~~~~~
-* Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
-*
-* References
-* ~~~~~~~~~~
-* [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
-* SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
-* [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
-* singular value decomposition on a vector computer.
-* SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
-* [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
-* [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
-* value computation in floating point arithmetic.
-* SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
-* [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
-* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
-* LAPACK Working note 169.
-* [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
-* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
-* LAPACK Working note 170.
-* [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
-* QSVD, (H,K)-SVD computations.
-* Department of Mathematics, University of Zagreb, 2008.
-*
-* Bugs, Examples and Comments
-* ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-* Please report all bugs and send interesting test examples and comments to
-* drmac@math.hr. Thank you.
-*
* Arguments
* =========
*
* of sweeps. The output may still be useful. See the
* description of WORK.
*
+* Further Details
+* ===============
+*
+* The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
+* rotations. The rotations are implemented as fast scaled rotations of
+* Anda and Park [1]. In the case of underflow of the Jacobi angle, a
+* modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
+* column interchanges of de Rijk [2]. The relative accuracy of the computed
+* singular values and the accuracy of the computed singular vectors (in
+* angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
+* The condition number that determines the accuracy in the full rank case
+* is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
+* spectral condition number. The best performance of this Jacobi SVD
+* procedure is achieved if used in an accelerated version of Drmac and
+* Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
+* Some tunning parameters (marked with [TP]) are available for the
+* implementer.
+* The computational range for the nonzero singular values is the machine
+* number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
+* denormalized singular values can be computed with the corresponding
+* gradual loss of accurate digits.
+*
+* Contributors
+* ============
+*
+* Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*
+* References
+* ==========
+*
+* [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
+* SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
+* [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
+* singular value decomposition on a vector computer.
+* SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
+* [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
+* [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
+* value computation in floating point arithmetic.
+* SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
+* [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
+* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
+* LAPACK Working note 169.
+* [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
+* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
+* LAPACK Working note 170.
+* [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
+* QSVD, (H,K)-SVD computations.
+* Department of Mathematics, University of Zagreb, 2008.
+*
+* Bugs, Examples and Comments
+* ===========================
+* Please report all bugs and send interesting test examples and comments to
+* drmac@math.hr. Thank you.
+*
* =====================================================================
*
* .. Local Parameters ..
$ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK,
$ INFO )
*
-* -- LAPACK driver routine (version 3.2.2) --
+* -- LAPACK driver routine (version 3.2.2) --
* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
* -- Jason Riedy of Univ. of California Berkeley. --
* -- June 2010 --
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* DGESVXX uses the LU factorization to compute the solution to a
* double precision system of linear equations A * X = B, where A is an
* user-provided factorizations and equilibration factors if they
* differ from what DGESVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by DGGBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array B. LDB >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* SELCTG is true for either eigenvalue count as 2.)
*
* ALPHAR (output) DOUBLE PRECISION array, dimension (N)
+*
* ALPHAI (output) DOUBLE PRECISION array, dimension (N)
+*
* BETA (output) DOUBLE PRECISION array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i,
* SELCTG is true for either eigenvalue count as 2.)
*
* ALPHAR (output) DOUBLE PRECISION array, dimension (N)
+*
* ALPHAI (output) DOUBLE PRECISION array, dimension (N)
+*
* BETA (output) DOUBLE PRECISION array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHAR (output) DOUBLE PRECISION array, dimension (N)
+*
* ALPHAI (output) DOUBLE PRECISION array, dimension (N)
+*
* BETA (output) DOUBLE PRECISION array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. If ALPHAI(j) is zero, then
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHAR (output) DOUBLE PRECISION array, dimension (N)
+*
* ALPHAI (output) DOUBLE PRECISION array, dimension (N)
+*
* BETA (output) DOUBLE PRECISION array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. If ALPHAI(j) is zero, then
* if JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* The order of the matrices H, T, Q, and Z. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* ILO and IHI mark the rows and columns of H which are in
* Hessenberg form. It is assumed that A is already upper
DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
$ Z( LDZ, * )
* ..
-* Purpose
+* Purpose
* =======
*
* DHSEQR computes the eigenvalues of a Hessenberg matrix H
* of a matrix A which has been reduced to the Hessenberg form H
* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* JOB (input) CHARACTER*1
* = 'E': compute eigenvalues only;
* If INFO .GT. 0 and COMPZ = 'N', then Z is not
* accessed.
*
-* ================================================================
+* ================================================================
* Default values supplied by
* ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
* It is suggested that these defaults be adjusted in order
* for ISPEC=16 is 0. Otherwise the default for
* ISPEC=16 is 2.
*
-* ================================================================
+* ================================================================
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
* or vector Z.
*
* Arguments
-* ==========
+* =========
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* N must be at least zero.
* Unchanged on exit.
*
-* ALPHA - DOUBLE PRECISION .
+* ALPHA (input) DOUBLE PRECISION .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
-* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
+* A (input) DOUBLE PRECISION array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
* X. INCX must not be zero.
* Unchanged on exit.
*
-* BETA - DOUBLE PRECISION .
+* BETA (input) DOUBLE PRECISION .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
DOUBLE PRECISION X( * ), Y( * ), W( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DLA_WWADDW adds a vector W into a doubled-single vector (X, Y).
* This works for all extant IBM's hex and binary floating point
* arithmetics, but not for decimal.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* N (input) INTEGER
* The length of vectors X, Y, and W.
* ISAVE is used to save variables between calls to DLACN2
*
* Further Details
-* ======= =======
+* ===============
*
* Contributed by Nick Higham, University of Manchester.
* Originally named SONEST, dated March 16, 1988.
* =========
*
* A (input) DOUBLE PRECISION
+*
* B (input) DOUBLE PRECISION
+*
* C (input) DOUBLE PRECISION
+*
* D (input) DOUBLE PRECISION
* The scalars a, b, c, and d in the above expression.
*
* P (output) DOUBLE PRECISION
+*
* Q (output) DOUBLE PRECISION
* The scalars p and q in the above expression.
*
*
* PIVMIN (input) DOUBLE PRECISION
* The minimum absolute value of a "pivot" in the Sturm
-* sequence loop. This *must* be at least max |e(j)**2| *
-* safe_min and at least safe_min, where safe_min is at least
+* sequence loop.
+* This must be at least max |e(j)**2|*safe_min and at
+* least safe_min, where safe_min is at least
* the smallest number that can divide one without overflow.
*
* D (input) DOUBLE PRECISION array, dimension (N)
* DLAED4. K >= 0.
*
* KSTART (input) INTEGER
+*
* KSTOP (input) INTEGER
* The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
* are to be computed. 1 <= KSTART <= KSTOP <= K.
* The leading dimension of the array H. LDH >= max(1,N).
*
* WR (input) DOUBLE PRECISION
+*
* WI (input) DOUBLE PRECISION
* The real and imaginary parts of the eigenvalue of H whose
* corresponding right or left eigenvector is to be computed.
*
* VR (input/output) DOUBLE PRECISION array, dimension (N)
+*
* VI (input/output) DOUBLE PRECISION array, dimension (N)
* On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
* a real starting vector for inverse iteration using the real
* The eigenvalue of smaller absolute value.
*
* CS1 (output) DOUBLE PRECISION
+*
* SN1 (output) DOUBLE PRECISION
* The vector (CS1, SN1) is a unit right eigenvector for RT1.
*
* = .FALSE.: the input matrices A and B are lower triangular.
*
* A1 (input) DOUBLE PRECISION
+*
* A2 (input) DOUBLE PRECISION
+*
* A3 (input) DOUBLE PRECISION
* On entry, A1, A2 and A3 are elements of the input 2-by-2
* upper (lower) triangular matrix A.
*
* B1 (input) DOUBLE PRECISION
+*
* B2 (input) DOUBLE PRECISION
+*
* B3 (input) DOUBLE PRECISION
* On entry, B1, B2 and B3 are elements of the input 2-by-2
* upper (lower) triangular matrix B.
*
* CSU (output) DOUBLE PRECISION
+*
* SNU (output) DOUBLE PRECISION
* The desired orthogonal matrix U.
*
* CSV (output) DOUBLE PRECISION
+*
* SNV (output) DOUBLE PRECISION
* The desired orthogonal matrix V.
*
* CSQ (output) DOUBLE PRECISION
+*
* SNQ (output) DOUBLE PRECISION
* The desired orthogonal matrix Q.
*
* THe leading dimension of the array B. LDB >= 2.
*
* ALPHAR (output) DOUBLE PRECISION array, dimension (2)
+*
* ALPHAI (output) DOUBLE PRECISION array, dimension (2)
+*
* BETA (output) DOUBLE PRECISION array, dimension (2)
* (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
* pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may
DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DLAHQR is an auxiliary routine called by DHSEQR to update the
* dealing with the Hessenberg submatrix in rows and columns ILO to
* IHI.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper quasi-triangular in
* rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
* The leading dimension of the array H. LDH >= max(1,N).
*
* WR (output) DOUBLE PRECISION array, dimension (N)
+*
* WI (output) DOUBLE PRECISION array, dimension (N)
* The real and imaginary parts, respectively, of the computed
* eigenvalues ILO to IHI are stored in the corresponding
* WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* where U is the orthogonal matrix in (*)
* (regardless of the value of WANTT.)
*
-* Further Details
-* ===============
+* Further Details
+* ===============
*
* 02-96 Based on modifications by
* David Day, Sandia National Laboratory, USA
* (2) adopts the more conservative Ahues & Tisseur stopping
* criterion (LAWN 122, 1997).
*
-* =========================================================
+* =========================================================
*
* .. Parameters ..
INTEGER ITMAX
* =========
*
* N1 (input) INTEGER
+*
* N2 (input) INTEGER
* These arguements contain the respective lengths of the two
* sorted lists to be merged.
* for the final N2 elements.
*
* DTRD1 (input) INTEGER
+*
* DTRD2 (input) INTEGER
* These are the strides to be taken through the array A.
* Allowable strides are 1 and -1. They indicate whether a
* =========
*
* A (input/output) DOUBLE PRECISION
+*
* B (input/output) DOUBLE PRECISION
+*
* C (input/output) DOUBLE PRECISION
+*
* D (input/output) DOUBLE PRECISION
* On entry, the elements of the input matrix.
* On exit, they are overwritten by the elements of the
* standardised Schur form.
*
* RT1R (output) DOUBLE PRECISION
+*
* RT1I (output) DOUBLE PRECISION
+*
* RT2R (output) DOUBLE PRECISION
+*
* RT2I (output) DOUBLE PRECISION
* The real and imaginary parts of the eigenvalues. If the
* eigenvalues are a complex conjugate pair, RT1I > 0.
*
* CS (output) DOUBLE PRECISION
+*
* SN (output) DOUBLE PRECISION
* Parameters of the rotation matrix.
*
* =========
*
* X (input) DOUBLE PRECISION
+*
* Y (input) DOUBLE PRECISION
* X and Y specify the values x and y.
*
* =========
*
* X (input) DOUBLE PRECISION
+*
* Y (input) DOUBLE PRECISION
+*
* Z (input) DOUBLE PRECISION
* X, Y and Z specify the values x, y and z.
*
$ Z( LDZ, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* DLAQR0 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* of a matrix A which has been reduced to the Hessenberg form H
* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* The leading dimension of the array H. LDH .GE. max(1,N).
*
* WR (output) DOUBLE PRECISION array, dimension (IHI)
+*
* WI (output) DOUBLE PRECISION array, dimension (IHI)
* The real and imaginary parts, respectively, of the computed
* eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
* WI(i+1) = -WI(i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Arguments
+* =========
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
DOUBLE PRECISION H( LDH, * ), V( * )
* ..
*
+* Purpose
+* =======
+*
* Given a 2-by-2 or 3-by-3 matrix H, DLAQR1 sets v to a
* scalar multiple of the first column of the product
*
* This is useful for starting double implicit shift bulges
* in the QR algorithm.
*
+* Arguments
+* =========
*
* N (input) integer
* Order of the matrix H. N must be either 2 or 3.
* the calling procedure. LDH.GE.N
*
* SR1 (input) DOUBLE PRECISION
-* SI1 The shifts in (*).
-* SR2
-* SI2
+*
+* SI1 (input) DOUBLE PRECISION
+*
+* SR2 (input) DOUBLE PRECISION
+*
+* SI2 (input) DOUBLE PRECISION
+* The shifts in (*).
*
* V (output) DOUBLE PRECISION array of dimension N
* A scalar multiple of the first column of the
* matrix K in (*).
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO
$ Z( LDZ, * )
* ..
*
-* This subroutine is identical to DLAQR3 except that it avoids
-* recursion by calling DLAHQR instead of DLAQR4.
+* Purpose
+* =======
*
+* DLAQR2 is identical to DLAQR3 except that it avoids
+* recursion by calling DLAHQR instead of DLAQR4.
*
-* ******************************************************************
* Aggressive early deflation:
*
* This subroutine accepts as input an upper Hessenberg matrix
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the quasi-triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* subroutine.
*
* SR (output) DOUBLE PRECISION array, dimension (KBOT)
+*
* SI (output) DOUBLE PRECISION array, dimension (KBOT)
* On output, the real and imaginary parts of approximate
* eigenvalues that may be used for shifts are stored in
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
$ Z( LDZ, * )
* ..
*
-* ******************************************************************
+* Purpose
+* =======
+*
* Aggressive early deflation:
*
-* This subroutine accepts as input an upper Hessenberg matrix
+* DLAQR3 accepts as input an upper Hessenberg matrix
* H and performs an orthogonal similarity transformation
* designed to detect and deflate fully converged eigenvalues from
* a trailing principal submatrix. On output H has been over-
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the quasi-triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* subroutine.
*
* SR (output) DOUBLE PRECISION array, dimension (KBOT)
+*
* SI (output) DOUBLE PRECISION array, dimension (KBOT)
* On output, the real and imaginary parts of approximate
* eigenvalues that may be used for shifts are stored in
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
$ Z( LDZ, * )
* ..
*
-* This subroutine implements one level of recursion for DLAQR0.
+* Purpose
+* =======
+*
+* DLAQR4 implements one level of recursion for DLAQR0.
* It is a complete implementation of the small bulge multi-shift
* QR algorithm. It may be called by DLAQR0 and, for large enough
* deflation window size, it may be called by DLAQR3. This
* subroutine is identical to DLAQR0 except that it calls DLAQR2
* instead of DLAQR3.
*
-* Purpose
-* =======
-*
* DLAQR4 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* H = Z T Z**T, where T is an upper quasi-triangular matrix (the
* of a matrix A which has been reduced to the Hessenberg form H
* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* The leading dimension of the array H. LDH .GE. max(1,N).
*
* WR (output) DOUBLE PRECISION array, dimension (IHI)
+*
* WI (output) DOUBLE PRECISION array, dimension (IHI)
* The real and imaginary parts, respectively, of the computed
* eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
* WI(i+1) = -WI(i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
$ Z( LDZ, * )
* ..
*
-* This auxiliary subroutine called by DLAQR0 performs a
+* Purpose
+* =======
+*
+* DLAQR5, called by DLAQR0, performs a
* single small-bulge multi-shift QR sweep.
*
+* Arguments
+* =========
+*
* WANTT (input) logical scalar
* WANTT = .true. if the quasi-triangular Schur factor
* is being computed. WANTT is set to .false. otherwise.
* subroutine operates.
*
* KTOP (input) integer scalar
+*
* KBOT (input) integer scalar
* These are the first and last rows and columns of an
* isolated diagonal block upon which the QR sweep is to be
* must be positive and even.
*
* SR (input/output) DOUBLE PRECISION array of size (NSHFTS)
+*
* SI (input/output) DOUBLE PRECISION array of size (NSHFTS)
* SR contains the real parts and SI contains the imaginary
* parts of the NSHFTS shifts of origin that define the
* calling procedure. LDH.GE.MAX(1,N).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
* LDWV is the leading dimension of WV as declared in the
* in the calling subroutine. LDWV.GE.NV.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* Reference:
*
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Level 3 Performance, SIAM Journal of Matrix Analysis,
* volume 23, pages 929--947, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
* etc., and the NSPLIT-th consists of rows/columns
* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
*
-*
* INFO (output) INTEGER
* = 0: successful exit
*
* The index of the last eigenvalue to be computed.
*
* RTOL1 (input) DOUBLE PRECISION
+*
* RTOL2 (input) DOUBLE PRECISION
* Tolerance for the convergence of the bisection intervals.
* An interval [LEFT,RIGHT] has converged if
* The order of the matrix. N > 0.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* The lower and upper bounds for the eigenvalues.
*
* that are in the interval (VL,VU]
*
* LCNT (output) INTEGER
+*
* RCNT (output) INTEGER
* The left and right negcounts of the interval.
*
* The order of the tridiagonal matrix T. N >= 0.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. Eigenvalues less than or equal
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* in W.
*
* WL (output) DOUBLE PRECISION
+*
* WU (output) DOUBLE PRECISION
* The interval (WL, WU] contains all the wanted eigenvalues.
* If RANGE='V', then WL=VL and WU=VU.
* The order of the matrix. N > 0.
*
* VL (input/output) DOUBLE PRECISION
+*
* VU (input/output) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds for the eigenvalues.
* Eigenvalues less than or equal to VL, or greater than VU,
* part of the spectrum.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* 1 <= I <= NSPLIT, have been set to zero
*
* RTOL1 (input) DOUBLE PRECISION
+*
* RTOL2 (input) DOUBLE PRECISION
* Parameters for bisection.
* An interval [LEFT,RIGHT] has converged if
* The index of the eigenvalues to be returned.
*
* GL (input) DOUBLE PRECISION
+*
* GU (input) DOUBLE PRECISION
* An upper and a lower bound on the eigenvalue.
*
* The order of the matrix. N >= 0.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* Lower and upper bounds of the interval that contains the desired
* eigenvalues. VL < VU. Needed to compute gaps on the left or right
* The total number of input eigenvalues. 0 <= M <= N.
*
* DOL (input) INTEGER
+*
* DOU (input) INTEGER
* If the user wants to compute only selected eigenvectors from all
* the eigenvalues supplied, he can specify an index range DOL:DOU.
* MINRGP (input) DOUBLE PRECISION
*
* RTOL1 (input) DOUBLE PRECISION
+*
* RTOL2 (input) DOUBLE PRECISION
* Parameters for bisection.
* An interval [LEFT,RIGHT] has converged if
* 'Q' or 'Z'.
*
* CFROM (input) DOUBLE PRECISION
+*
* CTO (input) DOUBLE PRECISION
* The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
* without over/underflow if the final result CTO*A(I,J)/CFROM
*
* Further Details
* ===============
+*
* Local Variables: I0:N0 defines a current unreduced segment of Z.
* The shifts are accumulated in SIGMA. Iteration count is in ITER.
* Ping-pong is controlled by PP (alternates between 0 and 1).
*
* Further Details
* ===============
+*
* CNST1 = 9/16
*
* =====================================================================
* abs(SSMAX) is the larger singular value.
*
* SNL (output) DOUBLE PRECISION
+*
* CSL (output) DOUBLE PRECISION
* The vector (CSL, SNL) is a unit left singular vector for the
* singular value abs(SSMAX).
*
* SNR (output) DOUBLE PRECISION
+*
* CSR (output) DOUBLE PRECISION
* The vector (CSR, SNR) is a unit right singular vector for the
* singular value abs(SSMAX).
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, DTBSV
* is called, otherwise, specific code is used which checks for possible
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, DTPSV
* is called, otherwise, specific code is used which checks for possible
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, DTRSV
* is called, otherwise, specific code is used which checks for possible
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DPORFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T
* user-provided factorizations and equilibration factors if they
* differ from what DPOSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(S) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* LDQ >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* corresponding elements of A.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found.
* = 'I': the IL-th through IU-th eigenvalues will be found.
-********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
-********** DSTEIN are called
+* For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
+* DSTEIN are called
*
* N (input) INTEGER
* The order of the matrix. N >= 0.
* to avoid over/underflow in computing the eigenvalues.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* indicating the nonzero elements in Z. The i-th eigenvector
* is nonzero only in elements ISUPPZ( 2*i-1 ) through
* ISUPPZ( 2*i ).
-********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
+* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
*
* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal (and
* to avoid over/underflow in computing the eigenvalues.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
-*
+*
* WAY (input) CHARACTER*1
* = 'C': Convert
* = 'R': Revert
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found.
* = 'I': the IL-th through IU-th eigenvalues will be found.
-********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
-********** DSTEIN are called
+* For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
+* DSTEIN are called
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* indicating the nonzero elements in Z. The i-th eigenvector
* is nonzero only in elements ISUPPZ( 2*i-1 ) through
* ISUPPZ( 2*i ).
-********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
+* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
*
* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* The leading dimension of the array B. LDB >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DSYRFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* DSYSVXX uses the diagonal pivoting factorization to compute the
* user-provided factorizations and equilibration factors if they
* differ from what DSYSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(R) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* If WANTZ = .TRUE., LDZ >= N.
*
* IFST (input/output) INTEGER
+*
* ILST (input/output) INTEGER
* Specify the reordering of the diagonal blocks of (A, B).
* The block with row index IFST is moved to row ILST, by a
* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
*
* Where
-* K-1 L-1
+* K-1 T L-1
* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
* I=1 J=1
*
* N (input) INTEGER
* The number of columns of the matrix A.
*
-* A (input) COMPLEX array, dimension (LDA,N)
+* A (input)COMPLEX array, dimension (LDA,N)
* The m by n matrix A.
*
* LDA (input) INTEGER
*
* Arguments
* =========
+*
* DIAG (input) CHARACTER*1
* = 'N': A is non-unit triangular;
* = 'U': A is unit triangular.
* be specified as OPTS = 'UTN'.
*
* N1 (input) INTEGER
+*
* N2 (input) INTEGER
+*
* N3 (input) INTEGER
+*
* N4 (input) INTEGER
* Problem dimensions for the subroutine NAME; these may not all
* be required.
*
* Arguments
* =========
+*
* PREC (input) CHARACTER
* Specifies the form of the system of equations:
* = 'S': Single
*
* Arguments
* =========
+*
* TRANS (input) CHARACTER*1
* Specifies the form of the system of equations:
* = 'N': No transpose
*
* Arguments
* =========
+*
* UPLO (input) CHARACTER
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* Arguments
* =========
+*
* VERS_MAJOR (output) INTEGER
* return the lapack major version
+*
* VERS_MINOR (output) INTEGER
* return the lapack minor version from the major version
+*
* VERS_PATCH (output) INTEGER
* return the lapack patch version from the minor version
+*
* =====================================================================
*
INTEGER VERS_MAJOR, VERS_MINOR, VERS_PATCH
* N is the order of the Hessenberg matrix H.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular
* in rows and columns 1:ILO-1 and IHI+1:N.
* (See ISPEC=16 above for details.)
* Default: 3.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
INTEGER INMIN, INWIN, INIBL, ISHFTS, IACC22
PARAMETER ( INMIN = 12, INWIN = 13, INIBL = 14,
* The number of characters in CA and CB to be compared.
*
* CA (input) CHARACTER*(*)
+*
* CB (input) CHARACTER*(*)
* CA and CB specify two character strings of length at least N.
* Only the first N characters of each string will be accessed.
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SGBRFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SGBSVXX uses the LU factorization to compute the solution to a
* user-provided factorizations and equilibration factors if they
* differ from what SGBSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by SGEBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array A. LDA >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
* eigenvalue count as 2.)
*
* WR (output) REAL array, dimension (N)
+*
* WI (output) REAL array, dimension (N)
* WR and WI contain the real and imaginary parts,
* respectively, of the computed eigenvalues in the same order
* eigenvalue count as 2.)
*
* WR (output) REAL array, dimension (N)
+*
* WI (output) REAL array, dimension (N)
* WR and WI contain the real and imaginary parts, respectively,
* of the computed eigenvalues, in the same order that they
* The leading dimension of the array A. LDA >= max(1,N).
*
* WR (output) REAL array, dimension (N)
+*
* WI (output) REAL array, dimension (N)
* WR and WI contain the real and imaginary parts,
* respectively, of the computed eigenvalues. Complex
* The leading dimension of the array A. LDA >= max(1,N).
*
* WR (output) REAL array, dimension (N)
+*
* WI (output) REAL array, dimension (N)
* WR and WI contain the real and imaginary parts,
* respectively, of the computed eigenvalues. Complex
* JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values determined when A was
* balanced. The balanced A(i,j) = 0 if I > J and
*
* are left eigenvectors of (A,B).
*
-* Note: this routine performs "full balancing" on A and B -- see
-* "Further Details", below.
+* Note: this routine performs "full balancing" on A and B
*
* Arguments
* =========
NFXD = NFXD - 1
*
* Factorize fixed columns
-* =======================
+* =======================
*
* Compute the QR factorization of fixed columns and update
* remaining columns.
END IF
*
* Factorize free columns
-* ======================
+* ======================
*
IF( NFXD.LT.MINMN ) THEN
*
* On exit, the elements on and above the diagonal of the array
* contain the min(M,N)-by-N upper trapezoidal matrix R (R is
* upper triangular if M >= N); the elements below the diagonal
-* are the columns of V. See below for further details.
+* are the columns of V.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
* The upper triangular block reflectors stored in compact form
* as a sequence of upper triangular blocks. See below
* for further details.
-*
+*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= NB.
*
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SGERFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SGESVXX uses the LU factorization to compute the solution to a
* user-provided factorizations and equilibration factors if they
* differ from what SGESVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by SGGBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array B. LDB >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* SELCTG is true for either eigenvalue count as 2.)
*
* ALPHAR (output) REAL array, dimension (N)
+*
* ALPHAI (output) REAL array, dimension (N)
+*
* BETA (output) REAL array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i,
* SELCTG is true for either eigenvalue count as 2.)
*
* ALPHAR (output) REAL array, dimension (N)
+*
* ALPHAI (output) REAL array, dimension (N)
+*
* BETA (output) REAL array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHAR (output) REAL array, dimension (N)
+*
* ALPHAI (output) REAL array, dimension (N)
+*
* BETA (output) REAL array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. If ALPHAI(j) is zero, then
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHAR (output) REAL array, dimension (N)
+*
* ALPHAI (output) REAL array, dimension (N)
+*
* BETA (output) REAL array, dimension (N)
* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
* be the generalized eigenvalues. If ALPHAI(j) is zero, then
* if JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* The order of the matrices H, T, Q, and Z. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* ILO and IHI mark the rows and columns of H which are in
* Hessenberg form. It is assumed that A is already upper
REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ),
$ Z( LDZ, * )
* ..
-* Purpose
+* Purpose
* =======
*
* SHSEQR computes the eigenvalues of a Hessenberg matrix H
* of a matrix A which has been reduced to the Hessenberg form H
* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* JOB (input) CHARACTER*1
* = 'E': compute eigenvalues only;
* If INFO .GT. 0 and COMPZ = 'N', then Z is not
* accessed.
*
-* ================================================================
+* ================================================================
* Default values supplied by
* ILAENV(ISPEC,'SHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
* It is suggested that these defaults be adjusted in order
* for ISPEC=16 is 0. Otherwise the default for
* ISPEC=16 is 2.
*
-* ================================================================
+* ================================================================
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
* or vector Z.
*
* Arguments
-* ==========
+* =========
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
-* A - REAL array of DIMENSION ( LDA, n ).
+* A (input) REAL array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
REAL X( * ), Y( * ), W( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SLA_WWADDW adds a vector W into a doubled-single vector (X, Y).
* This works for all extant IBM's hex and binary floating point
* arithmetics, but not for decimal.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* N (input) INTEGER
* The length of vectors X, Y, and W.
* ISAVE is used to save variables between calls to SLACN2
*
* Further Details
-* ======= =======
+* ===============
*
* Contributed by Nick Higham, University of Manchester.
* Originally named SONEST, dated March 16, 1988.
* =========
*
* A (input) REAL
+*
* B (input) REAL
+*
* C (input) REAL
+*
* D (input) REAL
* The scalars a, b, c, and d in the above expression.
*
* P (output) REAL
+*
* Q (output) REAL
* The scalars p and q in the above expression.
*
*
* PIVMIN (input) REAL
* The minimum absolute value of a "pivot" in the Sturm
-* sequence loop. This *must* be at least max |e(j)**2| *
-* safe_min and at least safe_min, where safe_min is at least
+* sequence loop.
+* This must be at least max |e(j)**2|*safe_min and at
+* least safe_min, where safe_min is at least
* the smallest number that can divide one without overflow.
*
* D (input) REAL array, dimension (N)
* SLAED4. K >= 0.
*
* KSTART (input) INTEGER
+*
* KSTOP (input) INTEGER
* The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
* are to be computed. 1 <= KSTART <= KSTOP <= K.
* The leading dimension of the array H. LDH >= max(1,N).
*
* WR (input) REAL
+*
* WI (input) REAL
* The real and imaginary parts of the eigenvalue of H whose
* corresponding right or left eigenvector is to be computed.
*
* VR (input/output) REAL array, dimension (N)
+*
* VI (input/output) REAL array, dimension (N)
* On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
* a real starting vector for inverse iteration using the real
* The eigenvalue of smaller absolute value.
*
* CS1 (output) REAL
+*
* SN1 (output) REAL
* The vector (CS1, SN1) is a unit right eigenvector for RT1.
*
* = .FALSE.: the input matrices A and B are lower triangular.
*
* A1 (input) REAL
+*
* A2 (input) REAL
+*
* A3 (input) REAL
* On entry, A1, A2 and A3 are elements of the input 2-by-2
* upper (lower) triangular matrix A.
*
* B1 (input) REAL
+*
* B2 (input) REAL
+*
* B3 (input) REAL
* On entry, B1, B2 and B3 are elements of the input 2-by-2
* upper (lower) triangular matrix B.
*
* CSU (output) REAL
+*
* SNU (output) REAL
* The desired orthogonal matrix U.
*
* CSV (output) REAL
+*
* SNV (output) REAL
* The desired orthogonal matrix V.
*
* CSQ (output) REAL
+*
* SNQ (output) REAL
* The desired orthogonal matrix Q.
*
* THe leading dimension of the array B. LDB >= 2.
*
* ALPHAR (output) REAL array, dimension (2)
+*
* ALPHAI (output) REAL array, dimension (2)
+*
* BETA (output) REAL array, dimension (2)
* (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
* pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may
REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SLAHQR is an auxiliary routine called by SHSEQR to update the
* dealing with the Hessenberg submatrix in rows and columns ILO to
* IHI.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper quasi-triangular in
* rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
* The leading dimension of the array H. LDH >= max(1,N).
*
* WR (output) REAL array, dimension (N)
+*
* WI (output) REAL array, dimension (N)
* The real and imaginary parts, respectively, of the computed
* eigenvalues ILO to IHI are stored in the corresponding
* WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* where U is the orthogonal matrix in (*)
* (regardless of the value of WANTT.)
*
-* Further Details
-* ===============
+* Further Details
+* ===============
*
* 02-96 Based on modifications by
* David Day, Sandia National Laboratory, USA
* (2) adopts the more conservative Ahues & Tisseur stopping
* criterion (LAWN 122, 1997).
*
-* =========================================================
+* =========================================================
*
* .. Parameters ..
INTEGER ITMAX
* =========
*
* N1 (input) INTEGER
+*
* N2 (input) INTEGER
* These arguements contain the respective lengths of the two
* sorted lists to be merged.
* for the final N2 elements.
*
* STRD1 (input) INTEGER
+*
* STRD2 (input) INTEGER
* These are the strides to be taken through the array A.
* Allowable strides are 1 and -1. They indicate whether a
* N (input) INTEGER
* The order of the matrix.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* The N diagonal elements of the diagonal matrix D.
*
-* LLD (input) REAL array, dimension (N-1)
+* LLD (input) REAL array, dimension (N-1)
* The (N-1) elements L(i)*L(i)*D(i).
*
* SIGMA (input) REAL
* =========
*
* A (input/output) REAL
+*
* B (input/output) REAL
+*
* C (input/output) REAL
+*
* D (input/output) REAL
* On entry, the elements of the input matrix.
* On exit, they are overwritten by the elements of the
* standardised Schur form.
*
* RT1R (output) REAL
+*
* RT1I (output) REAL
+*
* RT2R (output) REAL
+*
* RT2I (output) REAL
* The real and imaginary parts of the eigenvalues. If the
* eigenvalues are a complex conjugate pair, RT1I > 0.
*
* CS (output) REAL
+*
* SN (output) REAL
* Parameters of the rotation matrix.
*
* =========
*
* X (input) REAL
+*
* Y (input) REAL
* X and Y specify the values x and y.
*
* =========
*
* X (input) REAL
+*
* Y (input) REAL
+*
* Z (input) REAL
* X, Y and Z specify the values x, y and z.
*
$ Z( LDZ, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* SLAQR0 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* of a matrix A which has been reduced to the Hessenberg form H
* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* The leading dimension of the array H. LDH .GE. max(1,N).
*
* WR (output) REAL array, dimension (IHI)
+*
* WI (output) REAL array, dimension (IHI)
* The real and imaginary parts, respectively, of the computed
* eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
* WI(i+1) = -WI(i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
REAL H( LDH, * ), V( * )
* ..
*
+* Purpose
+* =======
+*
* Given a 2-by-2 or 3-by-3 matrix H, SLAQR1 sets v to a
* scalar multiple of the first column of the product
*
* This is useful for starting double implicit shift bulges
* in the QR algorithm.
*
+* Arguments
+* =========
*
* N (input) integer
* Order of the matrix H. N must be either 2 or 3.
* the calling procedure. LDH.GE.N
*
* SR1 (input) REAL
-* SI1 The shifts in (*).
-* SR2
-* SI2
+*
+* SI1 (input) REAL
+*
+* SR2 (input) REAL
+*
+* SI2 (input) REAL
+* The shifts in (*).
*
* V (output) REAL array of dimension N
* A scalar multiple of the first column of the
* matrix K in (*).
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
*
* .. Parameters ..
REAL ZERO
$ Z( LDZ, * )
* ..
*
-* This subroutine is identical to SLAQR3 except that it avoids
-* recursion by calling SLAHQR instead of SLAQR4.
+* Purpose
+* =======
*
+* SLAQR2 is identical to SLAQR3 except that it avoids
+* recursion by calling SLAHQR instead of SLAQR4.
*
-* ******************************************************************
* Aggressive early deflation:
*
* This subroutine accepts as input an upper Hessenberg matrix
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the quasi-triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* subroutine.
*
* SR (output) REAL array, dimension KBOT
+*
* SI (output) REAL array, dimension KBOT
* On output, the real and imaginary parts of approximate
* eigenvalues that may be used for shifts are stored in
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 )
$ Z( LDZ, * )
* ..
*
-* ******************************************************************
+* Purpose
+* =======
+*
* Aggressive early deflation:
*
-* This subroutine accepts as input an upper Hessenberg matrix
+* SLAQR3 accepts as input an upper Hessenberg matrix
* H and performs an orthogonal similarity transformation
* designed to detect and deflate fully converged eigenvalues from
* a trailing principal submatrix. On output H has been over-
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the quasi-triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* subroutine.
*
* SR (output) REAL array, dimension KBOT
+*
* SI (output) REAL array, dimension KBOT
* On output, the real and imaginary parts of approximate
* eigenvalues that may be used for shifts are stored in
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 )
$ Z( LDZ, * )
* ..
*
-* This subroutine implements one level of recursion for SLAQR0.
+* Purpose
+* =======
+*
+* SLAQR4 implements one level of recursion for SLAQR0.
* It is a complete implementation of the small bulge multi-shift
* QR algorithm. It may be called by SLAQR0 and, for large enough
* deflation window size, it may be called by SLAQR3. This
* subroutine is identical to SLAQR0 except that it calls SLAQR2
* instead of SLAQR3.
*
-* Purpose
-* =======
-*
* SLAQR4 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* H = Z T Z**T, where T is an upper quasi-triangular matrix (the
* of a matrix A which has been reduced to the Hessenberg form H
* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* The leading dimension of the array H. LDH .GE. max(1,N).
*
* WR (output) REAL array, dimension (IHI)
+*
* WI (output) REAL array, dimension (IHI)
* The real and imaginary parts, respectively, of the computed
* eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
* WI(i+1) = -WI(i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
$ Z( LDZ, * )
* ..
*
-* This auxiliary subroutine called by SLAQR0 performs a
+* Purpose
+* =======
+*
+* SLAQR5, called by SLAQR0, performs a
* single small-bulge multi-shift QR sweep.
*
+* Arguments
+* =========
+*
* WANTT (input) logical scalar
* WANTT = .true. if the quasi-triangular Schur factor
* is being computed. WANTT is set to .false. otherwise.
* subroutine operates.
*
* KTOP (input) integer scalar
+*
* KBOT (input) integer scalar
* These are the first and last rows and columns of an
* isolated diagonal block upon which the QR sweep is to be
* must be positive and even.
*
* SR (input/output) REAL array of size (NSHFTS)
+*
* SI (input/output) REAL array of size (NSHFTS)
* SR contains the real parts and SI contains the imaginary
* parts of the NSHFTS shifts of origin that define the
* calling procedure. LDH.GE.MAX(1,N).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
* LDWV is the leading dimension of WV as declared in the
* in the calling subroutine. LDWV.GE.NV.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* Reference:
*
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Level 3 Performance, SIAM Journal of Matrix Analysis,
* volume 23, pages 929--947, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 )
* N (input) INTEGER
* The order of the matrix. N > 0.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* On entry, the N diagonal elements of the tridiagonal
* matrix T.
*
-* E (input/output) REAL array, dimension (N)
+* E (input/output) REAL array, dimension (N)
* On entry, the first (N-1) entries contain the subdiagonal
* elements of the tridiagonal matrix T; E(N) need not be set.
* On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT,
* are set to zero, the other entries of E are untouched.
*
-* E2 (input/output) REAL array, dimension (N)
+* E2 (input/output) REAL array, dimension (N)
* On entry, the first (N-1) entries contain the SQUARES of the
* subdiagonal elements of the tridiagonal matrix T;
* E2(N) need not be set.
* N (input) INTEGER
* The order of the matrix.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* The N diagonal elements of the diagonal matrix D.
*
-* LLD (input) REAL array, dimension (N-1)
+* LLD (input) REAL array, dimension (N-1)
* The (N-1) elements L(i)*L(i)*D(i).
*
* IFIRST (input) INTEGER
* The index of the last eigenvalue to be computed.
*
* RTOL1 (input) REAL
+*
* RTOL2 (input) REAL
* Tolerance for the convergence of the bisection intervals.
* An interval [LEFT,RIGHT] has converged if
* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
* through ILAST-OFFSET elements of these arrays are to be used.
*
-* W (input/output) REAL array, dimension (N)
+* W (input/output) REAL array, dimension (N)
* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
* estimates of the eigenvalues of L D L^T indexed IFIRST throug
* ILAST.
* On output, these estimates are refined.
*
-* WGAP (input/output) REAL array, dimension (N-1)
+* WGAP (input/output) REAL array, dimension (N-1)
* On input, the (estimated) gaps between consecutive
* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
* then WGAP(IFIRST-OFFSET) must be set to ZERO.
* On output, these gaps are refined.
*
-* WERR (input/output) REAL array, dimension (N)
+* WERR (input/output) REAL array, dimension (N)
* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
* the errors in the estimates of the corresponding elements in W.
* On output, these errors are refined.
*
-* WORK (workspace) REAL array, dimension (2*N)
+* WORK (workspace) REAL array, dimension (2*N)
* Workspace.
*
* IWORK (workspace) INTEGER array, dimension (2*N)
* The order of the matrix. N > 0.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* The lower and upper bounds for the eigenvalues.
*
* that are in the interval (VL,VU]
*
* LCNT (output) INTEGER
+*
* RCNT (output) INTEGER
* The left and right negcounts of the interval.
*
* The order of the tridiagonal matrix T. N >= 0.
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. Eigenvalues less than or equal
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
* Not referenced if RANGE = 'A' or 'V'.
*
-* GERS (input) REAL array, dimension (2*N)
+* GERS (input) REAL array, dimension (2*N)
* The N Gerschgorin intervals (the i-th Gerschgorin interval
* is (GERS(2*i-1), GERS(2*i)).
*
* sufficiently small, i.e., converged. Note: this should
* always be at least radix*machine epsilon.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* The n diagonal elements of the tridiagonal matrix T.
*
-* E (input) REAL array, dimension (N-1)
+* E (input) REAL array, dimension (N-1)
* The (n-1) off-diagonal elements of the tridiagonal matrix T.
*
-* E2 (input) REAL array, dimension (N-1)
+* E2 (input) REAL array, dimension (N-1)
* The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
*
* PIVMIN (input) REAL
* The actual number of eigenvalues found. 0 <= M <= N.
* (See also the description of INFO=2,3.)
*
-* W (output) REAL array, dimension (N)
+* W (output) REAL array, dimension (N)
* On exit, the first M elements of W will contain the
* eigenvalue approximations. SLARRD computes an interval
* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
* WERR(j) = abs( a_j - b_j)/2
*
-* WERR (output) REAL array, dimension (N)
+* WERR (output) REAL array, dimension (N)
* The error bound on the corresponding eigenvalue approximation
* in W.
*
* WL (output) REAL
+*
* WU (output) REAL
* The interval (WL, WU] contains all the wanted eigenvalues.
* If RANGE='V', then WL=VL and WU=VU.
* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
* i-th eigenvalue W(i) is the j-th eigenvalue in block k.
*
-* WORK (workspace) REAL array, dimension (4*N)
+* WORK (workspace) REAL array, dimension (4*N)
*
* IWORK (workspace) INTEGER array, dimension (3*N)
*
* The order of the matrix. N > 0.
*
* VL (input/output) REAL
+*
* VU (input/output) REAL
* If RANGE='V', the lower and upper bounds for the eigenvalues.
* Eigenvalues less than or equal to VL, or greater than VU,
* part of the spectrum.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* 1 <= IL <= IU <= N.
*
-* D (input/output) REAL array, dimension (N)
+* D (input/output) REAL array, dimension (N)
* On entry, the N diagonal elements of the tridiagonal
* matrix T.
* On exit, the N diagonal elements of the diagonal
* matrices D_i.
*
-* E (input/output) REAL array, dimension (N)
+* E (input/output) REAL array, dimension (N)
* On entry, the first (N-1) entries contain the subdiagonal
* elements of the tridiagonal matrix T; E(N) need not be set.
* On exit, E contains the subdiagonal elements of the unit
* bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
* 1 <= I <= NSPLIT, contain the base points sigma_i on output.
*
-* E2 (input/output) REAL array, dimension (N)
+* E2 (input/output) REAL array, dimension (N)
* On entry, the first (N-1) entries contain the SQUARES of the
* subdiagonal elements of the tridiagonal matrix T;
* E2(N) need not be set.
* 1 <= I <= NSPLIT, have been set to zero
*
* RTOL1 (input) REAL
+*
* RTOL2 (input) REAL
* Parameters for bisection.
* An interval [LEFT,RIGHT] has converged if
* The total number of eigenvalues (of all L_i D_i L_i^T)
* found.
*
-* W (output) REAL array, dimension (N)
+* W (output) REAL array, dimension (N)
* The first M elements contain the eigenvalues. The
* eigenvalues of each of the blocks, L_i D_i L_i^T, are
* sorted in ascending order ( SLARRE may use the
* remaining N-M elements as workspace).
*
-* WERR (output) REAL array, dimension (N)
+* WERR (output) REAL array, dimension (N)
* The error bound on the corresponding eigenvalue in W.
*
-* WGAP (output) REAL array, dimension (N)
+* WGAP (output) REAL array, dimension (N)
* The separation from the right neighbor eigenvalue in W.
* The gap is only with respect to the eigenvalues of the same block
* as each block has its own representation tree.
* for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
* i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
*
-* GERS (output) REAL array, dimension (2*N)
+* GERS (output) REAL array, dimension (2*N)
* The N Gerschgorin intervals (the i-th Gerschgorin interval
* is (GERS(2*i-1), GERS(2*i)).
*
* PIVMIN (output) REAL
* The minimum pivot in the Sturm sequence for T.
*
-* WORK (workspace) REAL array, dimension (6*N)
+* WORK (workspace) REAL array, dimension (6*N)
* Workspace.
*
* IWORK (workspace) INTEGER array, dimension (5*N)
* =-6: Problem in SLASQ2.
*
* Further Details
+* ===============
+*
* The base representations are required to suffer very little
* element growth and consequently define all their eigenvalues to
* high relative accuracy.
-* ===============
*
* Based on contributions by
* Beresford Parlett, University of California, Berkeley, USA
* SIGMA (output) REAL
* The shift used to form L(+) D(+) L(+)^T.
*
-* DPLUS (output) REAL array, dimension (N)
+* DPLUS (output) REAL array, dimension (N)
* The N diagonal elements of the diagonal matrix D(+).
*
-* LPLUS (output) REAL array, dimension (N-1)
+* LPLUS (output) REAL array, dimension (N-1)
* The first (N-1) elements of LPLUS contain the subdiagonal
* elements of the unit bidiagonal matrix L(+).
*
* N (input) INTEGER
* The order of the matrix.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* The N diagonal elements of T.
*
-* E2 (input) REAL array, dimension (N-1)
+* E2 (input) REAL array, dimension (N-1)
* The Squares of the (N-1) subdiagonal elements of T.
*
* IFIRST (input) INTEGER
* Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
* through ILAST-OFFSET elements of these arrays are to be used.
*
-* W (input/output) REAL array, dimension (N)
+* W (input/output) REAL array, dimension (N)
* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
* estimates of the eigenvalues of L D L^T indexed IFIRST through
* ILAST.
* On output, these estimates are refined.
*
-* WERR (input/output) REAL array, dimension (N)
+* WERR (input/output) REAL array, dimension (N)
* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
* the errors in the estimates of the corresponding elements in W.
* On output, these errors are refined.
*
-* WORK (workspace) REAL array, dimension (2*N)
+* WORK (workspace) REAL array, dimension (2*N)
* Workspace.
*
* IWORK (workspace) INTEGER array, dimension (2*N)
* The index of the eigenvalues to be returned.
*
* GL (input) REAL
+*
* GU (input) REAL
* An upper and a lower bound on the eigenvalue.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* The n diagonal elements of the tridiagonal matrix T.
*
-* E2 (input) REAL array, dimension (N-1)
+* E2 (input) REAL array, dimension (N-1)
* The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
*
* PIVMIN (input) REAL
* N (input) INTEGER
* The order of the matrix. N > 0.
*
-* D (input) REAL array, dimension (N)
+* D (input) REAL array, dimension (N)
* The N diagonal elements of the tridiagonal matrix T.
*
* E (input/output) REAL array, dimension (N)
* The order of the matrix. N >= 0.
*
* VL (input) REAL
+*
* VU (input) REAL
* Lower and upper bounds of the interval that contains the desired
* eigenvalues. VL < VU. Needed to compute gaps on the left or right
* The total number of input eigenvalues. 0 <= M <= N.
*
* DOL (input) INTEGER
+*
* DOU (input) INTEGER
* If the user wants to compute only selected eigenvectors from all
* the eigenvalues supplied, he can specify an index range DOL:DOU.
* MINRGP (input) REAL
*
* RTOL1 (input) REAL
+*
* RTOL2 (input) REAL
* Parameters for bisection.
* An interval [LEFT,RIGHT] has converged if
* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*
-* W (input/output) REAL array, dimension (N)
+* W (input/output) REAL array, dimension (N)
* The first M elements of W contain the APPROXIMATE eigenvalues for
* which eigenvectors are to be computed. The eigenvalues
* should be grouped by split-off block and ordered from
* 'Q' or 'Z'.
*
* CFROM (input) REAL
+*
* CTO (input) REAL
* The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
* without over/underflow if the final result CTO*A(I,J)/CFROM
* abs(SSMAX) is the larger singular value.
*
* SNL (output) REAL
+*
* CSL (output) REAL
* The vector (CSL, SNL) is a unit left singular vector for the
* singular value abs(SSMAX).
*
* SNR (output) REAL
+*
* CSR (output) REAL
* The vector (CSR, SNR) is a unit right singular vector for the
* singular value abs(SSMAX).
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SPORFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T
* user-provided factorizations and equilibration factors if they
* differ from what SPOSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(S) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* LDQ >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* corresponding elements of A.
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found.
* = 'I': the IL-th through IU-th eigenvalues will be found.
-********** For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
-********** SSTEIN are called
+* For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
+* SSTEIN are called
*
* N (input) INTEGER
* The order of the matrix. N >= 0.
* to avoid over/underflow in computing the eigenvalues.
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* indicating the nonzero elements in Z. The i-th eigenvector
* is nonzero only in elements ISUPPZ( 2*i-1 ) through
* ISUPPZ( 2*i ).
-********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
+* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
*
* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal (and
* to avoid over/underflow in computing the eigenvalues.
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
SUBROUTINE SSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
*
* -- LAPACK PROTOTYPE routine (version 3.2.2) --
-*
* -- Written by Julie Langou of the Univ. of TN --
* May 2010
*
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
-*
+*
* WAY (input) CHARACTER*1
* = 'C': Convert
* = 'R': Revert
* = 'V': all eigenvalues in the half-open interval (VL,VU]
* will be found.
* = 'I': the IL-th through IU-th eigenvalues will be found.
-********** For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
-********** SSTEIN are called
+* For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
+* SSTEIN are called
*
* UPLO (input) CHARACTER*1
* = 'U': Upper triangle of A is stored;
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* indicating the nonzero elements in Z. The i-th eigenvector
* is nonzero only in elements ISUPPZ( 2*i-1 ) through
* ISUPPZ( 2*i ).
-********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
+* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
*
* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* The leading dimension of the array B. LDB >= max(1,N).
*
* VL (input) REAL
+*
* VU (input) REAL
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SSYRFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* SSYSVXX uses the diagonal pivoting factorization to compute the
* user-provided factorizations and equilibration factors if they
* differ from what SSYSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(R) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
REAL ZERO, ONE
* If WANTZ = .TRUE., LDZ >= N.
*
* IFST (input/output) INTEGER
+*
* ILST (input/output) INTEGER
* Specify the reordering of the diagonal blocks of (A, B).
* The block with row index IFST is moved to row ILST, by a
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZGBRFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* ZGBSVXX uses the LU factorization to compute the solution to a
* complex*16 system of linear equations A * X = B, where A is an
* user-provided factorizations and equilibration factors if they
* differ from what ZGBSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by ZGEBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array A. LDA >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
* JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values determined when A was
* balanced. The balanced A(i,j) = 0 if I > J and
* u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
* are left eigenvectors of (A,B).
*
-* Note: this routine performs "full balancing" on A and B -- see
-* "Further Details", below.
+* Note: this routine performs "full balancing" on A and B
*
* Arguments
* =========
* .. Executable Statements ..
*
* Test input arguments
-* ====================
+* ====================
*
INFO = 0
LQUERY = ( LWORK.EQ.-1 )
NFXD = NFXD - 1
*
* Factorize fixed columns
-* =======================
+* =======================
*
* Compute the QR factorization of fixed columns and update
* remaining columns.
END IF
*
* Factorize free columns
-* ======================
+* ======================
*
IF( NFXD.LT.MINMN ) THEN
*
* On exit, the elements on and above the diagonal of the array
* contain the min(M,N)-by-N upper trapezoidal matrix R (R is
* upper triangular if M >= N); the elements below the diagonal
-* are the columns of V. See below for further details.
+* are the columns of V.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
* The upper triangular block reflectors stored in compact form
* as a sequence of upper triangular blocks. See below
* for further details.
-*
+*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= NB.
*
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZGERFSX improves the computed solution to a system of linear
* and C below. In this case, the solution and error bounds returned
* are for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZGESVXX uses the LU factorization to compute the solution to a
* user-provided factorizations and equilibration factors if they
* differ from what ZGESVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
* that it solves the original system before equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* The integers ILO and IHI determined by ZGGBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
* The leading dimension of the array B. LDB >= max(1,N).
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* for which SELCTG is true.
*
* ALPHA (output) COMPLEX*16 array, dimension (N)
+*
* BETA (output) COMPLEX*16 array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
* for which SELCTG is true.
*
* ALPHA (output) COMPLEX*16 array, dimension (N)
+*
* BETA (output) COMPLEX*16 array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX*16 array, dimension (N)
+*
* BETA (output) COMPLEX*16 array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues.
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX*16 array, dimension (N)
+*
* BETA (output) COMPLEX*16 array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
* eigenvalues.
* if JOBVR = 'V', LDVR >= N.
*
* ILO (output) INTEGER
+*
* IHI (output) INTEGER
* ILO and IHI are integer values such that on exit
* A(i,j) = 0 and B(i,j) = 0 if i > j and
* For further explanation of the reciprocal condition numbers RCONDE
* and RCONDV, see section 4.11 of LAPACK User's Guide.
*
+* =====================================================================
+*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
* LDQ >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* The leading dimension of the array A. LDA >= max(1,N).
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
$ ERR_BNDS_NORM( NRHS, * ),
$ ERR_BNDS_COMP( NRHS, * )
*
-* Purpose
+* Purpose
* =======
*
* ZHERFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZHESVXX uses the diagonal pivoting factorization to compute the
* user-provided factorizations and equilibration factors if they
* differ from what ZHESVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(R) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* The order of the matrices H, T, Q, and Z. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* ILO and IHI mark the rows and columns of H which are in
* Hessenberg form. It is assumed that A is already upper
* corresponding elements of A.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* If RANGE='V', the lower and upper bounds of the interval to
* be searched for eigenvalues. VL < VU.
* Not referenced if RANGE = 'A' or 'I'.
*
* IL (input) INTEGER
+*
* IU (input) INTEGER
* If RANGE='I', the indices (in ascending order) of the
* smallest and largest eigenvalues to be returned.
* .. Array Arguments ..
COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
-* Purpose
+* Purpose
* =======
*
* ZHSEQR computes the eigenvalues of a Hessenberg matrix H
* of a matrix A which has been reduced to the Hessenberg form H
* by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* JOB (input) CHARACTER*1
* = 'E': compute eigenvalues only;
* If INFO .GT. 0 and COMPZ = 'N', then Z is not
* accessed.
*
-* ================================================================
+* ================================================================
* Default values supplied by
* ILAENV(ISPEC,'ZHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
* It is suggested that these defaults be adjusted in order
* for ISPEC=16 is 0. Otherwise the default for
* ISPEC=16 is 2.
*
-* ================================================================
+* ================================================================
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
* N must be at least zero.
* Unchanged on exit.
*
-* ALPHA - DOUBLE PRECISION .
+* ALPHA (input) DOUBLE PRECISION .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
-* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
+* A (input) COMPLEX*16 array, DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
* max( 1, n ).
* Unchanged on exit.
*
-* X - COMPLEX*16 array of DIMENSION at least
+* X (input) COMPLEX*16 array, DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCX ) )
* Before entry, the incremented array X must contain the
* vector x.
* X. INCX must not be zero.
* Unchanged on exit.
*
-* BETA - DOUBLE PRECISION .
+* BETA (input) DOUBLE PRECISION .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
* where abs(Z) is the componentwise absolute value of the matrix
* or vector Z.
*
+* Arguments
+* =========
+*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
* N must be at least zero.
* Unchanged on exit.
*
-* ALPHA - DOUBLE PRECISION .
+* ALPHA (input) DOUBLE PRECISION .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
-* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
+* A (input) COMPLEX*16 array, DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
* max( 1, n ).
* Unchanged on exit.
*
-* X - COMPLEX*16 array of DIMENSION at least
+* X (input) COMPLEX*16 array, DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCX ) )
* Before entry, the incremented array X must contain the
* vector x.
* X. INCX must not be zero.
* Unchanged on exit.
*
-* BETA - DOUBLE PRECISION .
+* BETA (input) DOUBLE PRECISION .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
COMPLEX*16 X( * ), Y( * ), W( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZLA_WWADDW adds a vector W into a doubled-single vector (X, Y).
* This works for all extant IBM's hex and binary floating point
* arithmetics, but not for decimal.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* N (input) INTEGER
* The length of vectors X, Y, and W.
* ISAVE is used to save variables between calls to ZLACN2
*
* Further Details
-* ======= =======
+* ===============
*
* Contributed by Nick Higham, University of Manchester.
* Originally named CONEST, dated March 16, 1988.
* On the final return from ZLACON, KASE will again be 0.
*
* Further Details
-* ======= =======
+* ===============
*
* Contributed by Nick Higham, University of Manchester.
* Originally named CONEST, dated March 16, 1988.
* The increment between successive values of CY. INCY <> 0.
*
* C (input) COMPLEX*16
+*
* S (input) COMPLEX*16
* C and S define the matrix
* [ C S ].
* =========
*
* X (input) COMPLEX*16
+*
* Y (input) COMPLEX*16
* The complex scalars X and Y.
*
* value THRESH (set below).
*
* CS1 (output) COMPLEX*16
+*
* SN1 (output) COMPLEX*16
* If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector
* for RT1.
* The eigenvalue of smaller absolute value.
*
* CS1 (output) DOUBLE PRECISION
+*
* SN1 (output) COMPLEX*16
* The vector (CS1, SN1) is a unit right eigenvector for RT1.
*
* = .FALSE.: the input matrices A and B are lower triangular.
*
* A1 (input) DOUBLE PRECISION
+*
* A2 (input) COMPLEX*16
+*
* A3 (input) DOUBLE PRECISION
* On entry, A1, A2 and A3 are elements of the input 2-by-2
* upper (lower) triangular matrix A.
*
* B1 (input) DOUBLE PRECISION
+*
* B2 (input) COMPLEX*16
+*
* B3 (input) DOUBLE PRECISION
* On entry, B1, B2 and B3 are elements of the input 2-by-2
* upper (lower) triangular matrix B.
*
* CSU (output) DOUBLE PRECISION
+*
* SNU (output) COMPLEX*16
* The desired unitary matrix U.
*
* CSV (output) DOUBLE PRECISION
+*
* SNV (output) COMPLEX*16
* The desired unitary matrix V.
*
* CSQ (output) DOUBLE PRECISION
+*
* SNQ (output) COMPLEX*16
* The desired unitary matrix Q.
*
COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* ZLAHQR is an auxiliary routine called by CHSEQR to update the
* eigenvalues and Schur decomposition already computed by CHSEQR, by
* dealing with the Hessenberg submatrix in rows and columns ILO to
* IHI.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N >= 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows and
* columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
* of the Schur form returned in H, with W(i) = H(i,i).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE..
* where U is the orthogonal matrix in (*)
* (regardless of the value of WANTT.)
*
-* Further Details
-* ===============
+* Further Details
+* ===============
*
* 02-96 Based on modifications by
* David Day, Sandia National Laboratory, USA
* (2) adopts the more conservative Ahues & Tisseur stopping
* criterion (LAWN 122, 1997).
*
-* =========================================================
+* =========================================================
*
* .. Parameters ..
INTEGER ITMAX
COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
*
-* Purpose
-* =======
+* Purpose
+* =======
*
* ZLAQR0 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* of a matrix A which has been reduced to the Hessenberg form H
* by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
COMPLEX*16 H( LDH, * ), V( * )
* ..
*
+* Purpose
+* =======
+*
* Given a 2-by-2 or 3-by-3 matrix H, ZLAQR1 sets v to a
* scalar multiple of the first column of the product
*
* This is useful for starting double implicit shift bulges
* in the QR algorithm.
*
+* Arguments
+* =========
*
* N (input) integer
* Order of the matrix H. N must be either 2 or 3.
* A scalar multiple of the first column of the
* matrix K in (*).
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
*
* .. Parameters ..
COMPLEX*16 ZERO
$ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
* ..
*
-* This subroutine is identical to ZLAQR3 except that it avoids
-* recursion by calling ZLAHQR instead of ZLAQR4.
+* Purpose
+* =======
*
+* ZLAQR2 is identical to ZLAQR3 except that it avoids
+* recursion by calling ZLAHQR instead of ZLAQR4.
*
-* ******************************************************************
* Aggressive early deflation:
*
-* This subroutine accepts as input an upper Hessenberg matrix
+* ZLAQR2 accepts as input an upper Hessenberg matrix
* H and performs an unitary similarity transformation
* designed to detect and deflate fully converged eigenvalues from
* a trailing principal submatrix. On output H has been over-
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+*
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the triangular Schur factor may be
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
COMPLEX*16 ZERO, ONE
PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
$ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
* ..
*
-* ******************************************************************
+* Purpose
+* =======
+*
* Aggressive early deflation:
*
-* This subroutine accepts as input an upper Hessenberg matrix
+* ZLAQR3 accepts as input an upper Hessenberg matrix
* H and performs an unitary similarity transformation
* designed to detect and deflate fully converged eigenvalues from
* a trailing principal submatrix. On output H has been over-
* hoped that the final version of H has many zero subdiagonal
* entries.
*
-* ******************************************************************
+*
+* Arguments
+* =========
+*
* WANTT (input) LOGICAL
* If .TRUE., then the Hessenberg matrix H is fully updated
* so that the triangular Schur factor may be
* subroutine. N .LE. LDH
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
* in WORK(1). No error message related to LWORK is issued
* by XERBLA. Neither H nor Z are accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
COMPLEX*16 ZERO, ONE
PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
*
-* This subroutine implements one level of recursion for ZLAQR0.
+* Purpose
+* =======
+*
+* ZLAQR4 implements one level of recursion for ZLAQR0.
* It is a complete implementation of the small bulge multi-shift
* QR algorithm. It may be called by ZLAQR0 and, for large enough
* deflation window size, it may be called by ZLAQR3. This
* subroutine is identical to ZLAQR0 except that it calls ZLAQR2
* instead of ZLAQR3.
*
-* Purpose
-* =======
-*
* ZLAQR4 computes the eigenvalues of a Hessenberg matrix H
* and, optionally, the matrices T and Z from the Schur decomposition
* H = Z T Z**H, where T is an upper triangular matrix (the
* of a matrix A which has been reduced to the Hessenberg form H
* by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* WANTT (input) LOGICAL
* = .TRUE. : the full Schur form T is required;
* The order of the matrix H. N .GE. 0.
*
* ILO (input) INTEGER
+*
* IHI (input) INTEGER
* It is assumed that H is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
* accessed.
*
-* ================================================================
+* Further Details
+* ===============
+*
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
* References:
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
* Algorithm Part II: Aggressive Early Deflation, SIAM Journal
* of Matrix Analysis, volume 23, pages 948--973, 2002.
*
-* ================================================================
+* ================================================================
+*
* .. Parameters ..
*
* ==== Matrices of order NTINY or smaller must be processed by
$ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
* ..
*
-* This auxiliary subroutine called by ZLAQR0 performs a
+* Purpose
+* =======
+*
+* ZLAQR5, called by ZLAQR0, performs a
* single small-bulge multi-shift QR sweep.
*
+* Arguments
+* =========
+*
* WANTT (input) logical scalar
* WANTT = .true. if the triangular Schur factor
* is being computed. WANTT is set to .false. otherwise.
* subroutine operates.
*
* KTOP (input) integer scalar
+*
* KBOT (input) integer scalar
* These are the first and last rows and columns of an
* isolated diagonal block upon which the QR sweep is to be
* calling procedure. LDH.GE.MAX(1,N).
*
* ILOZ (input) INTEGER
+*
* IHIZ (input) INTEGER
* Specify the rows of Z to which transformations must be
* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
* LDWV is the leading dimension of WV as declared in the
* in the calling subroutine. LDWV.GE.NV.
*
-* ================================================================
+* ================================================================
* Based on contributions by
* Karen Braman and Ralph Byers, Department of Mathematics,
* University of Kansas, USA
*
-* ================================================================
+* ================================================================
* Reference:
*
* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
* Level 3 Performance, SIAM Journal of Matrix Analysis,
* volume 23, pages 929--947, 2002.
*
-* ================================================================
+* ================================================================
* .. Parameters ..
COMPLEX*16 ZERO, ONE
PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
LASTV = 0
LASTC = 0
IF( TAU.NE.ZERO ) THEN
-! Set up variables for scanning V. LASTV begins pointing to the end
-! of V.
+* Set up variables for scanning V. LASTV begins pointing to the end
+* of V.
IF( APPLYLEFT ) THEN
LASTV = M
ELSE
ELSE
I = 1
END IF
-! Look for the last non-zero row in V.
+* Look for the last non-zero row in V.
DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
LASTV = LASTV - 1
I = I - INCV
END DO
IF( APPLYLEFT ) THEN
-! Scan for the last non-zero column in C(1:lastv,:).
+* Scan for the last non-zero column in C(1:lastv,:).
LASTC = ILAZLC(LASTV, N, C, LDC)
ELSE
-! Scan for the last non-zero row in C(:,1:lastv).
+* Scan for the last non-zero row in C(:,1:lastv).
LASTC = ILAZLR(M, LASTV, C, LDC)
END IF
END IF
-! Note that lastc.eq.0 renders the BLAS operations null; no special
-! case is needed at this level.
+* Note that lastc.eq.0 renders the BLAS operations null; no special
+* case is needed at this level.
IF( APPLYLEFT ) THEN
*
* Form H * C
* The increment between elements of C. INCC > 0.
*
* Further Details
-* ======= =======
+* ===============
*
* 6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
*
* The order of the matrix. N >= 0.
*
* VL (input) DOUBLE PRECISION
+*
* VU (input) DOUBLE PRECISION
* Lower and upper bounds of the interval that contains the desired
* eigenvalues. VL < VU. Needed to compute gaps on the left or right
* The total number of input eigenvalues. 0 <= M <= N.
*
* DOL (input) INTEGER
+*
* DOU (input) INTEGER
* If the user wants to compute only selected eigenvectors from all
* the eigenvalues supplied, he can specify an index range DOL:DOU.
* MINRGP (input) DOUBLE PRECISION
*
* RTOL1 (input) DOUBLE PRECISION
+*
* RTOL2 (input) DOUBLE PRECISION
* Parameters for bisection.
* An interval [LEFT,RIGHT] has converged if
* 'Q' or 'Z'.
*
* CFROM (input) DOUBLE PRECISION
+*
* CTO (input) DOUBLE PRECISION
* The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
* without over/underflow if the final result CTO*A(I,J)/CFROM
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, ZTBSV
* is called, otherwise, specific code is used which checks for possible
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, ZTPSV
* is called, otherwise, specific code is used which checks for possible
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
-* ======= =======
+* ===============
*
* A rough bound on x is computed; if that is less than overflow, ZTRSV
* is called, otherwise, specific code is used which checks for possible
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZPORFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T
* user-provided factorizations and equilibration factors if they
* differ from what ZPOSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(S) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
SUBROUTINE ZSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
*
* -- LAPACK PROTOTYPE routine (version 3.2.2) --
-*
* -- Written by Julie Langou of the Univ. of TN --
* May 2010
*
* as an upper or lower triangular matrix.
* = 'U': Upper triangular, form is A = U*D*U**T;
* = 'L': Lower triangular, form is A = L*D*L**T.
-*
+*
* WAY (input) CHARACTER*1
* = 'C': Convert
* = 'R': Revert
$ ERR_BNDS_COMP( NRHS, * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZSYRFSX improves the computed solution to a system of linear
* below. In this case, the solution and error bounds returned are
* for the original unequilibrated system.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
* ..
*
-* Purpose
+* Purpose
* =======
*
* ZSYSVXX uses the diagonal pivoting factorization to compute the
* user-provided factorizations and equilibration factors if they
* differ from what ZSYSVXX would itself produce.
*
-* Description
-* ===========
+* Description
+* ===========
*
* The following steps are performed:
*
* diag(R) so that it solves the original system before
* equilibration.
*
-* Arguments
-* =========
+* Arguments
+* =========
*
* Some optional parameters are bundled in the PARAMS array. These
* settings determine how refinement is performed, but often the
* about all of the right-hand sides check ERR_BNDS_NORM or
* ERR_BNDS_COMP.
*
-* ==================================================================
+* ==================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
* If WANTZ = .TRUE., LDZ >= N.
*
* IFST (input) INTEGER
+*
* ILST (input/output) INTEGER
* Specify the reordering of the diagonal blocks of (A, B).
* The block with row index IFST is moved to row ILST, by a