dbdsdc.f.html

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

      CALL <a name="DLASCL.286"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
      CALL <a name="DLASCL.287"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.289"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span>      MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
      SMLSZP = SMLSIZ + 1
<span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.1 ) THEN
         IU = 1
         IVT = 1 + SMLSIZ
         DIFL = IVT + SMLSZP
         DIFR = DIFL + MLVL
         Z = DIFR + MLVL*2
         IC = Z + MLVL
         IS = IC + 1
         POLES = IS + 1
         GIVNUM = POLES + 2*MLVL
<span class="comment">*</span><span class="comment">
</span>         K = 1
         GIVPTR = 2
         PERM = 3
         GIVCOL = PERM + MLVL
      END IF
<span class="comment">*</span><span class="comment">
</span>      DO 20 I = 1, N
         IF( ABS( D( I ) ).LT.EPS ) THEN
            D( I ) = SIGN( EPS, D( I ) )
         END IF
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      START = 1
      SQRE = 0
<span class="comment">*</span><span class="comment">
</span>      DO 30 I = 1, NM1
         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Subproblem found. First determine its size and then
</span><span class="comment">*</span><span class="comment">        apply divide and conquer on it.
</span><span class="comment">*</span><span class="comment">
</span>            IF( I.LT.NM1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A subproblem with E(I) small for I &lt; NM1.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - START + 1
            ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A subproblem with E(NM1) not too small but I = NM1.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = N - START + 1
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A subproblem with E(NM1) small. This implies an
</span><span class="comment">*</span><span class="comment">        1-by-1 subproblem at D(N). Solve this 1-by-1 problem
</span><span class="comment">*</span><span class="comment">        first.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - START + 1
               IF( ICOMPQ.EQ.2 ) THEN
                  U( N, N ) = SIGN( ONE, D( N ) )
                  VT( N, N ) = ONE
               ELSE IF( ICOMPQ.EQ.1 ) THEN
                  Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
                  Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
               END IF
               D( N ) = ABS( D( N ) )
            END IF
            IF( ICOMPQ.EQ.2 ) THEN
               CALL <a name="DLASD0.353"></a><a href="dlasd0.f.html#DLASD0.1">DLASD0</a>( NSIZE, SQRE, D( START ), E( START ),
     $                      U( START, START ), LDU, VT( START, START ),
     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
            ELSE
               CALL <a name="DLASDA.357"></a><a href="dlasda.f.html#DLASDA.1">DLASDA</a>( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
     $                      Q( START+( IVT+QSTART-2 )*N ),
     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
     $                      Q( START+( Z+QSTART-2 )*N ),
     $                      Q( START+( POLES+QSTART-2 )*N ),
     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
     $                      N, IQ( START+PERM*N ),
     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
     $                      Q( START+( IC+QSTART-2 )*N ),
     $                      Q( START+( IS+QSTART-2 )*N ),
     $                      WORK( WSTART ), IWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            START = I + 1
         END IF
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Unscale
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASCL.380"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Use Selection Sort to minimize swaps of singular vectors
</span><span class="comment">*</span><span class="comment">
</span>      DO 60 II = 2, N
         I = II - 1
         KK = I
         P = D( I )
         DO 50 J = II, N
            IF( D( J ).GT.P ) THEN
               KK = J
               P = D( J )
            END IF
   50    CONTINUE
         IF( KK.NE.I ) THEN
            D( KK ) = D( I )
            D( I ) = P
            IF( ICOMPQ.EQ.1 ) THEN
               IQ( I ) = KK
            ELSE IF( ICOMPQ.EQ.2 ) THEN
               CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
               CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
            END IF
         ELSE IF( ICOMPQ.EQ.1 ) THEN
            IQ( I ) = I
         END IF
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.1 ) THEN
         IF( IUPLO.EQ.1 ) THEN
            IQ( N ) = 1
         ELSE
            IQ( N ) = 0
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If B is lower bidiagonal, update U by those Givens rotations
</span><span class="comment">*</span><span class="comment">     which rotated B to be upper bidiagonal
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
     $   CALL <a name="DLASR.423"></a><a href="dlasr.f.html#DLASR.1">DLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, N, N, WORK( 1 ), WORK( N ), U, LDU )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DBDSDC.427"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>