dggsvp.f.html

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

</span><span class="comment">*</span><span class="comment">          If JOBQ = 'N', Q is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDQ     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array Q. LDQ &gt;= max(1,N) if
</span><span class="comment">*</span><span class="comment">          JOBQ = 'Q'; LDQ &gt;= 1 otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TAU     (workspace) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P))
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The subroutine uses LAPACK subroutine <a name="DGEQPF.142"></a><a href="dgeqpf.f.html#DGEQPF.1">DGEQPF</a> for the QR factorization
</span><span class="comment">*</span><span class="comment">  with column pivoting to detect the effective numerical rank of the
</span><span class="comment">*</span><span class="comment">  a matrix. It may be replaced by a better rank determination strategy.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            FORWRD, WANTQ, WANTU, WANTV
      INTEGER            I, J
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.157"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      EXTERNAL           <a name="LSAME.158"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DGEQPF.161"></a><a href="dgeqpf.f.html#DGEQPF.1">DGEQPF</a>, <a name="DGEQR2.161"></a><a href="dgeqr2.f.html#DGEQR2.1">DGEQR2</a>, <a name="DGERQ2.161"></a><a href="dgerq2.f.html#DGERQ2.1">DGERQ2</a>, <a name="DLACPY.161"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DLAPMT.161"></a><a href="dlapmt.f.html#DLAPMT.1">DLAPMT</a>, <a name="DLASET.161"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>,
     $                   <a name="DORG2R.162"></a><a href="dorg2r.f.html#DORG2R.1">DORG2R</a>, <a name="DORM2R.162"></a><a href="dorm2r.f.html#DORM2R.1">DORM2R</a>, <a name="DORMR2.162"></a><a href="dormr2.f.html#DORMR2.1">DORMR2</a>, <a name="XERBLA.162"></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          ABS, MAX, MIN
<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>      WANTU = <a name="LSAME.171"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBU, <span class="string">'U'</span> )
      WANTV = <a name="LSAME.172"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBV, <span class="string">'V'</span> )
      WANTQ = <a name="LSAME.173"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBQ, <span class="string">'Q'</span> )
      FORWRD = .TRUE.
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.( WANTU .OR. <a name="LSAME.177"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBU, <span class="string">'N'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( WANTV .OR. <a name="LSAME.179"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBV, <span class="string">'N'</span> ) ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( WANTQ .OR. <a name="LSAME.181"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBQ, <span class="string">'N'</span> ) ) ) THEN
         INFO = -3
      ELSE IF( M.LT.0 ) THEN
         INFO = -4
      ELSE IF( P.LT.0 ) THEN
         INFO = -5
      ELSE IF( N.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -8
      ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
         INFO = -10
      ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
         INFO = -16
      ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
         INFO = -18
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -20
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.201"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DGGSVP.201"></a><a href="dggsvp.f.html#DGGSVP.1">DGGSVP</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     QR with column pivoting of B: B*P = V*( S11 S12 )
</span><span class="comment">*</span><span class="comment">                                           (  0   0  )
</span><span class="comment">*</span><span class="comment">
</span>      DO 10 I = 1, N
         IWORK( I ) = 0
   10 CONTINUE
      CALL <a name="DGEQPF.211"></a><a href="dgeqpf.f.html#DGEQPF.1">DGEQPF</a>( P, N, B, LDB, IWORK, TAU, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update A := A*P
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLAPMT.215"></a><a href="dlapmt.f.html#DLAPMT.1">DLAPMT</a>( FORWRD, M, N, A, LDA, IWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine the effective rank of matrix B.
</span><span class="comment">*</span><span class="comment">
</span>      L = 0
      DO 20 I = 1, MIN( P, N )
         IF( ABS( B( I, I ) ).GT.TOLB )
     $      L = L + 1
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IF( WANTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy the details of V, and form V.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASET.229"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, P, P, ZERO, ZERO, V, LDV )
         IF( P.GT.1 )
     $      CALL <a name="DLACPY.231"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Lower'</span>, P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
     $                   LDV )
         CALL <a name="DORG2R.233"></a><a href="dorg2r.f.html#DORG2R.1">DORG2R</a>( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Clean up B
</span><span class="comment">*</span><span class="comment">
</span>      DO 40 J = 1, L - 1
         DO 30 I = J + 1, L
            B( I, J ) = ZERO
   30    CONTINUE
   40 CONTINUE
      IF( P.GT.L )
     $   CALL <a name="DLASET.244"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTQ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set Q = I and Update Q := Q*P
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASET.250"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, Q, LDQ )
         CALL <a name="DLAPMT.251"></a><a href="dlapmt.f.html#DLAPMT.1">DLAPMT</a>( FORWRD, N, N, Q, LDQ, IWORK )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( P.GE.L .AND. N.NE.L ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGERQ2.258"></a><a href="dgerq2.f.html#DGERQ2.1">DGERQ2</a>( L, N, B, LDB, TAU, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update A := A*Z'
</span><span class="comment">*</span><span class="comment">