clarrv.f

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

            DO 150 I = 1, OLDNCL
               IF( LSBDPT.EQ.0 ) THEN
                  OLDCLS = IINDC1
                  NEWCLS = IINDC2
               ELSE
                  OLDCLS = IINDC2
                  NEWCLS = IINDC1
               END IF
*
*              If NDEPTH > 1, retrieve the relatively robust
*              representation (RRR) and perform limited bisection
*              (if necessary) to get approximate eigenvalues.
*
               J = OLDCLS + 2*I
               OLDFST = IWORK( J-1 )
               OLDLST = IWORK( J )
               IF( NDEPTH.GT.0 ) THEN
                  J = OLDIEN + OLDFST
                  DO 45 K = 1, IN
                     D( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1,
     $                                 OLDIEN+OLDFST ) )
                     L( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1,
     $                                 OLDIEN+OLDFST+1 ) )
   45             CONTINUE
                  SIGMA = L( IEND )
               END IF
               K = IBEGIN
               OPS = OPS + REAL( 2*(IN-1) )
               DO 50 J = 1, IN - 1
                  WORK( INDLD+J ) = D( K )*L( K )
                  WORK( INDLLD+J ) = WORK( INDLD+J )*L( K )
                  K = K + 1
   50          CONTINUE
               IF( NDEPTH.GT.0 ) THEN
                  CALL SLARRB( IN, D( IBEGIN ), L( IBEGIN ),
     $                         WORK( INDLD+1 ), WORK( INDLLD+1 ),
     $                         OLDFST, OLDLST, SIGMA, RELTOL, WORK,
     $                         WORK( INDGAP+1 ), WORK( INDERR ),
     $                         WORK( INDWRK ), IWORK( IINDWK ), IINFO )
                  IF( IINFO.NE.0 ) THEN
                     INFO = 1
                     RETURN
                  END IF
               END IF
*
*              Classify eigenvalues of the current representation (RRR)
*              as (i) isolated, (ii) loosely clustered or (iii) tightly
*              clustered
*
               NEWFRS = OLDFST
               DO 140 J = OLDFST, OLDLST
                  OPS = OPS + REAL( 1 )
                  IF( J.EQ.OLDLST .OR. WORK( INDGAP+J ).GE.RELTOL*
     $                ABS( WORK( J ) ) ) THEN
                     NEWLST = J
                  ELSE
*
*                    continue (to the next loop)
*
                     OPS = OPS + REAL( 1 )
                     RELGAP = WORK( INDGAP+J ) / ABS( WORK( J ) )
                     IF( J.EQ.NEWFRS ) THEN
                        MINRGP = RELGAP
                     ELSE
                        MINRGP = MIN( MINRGP, RELGAP )
                     END IF
                     GO TO 140
                  END IF
                  NEWSIZ = NEWLST - NEWFRS + 1
                  MAXITR = 10
                  NEWFTT = OLDIEN + NEWFRS
                  IF( NEWSIZ.GT.1 ) THEN
                     MGSCLS = NEWSIZ.LE.MGSSIZ .AND. MINRGP.GE.MGSTOL
                     IF( .NOT.MGSCLS ) THEN
                        DO 55 K = 1, IN
                           WORK( INDIN1+K-1 ) = REAL( Z( IBEGIN+K-1,
     $                                          NEWFTT ) )
                           WORK( INDIN2+K-1 ) = REAL( Z( IBEGIN+K-1,
     $                                          NEWFTT+1 ) )
   55                   CONTINUE
                        CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ),
     $                               WORK( INDLD+1 ), WORK( INDLLD+1 ),
     $                               NEWFRS, NEWLST, WORK,
     $                               WORK( INDIN1 ), WORK( INDIN2 ),
     $                               WORK( INDWRK ), IWORK( IINDWK ),
     $                               INFO )
                        IF( INFO.EQ.0 ) THEN
                           NCLUS = NCLUS + 1
                           K = NEWCLS + 2*NCLUS
                           IWORK( K-1 ) = NEWFRS
                           IWORK( K ) = NEWLST
                        ELSE
                           INFO = 0
                           IF( MINRGP.GE.MGSTOL ) THEN
                              MGSCLS = .TRUE.
                           ELSE
