dtgevc.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 1,038 行 · 第 1/5 页

     $                   ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
               IF( TEMP.GT.BIGNUM*XSCALE ) THEN
                  DO 90 JW = 0, NW - 1
                     DO 80 JR = JE, J - 1
                        WORK( ( JW+2 )*N+JR ) = XSCALE*
     $                     WORK( ( JW+2 )*N+JR )
   80                CONTINUE
   90             CONTINUE
                  XMAX = XMAX*XSCALE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute dot products
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    j-1
</span><span class="comment">*</span><span class="comment">              SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
</span><span class="comment">*</span><span class="comment">                    k=je
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              To reduce the op count, this is done as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              _        j-1                  _        j-1
</span><span class="comment">*</span><span class="comment">              a*conjg( sum  S(k,j)*x(k) ) - b*conjg( sum  P(k,j)*x(k) )
</span><span class="comment">*</span><span class="comment">                       k=je                          k=je
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              which may cause underflow problems if A or B are close
</span><span class="comment">*</span><span class="comment">              to underflow.  (E.g., less than SMALL.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A series of compiler directives to defeat vectorization
</span><span class="comment">*</span><span class="comment">              for the next loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">$PL$ CMCHAR=' '
</span>CDIR$          NEXTSCALAR
C$DIR          SCALAR
CDIR$          NEXT SCALAR
CVD$L          NOVECTOR
CDEC$          NOVECTOR
CVD$           NOVECTOR
<span class="comment">*</span><span class="comment">VDIR          NOVECTOR
</span><span class="comment">*</span><span class="comment">VOCL          LOOP,SCALAR
</span>CIBM           PREFER SCALAR
<span class="comment">*</span><span class="comment">$PL$ CMCHAR='*'
</span><span class="comment">*</span><span class="comment">
</span>               DO 120 JW = 1, NW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">$PL$ CMCHAR=' '
</span>CDIR$             NEXTSCALAR
C$DIR             SCALAR
CDIR$             NEXT SCALAR
CVD$L             NOVECTOR
CDEC$             NOVECTOR
CVD$              NOVECTOR
<span class="comment">*</span><span class="comment">VDIR             NOVECTOR
</span><span class="comment">*</span><span class="comment">VOCL             LOOP,SCALAR
</span>CIBM              PREFER SCALAR
<span class="comment">*</span><span class="comment">$PL$ CMCHAR='*'
</span><span class="comment">*</span><span class="comment">
</span>                  DO 110 JA = 1, NA
                     SUMS( JA, JW ) = ZERO
                     SUMP( JA, JW ) = ZERO
<span class="comment">*</span><span class="comment">
</span>                     DO 100 JR = JE, J - 1
                        SUMS( JA, JW ) = SUMS( JA, JW ) +
     $                                   S( JR, J+JA-1 )*
     $                                   WORK( ( JW+1 )*N+JR )
                        SUMP( JA, JW ) = SUMP( JA, JW ) +
     $                                   P( JR, J+JA-1 )*
     $                                   WORK( ( JW+1 )*N+JR )
  100                CONTINUE
  110             CONTINUE
  120          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">$PL$ CMCHAR=' '
</span>CDIR$          NEXTSCALAR
C$DIR          SCALAR
CDIR$          NEXT SCALAR
CVD$L          NOVECTOR
CDEC$          NOVECTOR
CVD$           NOVECTOR
<span class="comment">*</span><span class="comment">VDIR          NOVECTOR
</span><span class="comment">*</span><span class="comment">VOCL          LOOP,SCALAR
</span>CIBM           PREFER SCALAR
<span class="comment">*</span><span class="comment">$PL$ CMCHAR='*'
</span><span class="comment">*</span><span class="comment">
</span>               DO 130 JA = 1, NA
                  IF( ILCPLX ) THEN
                     SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
     $                              BCOEFR*SUMP( JA, 1 ) -
     $                              BCOEFI*SUMP( JA, 2 )
                     SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
     $                              BCOEFR*SUMP( JA, 2 ) +
     $                              BCOEFI*SUMP( JA, 1 )
                  ELSE
                     SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
     $                              BCOEFR*SUMP( JA, 1 )
                  END IF
  130          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                                  T
</span><span class="comment">*</span><span class="comment">              Solve  ( a A - b B )  y = SUM(,)
</span><span class="comment">*</span><span class="comment">              with scaling and perturbation of the denominator
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLALN2.707"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
     $                      BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
     $                      BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
     $                      IINFO )
               IF( SCALE.LT.ONE ) THEN
                  DO 150 JW = 0, NW - 1
                     DO 140 JR = JE, J - 1
                        WORK( ( JW+2 )*N+JR ) = SCALE*
     $                     WORK( ( JW+2 )*N+JR )
  140                CONTINUE
  150             CONTINUE
                  XMAX = SCALE*XMAX
               END IF
               XMAX = MAX( XMAX, TEMP )
  160       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy eigenvector to VL, back transforming if
