slarrv.f
来自「计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．」· F 代码 · 共 440 行 · 第 1/2 页
F
440 行
OLDNCL = NCLUS NCLUS = 0 LSBDPT = 1 - LSBDPT 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 CALL SCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 ) CALL SCOPY( IN, Z( IBEGIN, J+1 ), 1, L( IBEGIN ), 1 ) 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 CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ), $ WORK( INDLD+1 ), WORK( INDLLD+1 ), $ NEWFRS, NEWLST, WORK, $ Z( IBEGIN, NEWFTT ), $ Z( IBEGIN, NEWFTT+1 ), $ 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 SSTEIN 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 SSTEIN( 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 SLAR1V( 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( IN ) CALL SSCAL( 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( 4*IN ) TMP1 = -SDOT( IN, Z( IBEGIN, P ), 1, $ Z( IBEGIN, Q ), 1 ) CALL SAXPY( IN, TMP1, Z( IBEGIN, Q ), 1, $ Z( IBEGIN, P ), 1 ) 110 CONTINUE OPS = OPS + REAL( 3*IN+1 ) TMP1 = ONE / SNRM2( IN, Z( IBEGIN, P ), 1 ) CALL SSCAL( 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 SLARRV* END⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?