*                           
*                             Call CSTEIN to process this tight cluster.
*                             This happens only if MINRGP <= MGSTOL
*                             and SLARRF returns INFO = 1. The latter
*                             means that a new RRR to "break" the
*                             cluster could not be found.
*
                              WORK( INDWRK ) = D( IBEGIN )
                              OPS = OPS + REAL( IN-1 )
                              DO 60 K = 1, IN - 1
                                 WORK( INDWRK+K ) = D( IBEGIN+K ) +
     $                                              WORK( INDLLD+K )
   60                         CONTINUE
                              DO 70 K = 1, NEWSIZ
                                 IWORK( IINDWK+K-1 ) = 1
   70                         CONTINUE
                              DO 80 K = NEWFRS, NEWLST
                                 ISUPPZ( 2*( IBEGIN+K )-3 ) = 1
                                 ISUPPZ( 2*( IBEGIN+K )-2 ) = IN
   80                         CONTINUE
                              TEMP( 1 ) = IN
                              CALL CSTEIN( IN, WORK( INDWRK ),
     $                                     WORK( INDLD+1 ), NEWSIZ,
     $                                     WORK( NEWFRS ),
     $                                     IWORK( IINDWK ), TEMP( 1 ),
     $                                     Z( IBEGIN, NEWFTT ), LDZ,
     $                                     WORK( INDWRK+IN ),
     $                                     IWORK( IINDWK+IN ),
     $                                     IWORK( IINDWK+2*IN ), IINFO )
                              IF( IINFO.NE.0 ) THEN
                                 INFO = 2
                                 RETURN
                              END IF
                              NDONE = NDONE + NEWSIZ
                           END IF
                        END IF
                     END IF
                  ELSE
                     MGSCLS = .FALSE.
                  END IF
                  IF( NEWSIZ.EQ.1 .OR. MGSCLS ) THEN
                     KTOT = NEWFTT
                     DO 100 K = NEWFRS, NEWLST
                        ITER = 0
   90                   CONTINUE
                        LAMBDA = WORK( K )
                        CALL CLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
     $                               L( IBEGIN ), WORK( INDLD+1 ),
     $                               WORK( INDLLD+1 ),
     $                               GERSCH( 2*OLDIEN+1 ),
     $                               Z( IBEGIN, KTOT ), ZTZ, MINGMA,
     $                               IWORK( IINDR+KTOT ),
     $                               ISUPPZ( 2*KTOT-1 ),
     $                               WORK( INDWRK ) )
                        OPS = OPS + REAL( 4 )
                        TMP1 = ONE / ZTZ
                        NRMINV = SQRT( TMP1 )
                        RESID = ABS( MINGMA )*NRMINV
                        RQCORR = MINGMA*TMP1
                        IF( K.EQ.IN ) THEN
                           GAP = WORK( INDGAP+K-1 )
                        ELSE IF( K.EQ.1 ) THEN
                           GAP = WORK( INDGAP+K )
                        ELSE
                           GAP = MIN( WORK( INDGAP+K-1 ),
     $                           WORK( INDGAP+K ) )
                        END IF
                        ITER = ITER + 1
                        OPS = OPS + REAL( 3 )
                        IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
     $                      FOUR*EPS*ABS( LAMBDA ) ) THEN
                           OPS = OPS + REAL( 1 )
                           WORK( K ) = LAMBDA + RQCORR
                           IF( ITER.LT.MAXITR ) THEN
                              GO TO 90
                           END IF
                        END IF
                        IWORK( KTOT ) = 1
                        IF( NEWSIZ.EQ.1 )
     $                     NDONE = NDONE + 1
                        OPS = OPS + REAL( 2*IN )
                        CALL CSSCAL( IN, NRMINV, Z( IBEGIN, KTOT ), 1 )
                        KTOT = KTOT + 1
  100                CONTINUE
                     IF( NEWSIZ.GT.1 ) THEN
                        ITMP1 = ISUPPZ( 2*NEWFTT-1 )
                        ITMP2 = ISUPPZ( 2*NEWFTT )
                        KTOT = OLDIEN + NEWLST
                        DO 120 P = NEWFTT + 1, KTOT
                           DO 110 Q = NEWFTT, P - 1
                              OPS = OPS + REAL( 10*IN )
                              TMP1 = -CDOTU( IN, Z( IBEGIN, P ), 1,
     $                               Z( IBEGIN, Q ), 1 )
                              CALL CAXPY( IN, CMPLX( TMP1, ZERO ),
     $                                    Z( IBEGIN, Q ), 1,
     $                                    Z( IBEGIN, P ), 1 )
  110                      CONTINUE
                           OPS = OPS + REAL( 8*IN+1 )
                           TMP1 = ONE / SCNRM2( IN, Z( IBEGIN, P ), 1 )
                           CALL CSSCAL( IN, TMP1, Z( IBEGIN, P ), 1 )
                           ITMP1 = MIN( ITMP1, ISUPPZ( 2*P-1 ) )
                           ITMP2 = MAX( ITMP2, ISUPPZ( 2*P ) )
  120                   CONTINUE
                        DO 130 P = NEWFTT, KTOT
                           ISUPPZ( 2*P-1 ) = ITMP1
                           ISUPPZ( 2*P ) = ITMP2
  130                   CONTINUE
                        NDONE = NDONE + NEWSIZ
                     END IF
                  END IF
                  NEWFRS = J + 1
  140          CONTINUE
  150       CONTINUE
            NDEPTH = NDEPTH + 1
            GO TO 40
         END IF
         J = 2*IBEGIN
         DO 160 I = IBEGIN, IEND
            ISUPPZ( J-1 ) = ISUPPZ( J-1 ) + OLDIEN
            ISUPPZ( J ) = ISUPPZ( J ) + OLDIEN
            J = J + 2
  160    CONTINUE
         IBEGIN = IEND + 1
  170 CONTINUE
*
      RETURN
*
*     End of CLARRV
*
      END