zlalsa.f.html

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

</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
     $                   JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML,
     $                   NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DGEMM, <a name="DLASDT.168"></a><a href="dlasdt.f.html#DLASDT.1">DLASDT</a>, <a name="XERBLA.168"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZCOPY, <a name="ZLALS0.168"></a><a href="zlals0.f.html#ZLALS0.1">ZLALS0</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          DBLE, DCMPLX, DIMAG
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
         INFO = -1
      ELSE IF( SMLSIZ.LT.3 ) THEN
         INFO = -2
      ELSE IF( N.LT.SMLSIZ ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.1 ) THEN
         INFO = -4
      ELSE IF( LDB.LT.N ) THEN
         INFO = -6
      ELSE IF( LDBX.LT.N ) THEN
         INFO = -8
      ELSE IF( LDU.LT.N ) THEN
         INFO = -10
      ELSE IF( LDGCOL.LT.N ) THEN
         INFO = -19
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.197"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZLALSA.197"></a><a href="zlalsa.f.html#ZLALSA.1">ZLALSA</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Book-keeping and  setting up the computation tree.
</span><span class="comment">*</span><span class="comment">
</span>      INODE = 1
      NDIML = INODE + N
      NDIMR = NDIML + N
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASDT.207"></a><a href="dlasdt.f.html#DLASDT.1">DLASDT</a>( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
     $             IWORK( NDIMR ), SMLSIZ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The following code applies back the left singular vector factors.
</span><span class="comment">*</span><span class="comment">     For applying back the right singular vector factors, go to 170.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.1 ) THEN
         GO TO 170
      END IF
<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.218"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>. The corresponding left and right singular vector
</span><span class="comment">*</span><span class="comment">     matrices are in explicit form. First apply back the left
</span><span class="comment">*</span><span class="comment">     singular vector matrices.
</span><span class="comment">*</span><span class="comment">
</span>      NDB1 = ( ND+1 ) / 2
      DO 130 I = NDB1, ND
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        IC : center row of each node
</span><span class="comment">*</span><span class="comment">        NL : number of rows of left  subproblem
</span><span class="comment">*</span><span class="comment">        NR : number of rows of right subproblem
</span><span class="comment">*</span><span class="comment">        NLF: starting row of the left   subproblem
</span><span class="comment">*</span><span class="comment">        NRF: starting row of the right  subproblem
</span><span class="comment">*</span><span class="comment">
</span>         I1 = I - 1
         IC = IWORK( INODE+I1 )
         NL = IWORK( NDIML+I1 )
         NR = IWORK( NDIMR+I1 )
         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
</span><span class="comment">*</span><span class="comment">        is 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', NL, NRHS, NL, ONE, U( 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 = NL*NRHS*2
         DO 20 JCOL = 1, NRHS
            DO 10 JROW = NLF, NLF + NL - 1
               J = J + 1
               RWORK( J ) = DBLE( B( JROW, JCOL ) )
   10       CONTINUE
   20    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
     $               RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1 ), NL )
         J = NL*NRHS*2
         DO 40 JCOL = 1, NRHS
            DO 30 JROW = NLF, NLF + NL - 1
               J = J + 1
               RWORK( J ) = DIMAG( B( JROW, JCOL ) )
   30       CONTINUE
   40    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
     $               RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1+NL*NRHS ),
     $               NL )
         JREAL = 0
         JIMAG = NL*NRHS
         DO 60 JCOL = 1, NRHS
            DO 50 JROW = NLF, NLF + NL - 1
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                            RWORK( JIMAG ) )
   50       CONTINUE
   60    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
</span><span class="comment">*</span><span class="comment">        is 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', NR, NRHS, NR, ONE, U( 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 = NR*NRHS*2
         DO 80 JCOL = 1, NRHS
            DO 70 JROW = NRF, NRF + NR - 1
               J = J + 1
               RWORK( J ) = DBLE( B( JROW, JCOL ) )
   70       CONTINUE
   80    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
     $               RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1 ), NR )
         J = NR*NRHS*2
         DO 100 JCOL = 1, NRHS
            DO 90 JROW = NRF, NRF + NR - 1
               J = J + 1
               RWORK( J ) = DIMAG( B( JROW, JCOL ) )
   90       CONTINUE
  100    CONTINUE
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
     $               RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1+NR*NRHS ),
     $               NR )
         JREAL = 0
         JIMAG = NR*NRHS
         DO 120 JCOL = 1, NRHS
            DO 110 JROW = NRF, NRF + NR - 1
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                            RWORK( JIMAG ) )
  110       CONTINUE
  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span>  130 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Next copy the rows of B that correspond to unchanged rows
</span><span class="comment">*</span><span class="comment">     in the bidiagonal matrix to BX.
</span><span class="comment">*</span><span class="comment">
</span>      DO 140 I = 1, ND
         IC = IWORK( INODE+I-1 )
         CALL ZCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
  140 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Finally go through the left singular vector matrices of all
</span><span class="comment">*</span><span class="comment">     the other subproblems bottom-up on the tree.
</span><span class="comment">*</span><span class="comment">
</span>      J = 2**NLVL
      SQRE = 0
<span class="comment">*</span><span class="comment">
</span>      DO 160 LVL = NLVL, 1, -1
         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