clalsd.f.html

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

         DO 10 I = 1, N - 1
            CALL <a name="SLARTG.190"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I ), E( I ), CS, SN, R )
            D( I ) = R
            E( I ) = SN*D( I+1 )
            D( I+1 ) = CS*D( I+1 )
            IF( NRHS.EQ.1 ) THEN
               CALL CSROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
            ELSE
               RWORK( I*2-1 ) = CS
               RWORK( I*2 ) = SN
            END IF
   10    CONTINUE
         IF( NRHS.GT.1 ) THEN
            DO 30 I = 1, NRHS
               DO 20 J = 1, N - 1
                  CS = RWORK( J*2-1 )
                  SN = RWORK( J*2 )
                  CALL CSROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
   20          CONTINUE
   30       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale.
</span><span class="comment">*</span><span class="comment">
</span>      NM1 = N - 1
      ORGNRM = <a name="SLANST.215"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( ORGNRM.EQ.ZERO ) THEN
         CALL <a name="CLASET.217"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, N, NRHS, CZERO, CZERO, B, LDB )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASCL.221"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
      CALL <a name="SLASCL.222"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If N is smaller than the minimum divide size SMLSIZ, then solve
</span><span class="comment">*</span><span class="comment">     the problem with another solver.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LE.SMLSIZ ) THEN
         IRWU = 1
         IRWVT = IRWU + N*N
         IRWWRK = IRWVT + N*N
         IRWRB = IRWWRK
         IRWIB = IRWRB + N*NRHS
         IRWB = IRWIB + N*NRHS
         CALL <a name="SLASET.234"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, RWORK( IRWU ), N )
         CALL <a name="SLASET.235"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, RWORK( IRWVT ), N )
         CALL <a name="SLASDQ.236"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
     $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
     $                RWORK( IRWWRK ), INFO )
         IF( INFO.NE.0 ) THEN
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        In the real version, B is passed to <a name="SLASDQ.243"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a> and multiplied
</span><span class="comment">*</span><span class="comment">        internally by Q'. Here B is complex and that product is
</span><span class="comment">*</span><span class="comment">        computed below in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span>         J = IRWB - 1
         DO 50 JCOL = 1, NRHS
            DO 40 JROW = 1, N
               J = J + 1
               RWORK( J ) = REAL( B( JROW, JCOL ) )
   40       CONTINUE
   50    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWU ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
         J = IRWB - 1
         DO 70 JCOL = 1, NRHS
            DO 60 JROW = 1, N
               J = J + 1
               RWORK( J ) = AIMAG( B( JROW, JCOL ) )
   60       CONTINUE
   70    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWU ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
         JREAL = IRWRB - 1
         JIMAG = IRWIB - 1
         DO 90 JCOL = 1, NRHS
            DO 80 JROW = 1, N
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
   80       CONTINUE
   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
         DO 100 I = 1, N
            IF( D( I ).LE.TOL ) THEN
               CALL <a name="CLASET.278"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
            ELSE
               CALL <a name="CLASCL.280"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
     $                      LDB, INFO )
               RANK = RANK + 1
            END IF
  100    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Since B is complex, the following call to SGEMM is performed
</span><span class="comment">*</span><span class="comment">        in two steps (real and imaginary parts). That is for V * B
</span><span class="comment">*</span><span class="comment">        (in the real version of the code V' is stored in WORK).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
</span><span class="comment">*</span><span class="comment">    $               WORK( NWORK ), N )
</span><span class="comment">*</span><span class="comment">
</span>         J = IRWB - 1
         DO 120 JCOL = 1, NRHS
            DO 110 JROW = 1, N
               J = J + 1
               RWORK( J ) = REAL( B( JROW, JCOL ) )
  110       CONTINUE
  120    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWVT ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
         J = IRWB - 1
         DO 140 JCOL = 1, NRHS
            DO 130 JROW = 1, N
               J = J + 1
               RWORK( J ) = AIMAG( B( JROW, JCOL ) )
  130       CONTINUE
  140    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWVT ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
         JREAL = IRWRB - 1
         JIMAG = IRWIB - 1
         DO 160 JCOL = 1, NRHS
            DO 150 JROW = 1, N
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
  150       CONTINUE
  160    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="SLASCL.323"></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.324"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, D, INFO )
         CALL <a name="CLASCL.325"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Book-keeping and setting up some constants.
</span><span class="comment">*</span><span class="comment">
</span>      NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
<span class="comment">*</span><span class="comment">
</span>      SMLSZP = SMLSIZ + 1
<span class="comment">*</span><span class="comment">
</span>      U = 1
      VT = 1 + SMLSIZ*N
      DIFL = VT + SMLSZP*N
      DIFR = DIFL + NLVL*N
      Z = DIFR + NLVL*N*2
      C = Z + NLVL*N
      S = C + N
      POLES = S + N
      GIVNUM = POLES + 2*NLVL*N
      NRWORK = GIVNUM + 2*NLVL*N
      BX = 1
<span class="comment">*</span><span class="comment">
</span>      IRWRB = NRWORK
      IRWIB = IRWRB + SMLSIZ*NRHS
      IRWB = IRWIB + SMLSIZ*NRHS
<span class="comment">*</span><span class="comment">
</span>      SIZEI = 1 + N
      K = SIZEI + N
      GIVPTR = K + N
      PERM = GIVPTR + N
      GIVCOL = PERM + NLVL*N
      IWK = GIVCOL + NLVL*N*2
<span class="comment">*</span><span class="comment">
</span>      ST = 1
      SQRE = 0
      ICMPQ1 = 1
      ICMPQ2 = 0
      NSUB = 0
<span class="comment">*</span><span class="comment">
</span>      DO 170 I = 1, N
         IF( ABS( D( I ) ).LT.EPS ) THEN
            D( I ) = SIGN( EPS, D( I ) )
         END IF
  170 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 240 I = 1, NM1
         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.