dlaeda.f.html

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

</span>      EXTERNAL           DCOPY, DGEMV, DROT, <a name="XERBLA.103"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          DBLE, INT, SQRT
<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( N.LT.0 ) THEN
         INFO = -1
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.118"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLAEDA.118"></a><a href="dlaeda.f.html#DLAEDA.1">DLAEDA</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine location of first number in second half.
</span><span class="comment">*</span><span class="comment">
</span>      MID = N / 2 + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Gather last/first rows of appropriate eigenblocks into center of Z
</span><span class="comment">*</span><span class="comment">
</span>      PTR = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine location of lowest level subproblem in the full storage
</span><span class="comment">*</span><span class="comment">     scheme
</span><span class="comment">*</span><span class="comment">
</span>      CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine size of these matrices.  We add HALF to the value of
</span><span class="comment">*</span><span class="comment">     the SQRT in case the machine underestimates one of these square
</span><span class="comment">*</span><span class="comment">     roots.
</span><span class="comment">*</span><span class="comment">
</span>      BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
      BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
      DO 10 K = 1, MID - BSIZ1 - 1
         Z( K ) = ZERO
   10 CONTINUE
      CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
     $            Z( MID-BSIZ1 ), 1 )
      CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
      DO 20 K = MID + BSIZ2, N
         Z( K ) = ZERO
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Loop thru remaining levels 1 -&gt; CURLVL applying the Givens
</span><span class="comment">*</span><span class="comment">     rotations and permutation and then multiplying the center matrices
</span><span class="comment">*</span><span class="comment">     against the current Z.
</span><span class="comment">*</span><span class="comment">
</span>      PTR = 2**TLVLS + 1
      DO 70 K = 1, CURLVL - 1
         CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
         PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
         PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
         ZPTR1 = MID - PSIZ1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       Apply Givens at CURR and CURR+1
</span><span class="comment">*</span><span class="comment">
</span>         DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
            CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
     $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
     $                 GIVNUM( 2, I ) )
   30    CONTINUE
         DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
            CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
     $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
     $                 GIVNUM( 2, I ) )
   40    CONTINUE
         PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
         PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
         DO 50 I = 0, PSIZ1 - 1
            ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
   50    CONTINUE
         DO 60 I = 0, PSIZ2 - 1
            ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply Blocks at CURR and CURR+1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine size of these matrices.  We add HALF to the value of
</span><span class="comment">*</span><span class="comment">        the SQRT in case the machine underestimates one of these
</span><span class="comment">*</span><span class="comment">        square roots.
</span><span class="comment">*</span><span class="comment">
</span>         BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
         BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
     $           1 ) ) ) )
         IF( BSIZ1.GT.0 ) THEN
            CALL DGEMV( <span class="string">'T'</span>, BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
     $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
         END IF
         CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
     $               1 )
         IF( BSIZ2.GT.0 ) THEN
            CALL DGEMV( <span class="string">'T'</span>, BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
     $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
         END IF
         CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
     $               Z( MID+BSIZ2 ), 1 )
<span class="comment">*</span><span class="comment">
</span>         PTR = PTR + 2**( TLVLS-K )
   70 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="DLAEDA.215"></a><a href="dlaeda.f.html#DLAEDA.1">DLAEDA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>