zlalsa.f.html

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

            LL = 1
         ELSE
            LF = 2**( LVL-1 )
            LL = 2*LF - 1
         END IF
         DO 150 I = LF, LL
            IM1 = I - 1
            IC = IWORK( INODE+IM1 )
            NL = IWORK( NDIML+IM1 )
            NR = IWORK( NDIMR+IM1 )
            NLF = IC - NL
            NRF = IC + 1
            J = J - 1
            CALL <a name="ZLALS0.347"></a><a href="zlals0.f.html#ZLALS0.1">ZLALS0</a>( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
     $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
     $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
     $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
     $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
     $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK,
     $                   INFO )
  150    CONTINUE
  160 CONTINUE
      GO TO 330
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ICOMPQ = 1: applying back the right singular vector factors.
</span><span class="comment">*</span><span class="comment">
</span>  170 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     First now go through the right singular vector matrices of all
</span><span class="comment">*</span><span class="comment">     the tree nodes top-down.
</span><span class="comment">*</span><span class="comment">
</span>      J = 0
      DO 190 LVL = 1, NLVL
         LVL2 = 2*LVL - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Find the first node LF and last node LL on
</span><span class="comment">*</span><span class="comment">        the current level LVL.
</span><span class="comment">*</span><span class="comment">
</span>         IF( LVL.EQ.1 ) THEN
            LF = 1
            LL = 1
         ELSE
            LF = 2**( LVL-1 )
            LL = 2*LF - 1
         END IF
         DO 180 I = LL, LF, -1
            IM1 = I - 1
            IC = IWORK( INODE+IM1 )
            NL = IWORK( NDIML+IM1 )
            NR = IWORK( NDIMR+IM1 )
            NLF = IC - NL
            NRF = IC + 1
            IF( I.EQ.LL ) THEN
               SQRE = 0
            ELSE
               SQRE = 1
            END IF
            J = J + 1
            CALL <a name="ZLALS0.392"></a><a href="zlals0.f.html#ZLALS0.1">ZLALS0</a>( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
     $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
     $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
     $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
     $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
     $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK,
     $                   INFO )
  180    CONTINUE
  190 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The nodes on the bottom level of the tree were solved
</span><span class="comment">*</span><span class="comment">     by <a name="DLASDQ.403"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>. The corresponding right singular vector
</span><span class="comment">*</span><span class="comment">     matrices are in explicit form. Apply them back.
</span><span class="comment">*</span><span class="comment">
</span>      NDB1 = ( ND+1 ) / 2
      DO 320 I = NDB1, ND
         I1 = I - 1
         IC = IWORK( INODE+I1 )
         NL = IWORK( NDIML+I1 )
         NR = IWORK( NDIMR+I1 )
         NLP1 = NL + 1
         IF( I.EQ.ND ) THEN
            NRP1 = NR
         ELSE
            NRP1 = NR + 1
         END IF
         NLF = IC - NL
         NRF = IC + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Since B and BX are complex, the following call to DGEMM is
</span><span class="comment">*</span><span class="comment">        performed in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
</span><span class="comment">*</span><span class="comment">    $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
</span><span class="comment">*</span><span class="comment">
</span>         J = NLP1*NRHS*2
         DO 210 JCOL = 1, NRHS
            DO 200 JROW = NLF, NLF + NLP1 - 1
               J = J + 1
               RWORK( J ) = DBLE( B( JROW, JCOL ) )
  200       CONTINUE
  210    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
     $               RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO, RWORK( 1 ),
     $               NLP1 )
         J = NLP1*NRHS*2
         DO 230 JCOL = 1, NRHS
            DO 220 JROW = NLF, NLF + NLP1 - 1
               J = J + 1
               RWORK( J ) = DIMAG( B( JROW, JCOL ) )
  220       CONTINUE
  230    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
     $               RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO,
     $               RWORK( 1+NLP1*NRHS ), NLP1 )
         JREAL = 0
         JIMAG = NLP1*NRHS
         DO 250 JCOL = 1, NRHS
            DO 240 JROW = NLF, NLF + NLP1 - 1
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                            RWORK( JIMAG ) )
  240       CONTINUE
  250    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Since B and BX are complex, the following call to DGEMM is
</span><span class="comment">*</span><span class="comment">        performed in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
</span><span class="comment">*</span><span class="comment">    $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
</span><span class="comment">*</span><span class="comment">
</span>         J = NRP1*NRHS*2
         DO 270 JCOL = 1, NRHS
            DO 260 JROW = NRF, NRF + NRP1 - 1
               J = J + 1
               RWORK( J ) = DBLE( B( JROW, JCOL ) )
  260       CONTINUE
  270    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
     $               RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO, RWORK( 1 ),
     $               NRP1 )
         J = NRP1*NRHS*2
         DO 290 JCOL = 1, NRHS
            DO 280 JROW = NRF, NRF + NRP1 - 1
               J = J + 1
               RWORK( J ) = DIMAG( B( JROW, JCOL ) )
  280       CONTINUE
  290    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
     $               RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO,
     $               RWORK( 1+NRP1*NRHS ), NRP1 )
         JREAL = 0
         JIMAG = NRP1*NRHS
         DO 310 JCOL = 1, NRHS
            DO 300 JROW = NRF, NRF + NRP1 - 1
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                            RWORK( JIMAG ) )
  300       CONTINUE
  310    CONTINUE
<span class="comment">*</span><span class="comment">
</span>  320 CONTINUE
<span class="comment">*</span><span class="comment">
</span>  330 CONTINUE
<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="ZLALSA.501"></a><a href="zlalsa.f.html#ZLALSA.1">ZLALSA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>