--- rpl/lapack/lapack/dstein.f 2010/12/21 13:53:38 1.7
+++ rpl/lapack/lapack/dstein.f 2023/08/07 08:39:07 1.18
@@ -1,10 +1,180 @@
+*> \brief \b DSTEIN
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DSTEIN + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
+* IWORK, IFAIL, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, LDZ, M, N
+* ..
+* .. Array Arguments ..
+* INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
+* $ IWORK( * )
+* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DSTEIN computes the eigenvectors of a real symmetric tridiagonal
+*> matrix T corresponding to specified eigenvalues, using inverse
+*> iteration.
+*>
+*> The maximum number of iterations allowed for each eigenvector is
+*> specified by an internal parameter MAXITS (currently set to 5).
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] D
+*> \verbatim
+*> D is DOUBLE PRECISION array, dimension (N)
+*> The n diagonal elements of the tridiagonal matrix T.
+*> \endverbatim
+*>
+*> \param[in] E
+*> \verbatim
+*> E is DOUBLE PRECISION array, dimension (N-1)
+*> The (n-1) subdiagonal elements of the tridiagonal matrix
+*> T, in elements 1 to N-1.
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*> M is INTEGER
+*> The number of eigenvectors to be found. 0 <= M <= N.
+*> \endverbatim
+*>
+*> \param[in] W
+*> \verbatim
+*> W is DOUBLE PRECISION array, dimension (N)
+*> The first M elements of W contain the eigenvalues for
+*> which eigenvectors are to be computed. The eigenvalues
+*> should be grouped by split-off block and ordered from
+*> smallest to largest within the block. ( The output array
+*> W from DSTEBZ with ORDER = 'B' is expected here. )
+*> \endverbatim
+*>
+*> \param[in] IBLOCK
+*> \verbatim
+*> IBLOCK is INTEGER array, dimension (N)
+*> The submatrix indices associated with the corresponding
+*> eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
+*> the first submatrix from the top, =2 if W(i) belongs to
+*> the second submatrix, etc. ( The output array IBLOCK
+*> from DSTEBZ is expected here. )
+*> \endverbatim
+*>
+*> \param[in] ISPLIT
+*> \verbatim
+*> ISPLIT is INTEGER array, dimension (N)
+*> The splitting points, at which T breaks up into submatrices.
+*> The first submatrix consists of rows/columns 1 to
+*> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
+*> through ISPLIT( 2 ), etc.
+*> ( The output array ISPLIT from DSTEBZ is expected here. )
+*> \endverbatim
+*>
+*> \param[out] Z
+*> \verbatim
+*> Z is DOUBLE PRECISION array, dimension (LDZ, M)
+*> The computed eigenvectors. The eigenvector associated
+*> with the eigenvalue W(i) is stored in the i-th column of
+*> Z. Any vector which fails to converge is set to its current
+*> iterate after MAXITS iterations.
+*> \endverbatim
+*>
+*> \param[in] LDZ
+*> \verbatim
+*> LDZ is INTEGER
+*> The leading dimension of the array Z. LDZ >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (5*N)
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] IFAIL
+*> \verbatim
+*> IFAIL is INTEGER array, dimension (M)
+*> On normal exit, all elements of IFAIL are zero.
+*> If one or more eigenvectors fail to converge after
+*> MAXITS iterations, then their indices are stored in
+*> array IFAIL.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit.
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> > 0: if INFO = i, then i eigenvectors failed to converge
+*> in MAXITS iterations. Their indices are stored in
+*> array IFAIL.
+*> \endverbatim
+*
+*> \par Internal Parameters:
+* =========================
+*>
+*> \verbatim
+*> MAXITS INTEGER, default = 5
+*> The maximum number of iterations performed.
+*>
+*> EXTRA INTEGER, default = 2
+*> The number of iterations performed after norm growth
+*> criterion is satisfied, should be at least 1.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleOTHERcomputational
+*
+* =====================================================================
SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
$ IWORK, IFAIL, INFO )
*
-* -- LAPACK routine (version 3.2) --
+* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDZ, M, N
@@ -15,89 +185,6 @@
DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
* ..
*
-* Purpose
-* =======
-*
-* DSTEIN computes the eigenvectors of a real symmetric tridiagonal
-* matrix T corresponding to specified eigenvalues, using inverse
-* iteration.
-*
-* The maximum number of iterations allowed for each eigenvector is
-* specified by an internal parameter MAXITS (currently set to 5).
-*
-* Arguments
-* =========
-*
-* N (input) INTEGER
-* The order of the matrix. N >= 0.
-*
-* D (input) DOUBLE PRECISION array, dimension (N)
-* The n diagonal elements of the tridiagonal matrix T.
-*
-* E (input) DOUBLE PRECISION array, dimension (N-1)
-* The (n-1) subdiagonal elements of the tridiagonal matrix
-* T, in elements 1 to N-1.
-*
-* M (input) INTEGER
-* The number of eigenvectors to be found. 0 <= M <= N.
-*
-* W (input) DOUBLE PRECISION array, dimension (N)
-* The first M elements of W contain the eigenvalues for
-* which eigenvectors are to be computed. The eigenvalues
-* should be grouped by split-off block and ordered from
-* smallest to largest within the block. ( The output array
-* W from DSTEBZ with ORDER = 'B' is expected here. )
-*
-* IBLOCK (input) INTEGER array, dimension (N)
-* The submatrix indices associated with the corresponding
-* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
-* the first submatrix from the top, =2 if W(i) belongs to
-* the second submatrix, etc. ( The output array IBLOCK
-* from DSTEBZ is expected here. )
-*
-* ISPLIT (input) INTEGER array, dimension (N)
-* The splitting points, at which T breaks up into submatrices.
-* The first submatrix consists of rows/columns 1 to
-* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
-* through ISPLIT( 2 ), etc.
-* ( The output array ISPLIT from DSTEBZ is expected here. )
-*
-* Z (output) DOUBLE PRECISION array, dimension (LDZ, M)
-* The computed eigenvectors. The eigenvector associated
-* with the eigenvalue W(i) is stored in the i-th column of
-* Z. Any vector which fails to converge is set to its current
-* iterate after MAXITS iterations.
-*
-* LDZ (input) INTEGER
-* The leading dimension of the array Z. LDZ >= max(1,N).
-*
-* WORK (workspace) DOUBLE PRECISION array, dimension (5*N)
-*
-* IWORK (workspace) INTEGER array, dimension (N)
-*
-* IFAIL (output) INTEGER array, dimension (M)
-* On normal exit, all elements of IFAIL are zero.
-* If one or more eigenvectors fail to converge after
-* MAXITS iterations, then their indices are stored in
-* array IFAIL.
-*
-* INFO (output) INTEGER
-* = 0: successful exit.
-* < 0: if INFO = -i, the i-th argument had an illegal value
-* > 0: if INFO = i, then i eigenvectors failed to converge
-* in MAXITS iterations. Their indices are stored in
-* array IFAIL.
-*
-* Internal Parameters
-* ===================
-*
-* MAXITS INTEGER, default = 5
-* The maximum number of iterations performed.
-*
-* EXTRA INTEGER, default = 2
-* The number of iterations performed after norm growth
-* criterion is satisfied, should be at least 1.
-*
* =====================================================================
*
* .. Parameters ..
@@ -119,8 +206,8 @@
* ..
* .. External Functions ..
INTEGER IDAMAX
- DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2
- EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DNRM2
+ DOUBLE PRECISION DDOT, DLAMCH, DNRM2
+ EXTERNAL IDAMAX, DDOT, DLAMCH, DNRM2
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL,
@@ -207,7 +294,7 @@
BLKSIZ = BN - B1 + 1
IF( BLKSIZ.EQ.1 )
$ GO TO 60
- GPIND = B1
+ GPIND = J1
*
* Compute reorthogonalization criterion and stopping criterion.
*
@@ -280,9 +367,10 @@
*
* Normalize and scale the righthand side vector Pb.
*
+ JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
SCL = BLKSIZ*ONENRM*MAX( EPS,
$ ABS( WORK( INDRV4+BLKSIZ ) ) ) /
- $ DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
+ $ ABS( WORK( INDRV1+JMAX ) )
CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
*
* Solve the system LU = Pb.