clarrv.f

计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．

      SUBROUTINE CLARRV( N, D, L, ISPLIT, M, W, IBLOCK, GERSCH, TOL, Z,
     $                   LDZ, ISUPPZ, WORK, IWORK, INFO )
*
*  -- LAPACK auxiliary routine (instru to count ops, version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDZ, M, N
      REAL               TOL
*     ..
*     .. Array Arguments ..
      INTEGER            IBLOCK( * ), ISPLIT( * ), ISUPPZ( * ),
     $                   IWORK( * )
      REAL               D( * ), GERSCH( * ), L( * ), W( * ), WORK( * )
      COMPLEX            Z( LDZ, * )
*     ..
*     Common block to return operation count ..
*     .. Common blocks ..
      COMMON             / LATIME / OPS, ITCNT
*     ..
*     .. Scalars in Common ..
      REAL               ITCNT, OPS
*     ..
*
*  Purpose
*  =======
*
*  CLARRV computes the eigenvectors of the tridiagonal matrix
*  T = L D L^T given L, D and the eigenvalues of L D L^T.
*  The input eigenvalues should have high relative accuracy with
*  respect to the entries of L and D. The desired accuracy of the
*  output can be specified by the input parameter TOL.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix.  N >= 0.
*
*  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-1)
*          On entry, the (n-1) subdiagonal elements of the unit
*          bidiagonal matrix L in elements 1 to N-1 of L. L(N) need
*          not be set. On exit, L is overwritten.
*
*  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.
*
*  TOL     (input) REAL
*          The absolute error tolerance for the
*          eigenvalues/eigenvectors.
*          Errors in the input eigenvalues must be bounded by TOL.
*          The eigenvectors output have residual norms
*          bounded by TOL, and the dot products between different
*          eigenvectors are bounded by TOL. TOL must be at least
*          N*EPS*|T|, where EPS is the machine precision and |T| is
*          the 1-norm of the tridiagonal matrix.
*
*  M       (input) INTEGER
*          The total number of eigenvalues found.  0 <= M <= N.
*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
*
*  W       (input) REAL 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 SLARRE is expected here ).
*          Errors in W must be bounded by TOL (see above).
*
*  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. 
*
*  Z       (output) COMPLEX array, dimension (LDZ, max(1,M) )
*          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
*          contain the orthonormal eigenvectors of the matrix T
*          corresponding to the selected eigenvalues, with the i-th
*          column of Z holding the eigenvector associated with W(i).
*          If JOBZ = 'N', then Z is not referenced.
*          Note: the user must ensure that at least max(1,M) columns are
*          supplied in the array Z; if RANGE = 'V', the exact value of M
*          is not known in advance and an upper bound must be used.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.  LDZ >= 1, and if
*          JOBZ = 'V', LDZ >= max(1,N).
*
*  ISUPPZ  (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
*          The support of the eigenvectors in Z, i.e., the indices
*          indicating the nonzero elements in Z. The i-th eigenvector
*          is nonzero only in elements ISUPPZ( 2*i-1 ) through
*          ISUPPZ( 2*i ).
*
*  WORK    (workspace) REAL array, dimension (13*N)
*
*  IWORK   (workspace) INTEGER array, dimension (6*N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = 1, internal error in SLARRB
*                if INFO = 2, internal error in CSTEIN
*
*  Further Details
*  ===============
*
*  Based on contributions by
*     Inderjit Dhillon, IBM Almaden, USA
*     Osni Marques, LBNL/NERSC, USA
*     Ken Stanley, Computer Science Division, University of
*       California at Berkeley, USA
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            MGSSIZ
      PARAMETER          ( MGSSIZ = 20 )
      REAL               ZERO, ONE, FOUR
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, FOUR = 4.0E0 )
      COMPLEX            CZERO
      PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            MGSCLS
      INTEGER            I, IBEGIN, IEND, IINDC1, IINDC2, IINDR, IINDWK, 
     $                   IINFO, IM, IN, INDERR, INDGAP, INDIN1, INDIN2, 
     $                   INDLD, INDLLD, INDWRK, ITER, ITMP1, ITMP2, J, 
     $                   JBLK, K, KTOT, LSBDPT, MAXITR, NCLUS, NDEPTH, 
     $                   NDONE, NEWCLS, NEWFRS, NEWFTT, NEWLST, NEWSIZ, 
     $                   NSPLIT, OLDCLS, OLDFST, OLDIEN, OLDLST, OLDNCL, 
     $                   P, Q, TEMP( 1 )
      REAL               EPS, GAP, LAMBDA, MGSTOL, MINGMA, MINRGP, 
     $                   NRMINV, RELGAP, RELTOL, RESID, RQCORR, SIGMA, 
     $                   TMP1, ZTZ
*     ..
*     .. External Functions ..
      REAL               SCNRM2, SLAMCH
      COMPLEX            CDOTU
      EXTERNAL           CDOTU, SCNRM2, SLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           CAXPY, CLAR1V, CLASET, CSTEIN, SCOPY, SLARRB
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, CMPLX, MAX, MIN, REAL, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INDERR = N + 1
      INDLD = 2*N
      INDLLD = 3*N
      INDGAP = 4*N
      INDIN1 = 5*N + 1
      INDIN2 = 6*N + 1
      INDWRK = 7*N + 1
*
      IINDR = N
      IINDC1 = 2*N
      IINDC2 = 3*N
      IINDWK = 4*N + 1
*
      EPS = SLAMCH( 'Precision' )
*
      DO 10 I = 1, 2*N
         IWORK( I ) = 0
   10 CONTINUE
      OPS = OPS + REAL( M+1 )
      DO 20 I = 1, M
         WORK( INDERR+I-1 ) = EPS*ABS( W( I ) )
   20 CONTINUE
      CALL CLASET( 'Full', N, N, CZERO, CZERO, Z, LDZ )
      MGSTOL = 5.0E0*EPS
*
      NSPLIT = IBLOCK( M )
      IBEGIN = 1
      DO 170 JBLK = 1, NSPLIT
         IEND = ISPLIT( JBLK )
*
*        Find the eigenvectors of the submatrix indexed IBEGIN
*        through IEND.
*
         IF( IBEGIN.EQ.IEND ) THEN
            Z( IBEGIN, IBEGIN ) = ONE
            ISUPPZ( 2*IBEGIN-1 ) = IBEGIN
            ISUPPZ( 2*IBEGIN ) = IBEGIN
            IBEGIN = IEND + 1
            GO TO 170
         END IF
         OLDIEN = IBEGIN - 1
         IN = IEND - OLDIEN
         OPS = OPS + REAL( 1 )
         RELTOL = MIN( 1.0E-2, ONE / REAL( IN ) )
         IM = IN
         CALL SCOPY( IM, W( IBEGIN ), 1, WORK, 1 )
         OPS = OPS + REAL( IN-1 )
         DO 30 I = 1, IN - 1
            WORK( INDGAP+I ) = WORK( I+1 ) - WORK( I )
   30    CONTINUE
         WORK( INDGAP+IN ) = MAX( ABS( WORK( IN ) ), EPS )
         NDONE = 0
*
         NDEPTH = 0
         LSBDPT = 1
         NCLUS = 1
         IWORK( IINDC1+1 ) = 1
         IWORK( IINDC1+2 ) = IN
*
*        While( NDONE.LT.IM ) do
*
   40    CONTINUE
         IF( NDONE.LT.IM ) THEN
            OLDNCL = NCLUS
            NCLUS = 0
            LSBDPT = 1 - LSBDPT