zlaed0.f.html

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

      SUBPBS = 1
      TLVLS = 0
   10 CONTINUE
      IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
         DO 20 J = SUBPBS, 1, -1
            IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
            IWORK( 2*J-1 ) = IWORK( J ) / 2
   20    CONTINUE
         TLVLS = TLVLS + 1
         SUBPBS = 2*SUBPBS
         GO TO 10
      END IF
      DO 30 J = 2, SUBPBS
         IWORK( J ) = IWORK( J ) + IWORK( J-1 )
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
</span><span class="comment">*</span><span class="comment">     using rank-1 modifications (cuts).
</span><span class="comment">*</span><span class="comment">
</span>      SPM1 = SUBPBS - 1
      DO 40 I = 1, SPM1
         SUBMAT = IWORK( I ) + 1
         SMM1 = SUBMAT - 1
         D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
         D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      INDXQ = 4*N + 3
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set up workspaces for eigenvalues only/accumulate new vectors
</span><span class="comment">*</span><span class="comment">     routine
</span><span class="comment">*</span><span class="comment">
</span>      TEMP = LOG( DBLE( N ) ) / LOG( TWO )
      LGN = INT( TEMP )
      IF( 2**LGN.LT.N )
     $   LGN = LGN + 1
      IF( 2**LGN.LT.N )
     $   LGN = LGN + 1
      IPRMPT = INDXQ + N + 1
      IPERM = IPRMPT + N*LGN
      IQPTR = IPERM + N*LGN
      IGIVPT = IQPTR + N + 2
      IGIVCL = IGIVPT + N*LGN
<span class="comment">*</span><span class="comment">
</span>      IGIVNM = 1
      IQ = IGIVNM + 2*N*LGN
      IWREM = IQ + N**2 + 1
<span class="comment">*</span><span class="comment">     Initialize pointers
</span>      DO 50 I = 0, SUBPBS
         IWORK( IPRMPT+I ) = 1
         IWORK( IGIVPT+I ) = 1
   50 CONTINUE
      IWORK( IQPTR ) = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve each submatrix eigenproblem at the bottom of the divide and
</span><span class="comment">*</span><span class="comment">     conquer tree.
</span><span class="comment">*</span><span class="comment">
</span>      CURR = 0
      DO 70 I = 0, SPM1
         IF( I.EQ.0 ) THEN
            SUBMAT = 1
            MATSIZ = IWORK( 1 )
         ELSE
            SUBMAT = IWORK( I ) + 1
            MATSIZ = IWORK( I+1 ) - IWORK( I )
         END IF
         LL = IQ - 1 + IWORK( IQPTR+CURR )
         CALL <a name="DSTEQR.206"></a><a href="dsteqr.f.html#DSTEQR.1">DSTEQR</a>( <span class="string">'I'</span>, MATSIZ, D( SUBMAT ), E( SUBMAT ),
     $                RWORK( LL ), MATSIZ, RWORK, INFO )
         CALL <a name="ZLACRM.208"></a><a href="zlacrm.f.html#ZLACRM.1">ZLACRM</a>( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
     $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
     $                RWORK( IWREM ) )
         IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
         CURR = CURR + 1
         IF( INFO.GT.0 ) THEN
            INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
            RETURN
         END IF
         K = 1
         DO 60 J = SUBMAT, IWORK( I+1 )
            IWORK( INDXQ+J ) = K
            K = K + 1
   60    CONTINUE
   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Successively merge eigensystems of adjacent submatrices
</span><span class="comment">*</span><span class="comment">     into eigensystem for the corresponding larger matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     while ( SUBPBS &gt; 1 )
</span><span class="comment">*</span><span class="comment">
</span>      CURLVL = 1
   80 CONTINUE
      IF( SUBPBS.GT.1 ) THEN
         SPM2 = SUBPBS - 2
         DO 90 I = 0, SPM2, 2
            IF( I.EQ.0 ) THEN
               SUBMAT = 1
               MATSIZ = IWORK( 2 )
               MSD2 = IWORK( 1 )
               CURPRB = 0
            ELSE
               SUBMAT = IWORK( I ) + 1
               MATSIZ = IWORK( I+2 ) - IWORK( I )
               MSD2 = MATSIZ / 2
               CURPRB = CURPRB + 1
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
</span><span class="comment">*</span><span class="comment">     into an eigensystem of size MATSIZ.  <a name="ZLAED7.247"></a><a href="zlaed7.f.html#ZLAED7.1">ZLAED7</a> handles the case
</span><span class="comment">*</span><span class="comment">     when the eigenvectors of a full or band Hermitian matrix (which
</span><span class="comment">*</span><span class="comment">     was reduced to tridiagonal form) are desired.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     I am free to use Q as a valuable working space until Loop 150.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="ZLAED7.253"></a><a href="zlaed7.f.html#ZLAED7.1">ZLAED7</a>( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
     $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
     $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
     $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
     $                   IWORK( IPERM ), IWORK( IGIVPT ),
     $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
     $                   Q( 1, SUBMAT ), RWORK( IWREM ),
     $                   IWORK( SUBPBS+1 ), INFO )
            IF( INFO.GT.0 ) THEN
               INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
               RETURN
            END IF
            IWORK( I / 2+1 ) = IWORK( I+2 )
   90    CONTINUE
         SUBPBS = SUBPBS / 2
         CURLVL = CURLVL + 1
         GO TO 80
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     end while
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Re-merge the eigenvalues/vectors which were deflated at the final
</span><span class="comment">*</span><span class="comment">     merge step.
</span><span class="comment">*</span><span class="comment">
</span>      DO 100 I = 1, N
         J = IWORK( INDXQ+I )
         RWORK( I ) = D( J )
         CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
  100 CONTINUE
      CALL DCOPY( N, RWORK, 1, D, 1 )
<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="ZLAED0.286"></a><a href="zlaed0.f.html#ZLAED0.1">ZLAED0</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>