floraachy/OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Refs #247. Included lapack source codes. Avoid downloading tar.gz from netlib.org

Browse Source 
Based on 3.4.2 version, apply patch.for_lapack-3.4.2.
This commit is contained in:
Zhang Xianyi
2013-07-09 17:00:02 +08:00
parent fbb75e58b1
commit 3eb5af1955
5320 changed files with 1448313 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

603
lapack-netlib/SRC/dtrsna.f Normal file
View File

@@ -0,0 +1,603 @@
*> \brief \b DTRSNA
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DTRSNA + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
* LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
* INFO )
*
* .. Scalar Arguments ..
* CHARACTER HOWMNY, JOB
* INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
* ..
* .. Array Arguments ..
* LOGICAL SELECT( * )
* INTEGER IWORK( * )
* DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
* $ VR( LDVR, * ), WORK( LDWORK, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTRSNA estimates reciprocal condition numbers for specified
*> eigenvalues and/or right eigenvectors of a real upper
*> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
*> orthogonal).
*>
*> T must be in Schur canonical form (as returned by DHSEQR), that is,
*> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
*> 2-by-2 diagonal block has its diagonal elements equal and its
*> off-diagonal elements of opposite sign.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] JOB
*> \verbatim
*> JOB is CHARACTER*1
*> Specifies whether condition numbers are required for
*> eigenvalues (S) or eigenvectors (SEP):
*> = 'E': for eigenvalues only (S);
*> = 'V': for eigenvectors only (SEP);
*> = 'B': for both eigenvalues and eigenvectors (S and SEP).
*> \endverbatim
*>
*> \param[in] HOWMNY
*> \verbatim
*> HOWMNY is CHARACTER*1
*> = 'A': compute condition numbers for all eigenpairs;
*> = 'S': compute condition numbers for selected eigenpairs
*> specified by the array SELECT.
*> \endverbatim
*>
*> \param[in] SELECT
*> \verbatim
*> SELECT is LOGICAL array, dimension (N)
*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
*> condition numbers are required. To select condition numbers
*> for the eigenpair corresponding to a real eigenvalue w(j),
*> SELECT(j) must be set to .TRUE.. To select condition numbers
*> corresponding to a complex conjugate pair of eigenvalues w(j)
*> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
*> set to .TRUE..
*> If HOWMNY = 'A', SELECT is not referenced.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The order of the matrix T. N >= 0.
*> \endverbatim
*>
*> \param[in] T
*> \verbatim
*> T is DOUBLE PRECISION array, dimension (LDT,N)
*> The upper quasi-triangular matrix T, in Schur canonical form.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*> LDT is INTEGER
*> The leading dimension of the array T. LDT >= max(1,N).
*> \endverbatim
*>
*> \param[in] VL
*> \verbatim
*> VL is DOUBLE PRECISION array, dimension (LDVL,M)
*> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
*> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
*> must be stored in consecutive columns of VL, as returned by
*> DHSEIN or DTREVC.
*> If JOB = 'V', VL is not referenced.
*> \endverbatim
*>
*> \param[in] LDVL
*> \verbatim
*> LDVL is INTEGER
*> The leading dimension of the array VL.
*> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
*> \endverbatim
*>
*> \param[in] VR
*> \verbatim
*> VR is DOUBLE PRECISION array, dimension (LDVR,M)
*> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
*> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
*> must be stored in consecutive columns of VR, as returned by
*> DHSEIN or DTREVC.
*> If JOB = 'V', VR is not referenced.
*> \endverbatim
*>
*> \param[in] LDVR
*> \verbatim
*> LDVR is INTEGER
*> The leading dimension of the array VR.
*> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*> S is DOUBLE PRECISION array, dimension (MM)
*> If JOB = 'E' or 'B', the reciprocal condition numbers of the
*> selected eigenvalues, stored in consecutive elements of the
*> array. For a complex conjugate pair of eigenvalues two
*> consecutive elements of S are set to the same value. Thus
*> S(j), SEP(j), and the j-th columns of VL and VR all
*> correspond to the same eigenpair (but not in general the
*> j-th eigenpair, unless all eigenpairs are selected).
*> If JOB = 'V', S is not referenced.
*> \endverbatim
*>
*> \param[out] SEP
*> \verbatim
*> SEP is DOUBLE PRECISION array, dimension (MM)
*> If JOB = 'V' or 'B', the estimated reciprocal condition
*> numbers of the selected eigenvectors, stored in consecutive
*> elements of the array. For a complex eigenvector two
*> consecutive elements of SEP are set to the same value. If
*> the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
*> is set to 0; this can only occur when the true value would be
*> very small anyway.
*> If JOB = 'E', SEP is not referenced.
*> \endverbatim
*>
*> \param[in] MM
*> \verbatim
*> MM is INTEGER
*> The number of elements in the arrays S (if JOB = 'E' or 'B')
*> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
*> \endverbatim
*>
*> \param[out] M
*> \verbatim
*> M is INTEGER
*> The number of elements of the arrays S and/or SEP actually
*> used to store the estimated condition numbers.
*> If HOWMNY = 'A', M is set to N.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6)
*> If JOB = 'E', WORK is not referenced.
*> \endverbatim
*>
*> \param[in] LDWORK
*> \verbatim
*> LDWORK is INTEGER
*> The leading dimension of the array WORK.
*> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (2*(N-1))
*> If JOB = 'E', IWORK is not referenced.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup doubleOTHERcomputational
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> The reciprocal of the condition number of an eigenvalue lambda is
*> defined as
*>
*> S(lambda) = |v**T*u| / (norm(u)*norm(v))
*>
*> where u and v are the right and left eigenvectors of T corresponding
*> to lambda; v**T denotes the transpose of v, and norm(u)
*> denotes the Euclidean norm. These reciprocal condition numbers always
*> lie between zero (very badly conditioned) and one (very well
*> conditioned). If n = 1, S(lambda) is defined to be 1.
*>
*> An approximate error bound for a computed eigenvalue W(i) is given by
*>
*> EPS * norm(T) / S(i)
*>
*> where EPS is the machine precision.
*>
*> The reciprocal of the condition number of the right eigenvector u
*> corresponding to lambda is defined as follows. Suppose
*>
*> T = ( lambda c )
*> ( 0 T22 )
*>
*> Then the reciprocal condition number is
*>
*> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
*>
*> where sigma-min denotes the smallest singular value. We approximate
*> the smallest singular value by the reciprocal of an estimate of the
*> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
*> defined to be abs(T(1,1)).
*>
*> An approximate error bound for a computed right eigenvector VR(i)
*> is given by
*>
*> EPS * norm(T) / SEP(i)
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
$ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
$ INFO )
*
* -- LAPACK computational routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
*
* .. Scalar Arguments ..
CHARACTER HOWMNY, JOB
INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
* ..
* .. Array Arguments ..
LOGICAL SELECT( * )
INTEGER IWORK( * )
DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
$ VR( LDVR, * ), WORK( LDWORK, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE, TWO
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
$ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
* ..
* .. Local Arrays ..
INTEGER ISAVE( 3 )
DOUBLE PRECISION DUMMY( 1 )
* ..
* .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
* ..
* .. External Subroutines ..
EXTERNAL DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT
* ..
* .. Executable Statements ..
*
* Decode and test the input parameters
*
WANTBH = LSAME( JOB, 'B' )
WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
*
SOMCON = LSAME( HOWMNY, 'S' )
*
INFO = 0
IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
INFO = -1
ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -4
ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
INFO = -6
ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
INFO = -8
ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
INFO = -10
ELSE
*
* Set M to the number of eigenpairs for which condition numbers
* are required, and test MM.
*
IF( SOMCON ) THEN
M = 0
PAIR = .FALSE.
DO 10 K = 1, N
IF( PAIR ) THEN
PAIR = .FALSE.
ELSE
IF( K.LT.N ) THEN
IF( T( K+1, K ).EQ.ZERO ) THEN
IF( SELECT( K ) )
$ M = M + 1
ELSE
PAIR = .TRUE.
IF( SELECT( K ) .OR. SELECT( K+1 ) )
$ M = M + 2
END IF
ELSE
IF( SELECT( N ) )
$ M = M + 1
END IF
END IF
10 CONTINUE
ELSE
M = N
END IF
*
IF( MM.LT.M ) THEN
INFO = -13
ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
INFO = -16
END IF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DTRSNA', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
IF( N.EQ.1 ) THEN
IF( SOMCON ) THEN
IF( .NOT.SELECT( 1 ) )
$ RETURN
END IF
IF( WANTS )
$ S( 1 ) = ONE
IF( WANTSP )
$ SEP( 1 ) = ABS( T( 1, 1 ) )
RETURN
END IF
*
* Get machine constants
*
EPS = DLAMCH( 'P' )
SMLNUM = DLAMCH( 'S' ) / EPS
BIGNUM = ONE / SMLNUM
CALL DLABAD( SMLNUM, BIGNUM )
*
KS = 0
PAIR = .FALSE.
DO 60 K = 1, N
*
* Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
*
IF( PAIR ) THEN
PAIR = .FALSE.
GO TO 60
ELSE
IF( K.LT.N )
$ PAIR = T( K+1, K ).NE.ZERO
END IF
*
* Determine whether condition numbers are required for the k-th
* eigenpair.
*
IF( SOMCON ) THEN
IF( PAIR ) THEN
IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
$ GO TO 60
ELSE
IF( .NOT.SELECT( K ) )
$ GO TO 60
END IF
END IF
*
KS = KS + 1
*
IF( WANTS ) THEN
*
* Compute the reciprocal condition number of the k-th
* eigenvalue.
*
IF( .NOT.PAIR ) THEN
*
* Real eigenvalue.
*
PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
RNRM = DNRM2( N, VR( 1, KS ), 1 )
LNRM = DNRM2( N, VL( 1, KS ), 1 )
S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
ELSE
*
* Complex eigenvalue.
*
PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
$ 1 )
PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
$ 1 )
RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
$ DNRM2( N, VR( 1, KS+1 ), 1 ) )
LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
$ DNRM2( N, VL( 1, KS+1 ), 1 ) )
COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
S( KS ) = COND
S( KS+1 ) = COND
END IF
END IF
*
IF( WANTSP ) THEN
*
* Estimate the reciprocal condition number of the k-th
* eigenvector.
*
* Copy the matrix T to the array WORK and swap the diagonal
* block beginning at T(k,k) to the (1,1) position.
*
CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
IFST = K
ILST = 1
CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
$ WORK( 1, N+1 ), IERR )
*
IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
*
* Could not swap because blocks not well separated
*
SCALE = ONE
EST = BIGNUM
ELSE
*
* Reordering successful
*
IF( WORK( 2, 1 ).EQ.ZERO ) THEN
*
* Form C = T22 - lambda*I in WORK(2:N,2:N).
*
DO 20 I = 2, N
WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
20 CONTINUE
N2 = 1
NN = N - 1
ELSE
*
* Triangularize the 2 by 2 block by unitary
* transformation U = [ cs i*ss ]
* [ i*ss cs ].
* such that the (1,1) position of WORK is complex
* eigenvalue lambda with positive imaginary part. (2,2)
* position of WORK is the complex eigenvalue lambda
* with negative imaginary part.
*
MU = SQRT( ABS( WORK( 1, 2 ) ) )*
$ SQRT( ABS( WORK( 2, 1 ) ) )
DELTA = DLAPY2( MU, WORK( 2, 1 ) )
CS = MU / DELTA
SN = -WORK( 2, 1 ) / DELTA
*
* Form
*
* C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
* [ mu ]
* [ .. ]
* [ .. ]
* [ mu ]
* where C**T is transpose of matrix C,
* and RWORK is stored starting in the N+1-st column of
* WORK.
*
DO 30 J = 3, N
WORK( 2, J ) = CS*WORK( 2, J )
WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
30 CONTINUE
WORK( 2, 2 ) = ZERO
*
WORK( 1, N+1 ) = TWO*MU
DO 40 I = 2, N - 1
WORK( I, N+1 ) = SN*WORK( 1, I+1 )
40 CONTINUE
N2 = 2
NN = 2*( N-1 )
END IF
*
* Estimate norm(inv(C**T))
*
EST = ZERO
KASE = 0
50 CONTINUE
CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
$ EST, KASE, ISAVE )
IF( KASE.NE.0 ) THEN
IF( KASE.EQ.1 ) THEN
IF( N2.EQ.1 ) THEN
*
* Real eigenvalue: solve C**T*x = scale*c.
*
CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
$ LDWORK, DUMMY, DUMM, SCALE,
$ WORK( 1, N+4 ), WORK( 1, N+6 ),
$ IERR )
ELSE
*
* Complex eigenvalue: solve
* C**T*(p+iq) = scale*(c+id) in real arithmetic.
*
CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
$ LDWORK, WORK( 1, N+1 ), MU, SCALE,
$ WORK( 1, N+4 ), WORK( 1, N+6 ),
$ IERR )
END IF
ELSE
IF( N2.EQ.1 ) THEN
*
* Real eigenvalue: solve C*x = scale*c.
*
CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
$ LDWORK, DUMMY, DUMM, SCALE,
$ WORK( 1, N+4 ), WORK( 1, N+6 ),
$ IERR )
ELSE
*
* Complex eigenvalue: solve
* C*(p+iq) = scale*(c+id) in real arithmetic.
*
CALL DLAQTR( .FALSE., .FALSE., N-1,
$ WORK( 2, 2 ), LDWORK,
$ WORK( 1, N+1 ), MU, SCALE,
$ WORK( 1, N+4 ), WORK( 1, N+6 ),
$ IERR )
*
END IF
END IF
*
GO TO 50
END IF
END IF
*
SEP( KS ) = SCALE / MAX( EST, SMLNUM )
IF( PAIR )
$ SEP( KS+1 ) = SEP( KS )
END IF
*
IF( PAIR )
$ KS = KS + 1
*
60 CONTINUE
RETURN
*
* End of DTRSNA
*
END