slalsd.f.html

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

         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
            NSUB = NSUB + 1
            IWORK( NSUB ) = ST
<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 - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
            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 - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
            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), which is not solved
</span><span class="comment">*</span><span class="comment">              explicitly.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
               NSUB = NSUB + 1
               IWORK( NSUB ) = N
               IWORK( SIZEI+NSUB-1 ) = 1
               CALL SCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
            END IF
            ST1 = ST - 1
            IF( NSIZE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              This is a 1-by-1 subproblem and is not solved
</span><span class="comment">*</span><span class="comment">              explicitly.
</span><span class="comment">*</span><span class="comment">
</span>               CALL SCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
            ELSE IF( NSIZE.LE.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              This is a small subproblem and is solved by <a name="SLASDQ.329"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASET.331"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, NSIZE, NSIZE, ZERO, ONE,
     $                      WORK( VT+ST1 ), N )
               CALL <a name="SLASDQ.333"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
     $                      E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
     $                      N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
               CALL <a name="SLACPY.339"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, NSIZE, NRHS, B( ST, 1 ), LDB,
     $                      WORK( BX+ST1 ), N )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A large problem. Solve it using divide and conquer.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASDA.345"></a><a href="slasda.f.html#SLASDA.1">SLASDA</a>( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
     $                      E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
     $                      IWORK( K+ST1 ), WORK( DIFL+ST1 ),
     $                      WORK( DIFR+ST1 ), WORK( Z+ST1 ),
     $                      WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
     $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
     $                      WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
     $                      WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
     $                      INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
               BXST = BX + ST1
               CALL <a name="SLALSA.358"></a><a href="slalsa.f.html#SLALSA.1">SLALSA</a>( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
     $                      LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
     $                      WORK( VT+ST1 ), IWORK( K+ST1 ),
     $                      WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
     $                      WORK( Z+ST1 ), WORK( POLES+ST1 ),
     $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
     $                      IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
     $                      WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
     $                      IWORK( IWK ), INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            ST = I + 1
         END IF
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply the singular values and treat the tiny ones as zero.
</span><span class="comment">*</span><span class="comment">
</span>      TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      DO 70 I = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Some of the elements in D can be negative because 1-by-1
</span><span class="comment">*</span><span class="comment">        subproblems were not solved explicitly.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( D( I ) ).LE.TOL ) THEN
            CALL <a name="SLASET.385"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
         ELSE
            RANK = RANK + 1
            CALL <a name="SLASCL.388"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, D( I ), ONE, 1, NRHS,
     $                   WORK( BX+I-1 ), N, INFO )
         END IF
         D( I ) = ABS( D( I ) )
   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Now apply back the right singular vectors.
</span><span class="comment">*</span><span class="comment">
</span>      ICMPQ2 = 1
      DO 80 I = 1, NSUB
         ST = IWORK( I )
         ST1 = ST - 1
         NSIZE = IWORK( SIZEI+I-1 )
         BXST = BX + ST1
         IF( NSIZE.EQ.1 ) THEN
            CALL SCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
         ELSE IF( NSIZE.LE.SMLSIZ ) THEN
            CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                  WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
     $                  B( ST, 1 ), LDB )
         ELSE
            CALL <a name="SLALSA.409"></a><a href="slalsa.f.html#SLALSA.1">SLALSA</a>( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
     $                   B( ST, 1 ), LDB, WORK( U+ST1 ), N,
     $                   WORK( VT+ST1 ), IWORK( K+ST1 ),
     $                   WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
     $                   WORK( Z+ST1 ), WORK( POLES+ST1 ),
     $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
     $                   IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
     $                   WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
     $                   IWORK( IWK ), INFO )
            IF( INFO.NE.0 ) THEN
               RETURN
            END IF
         END IF
   80 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Unscale and sort the singular values.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASCL.426"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
      CALL <a name="SLASRT.427"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, D, INFO )
      CALL <a name="SLASCL.428"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
<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="SLALSD.432"></a><a href="slalsd.f.html#SLALSD.1">SLALSD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>