</span><span class="comment">*</span><span class="comment">           HOWMNY='B'.
</span><span class="comment">*</span><span class="comment">
</span>            IEIG = IEIG + 1
            IF( ILBACK ) THEN
               DO 170 JW = 0, NW - 1
                  CALL DGEMV( <span class="string">'N'</span>, N, N+1-JE, ONE, VL( 1, JE ), LDVL,
     $                        WORK( ( JW+2 )*N+JE ), 1, ZERO,
     $                        WORK( ( JW+4 )*N+1 ), 1 )
  170          CONTINUE
               CALL <a name="DLACPY.733"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">' '</span>, N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
     $                      LDVL )
               IBEG = 1
            ELSE
               CALL <a name="DLACPY.737"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">' '</span>, N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
     $                      LDVL )
               IBEG = JE
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Scale eigenvector
</span><span class="comment">*</span><span class="comment">
</span>            XMAX = ZERO
            IF( ILCPLX ) THEN
               DO 180 J = IBEG, N
                  XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
     $                   ABS( VL( J, IEIG+1 ) ) )
  180          CONTINUE
            ELSE
               DO 190 J = IBEG, N
                  XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
  190          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( XMAX.GT.SAFMIN ) THEN
               XSCALE = ONE / XMAX
<span class="comment">*</span><span class="comment">
</span>               DO 210 JW = 0, NW - 1
                  DO 200 JR = IBEG, N
                     VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
  200             CONTINUE
  210          CONTINUE
            END IF
            IEIG = IEIG + NW - 1
<span class="comment">*</span><span class="comment">
</span>  220    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Right eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>      IF( COMPR ) THEN
         IEIG = IM + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Main loop over eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         ILCPLX = .FALSE.
         DO 500 JE = N, 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
</span><span class="comment">*</span><span class="comment">           (b) this would be the second of a complex pair.
</span><span class="comment">*</span><span class="comment">           Check for complex eigenvalue, so as to be sure of which
</span><span class="comment">*</span><span class="comment">           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
</span><span class="comment">*</span><span class="comment">           or SELECT(JE-1).
</span><span class="comment">*</span><span class="comment">           If this is a complex pair, the 2-by-2 diagonal block
</span><span class="comment">*</span><span class="comment">           corresponding to the eigenvalue is in rows/columns JE-1:JE
</span><span class="comment">*</span><span class="comment">
</span>            IF( ILCPLX ) THEN
               ILCPLX = .FALSE.
               GO TO 500
            END IF
            NW = 1
            IF( JE.GT.1 ) THEN
               IF( S( JE, JE-1 ).NE.ZERO ) THEN
                  ILCPLX = .TRUE.
                  NW = 2
               END IF
            END IF
            IF( ILALL ) THEN
               ILCOMP = .TRUE.
            ELSE IF( ILCPLX ) THEN
               ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
            ELSE
               ILCOMP = SELECT( JE )
            END IF
            IF( .NOT.ILCOMP )
     $         GO TO 500
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Decide if (a) singular pencil, (b) real eigenvalue, or
</span><span class="comment">*</span><span class="comment">           (c) complex eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>            IF( .NOT.ILCPLX ) THEN
               IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.