LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dggsvp.f
Go to the documentation of this file.
00001 *> \brief \b DGGSVP
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DGGSVP + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvp.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvp.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvp.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
00022 *                          TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
00023 *                          IWORK, TAU, WORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          JOBQ, JOBU, JOBV
00027 *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
00028 *       DOUBLE PRECISION   TOLA, TOLB
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       INTEGER            IWORK( * )
00032 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00033 *      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> DGGSVP computes orthogonal matrices U, V and Q such that
00043 *>
00044 *>                    N-K-L  K    L
00045 *>  U**T*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
00046 *>                 L ( 0     0   A23 )
00047 *>             M-K-L ( 0     0    0  )
00048 *>
00049 *>                  N-K-L  K    L
00050 *>         =     K ( 0    A12  A13 )  if M-K-L < 0;
00051 *>             M-K ( 0     0   A23 )
00052 *>
00053 *>                  N-K-L  K    L
00054 *>  V**T*B*Q =   L ( 0     0   B13 )
00055 *>             P-L ( 0     0    0  )
00056 *>
00057 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
00058 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
00059 *> otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
00060 *> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T. 
00061 *>
00062 *> This decomposition is the preprocessing step for computing the
00063 *> Generalized Singular Value Decomposition (GSVD), see subroutine
00064 *> DGGSVD.
00065 *> \endverbatim
00066 *
00067 *  Arguments:
00068 *  ==========
00069 *
00070 *> \param[in] JOBU
00071 *> \verbatim
00072 *>          JOBU is CHARACTER*1
00073 *>          = 'U':  Orthogonal matrix U is computed;
00074 *>          = 'N':  U is not computed.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] JOBV
00078 *> \verbatim
00079 *>          JOBV is CHARACTER*1
00080 *>          = 'V':  Orthogonal matrix V is computed;
00081 *>          = 'N':  V is not computed.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] JOBQ
00085 *> \verbatim
00086 *>          JOBQ is CHARACTER*1
00087 *>          = 'Q':  Orthogonal matrix Q is computed;
00088 *>          = 'N':  Q is not computed.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] M
00092 *> \verbatim
00093 *>          M is INTEGER
00094 *>          The number of rows of the matrix A.  M >= 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] P
00098 *> \verbatim
00099 *>          P is INTEGER
00100 *>          The number of rows of the matrix B.  P >= 0.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] N
00104 *> \verbatim
00105 *>          N is INTEGER
00106 *>          The number of columns of the matrices A and B.  N >= 0.
00107 *> \endverbatim
00108 *>
00109 *> \param[in,out] A
00110 *> \verbatim
00111 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00112 *>          On entry, the M-by-N matrix A.
00113 *>          On exit, A contains the triangular (or trapezoidal) matrix
00114 *>          described in the Purpose section.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] LDA
00118 *> \verbatim
00119 *>          LDA is INTEGER
00120 *>          The leading dimension of the array A. LDA >= max(1,M).
00121 *> \endverbatim
00122 *>
00123 *> \param[in,out] B
00124 *> \verbatim
00125 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
00126 *>          On entry, the P-by-N matrix B.
00127 *>          On exit, B contains the triangular matrix described in
00128 *>          the Purpose section.
00129 *> \endverbatim
00130 *>
00131 *> \param[in] LDB
00132 *> \verbatim
00133 *>          LDB is INTEGER
00134 *>          The leading dimension of the array B. LDB >= max(1,P).
00135 *> \endverbatim
00136 *>
00137 *> \param[in] TOLA
00138 *> \verbatim
00139 *>          TOLA is DOUBLE PRECISION
00140 *> \endverbatim
00141 *>
00142 *> \param[in] TOLB
00143 *> \verbatim
00144 *>          TOLB is DOUBLE PRECISION
00145 *>
00146 *>          TOLA and TOLB are the thresholds to determine the effective
00147 *>          numerical rank of matrix B and a subblock of A. Generally,
00148 *>          they are set to
00149 *>             TOLA = MAX(M,N)*norm(A)*MACHEPS,
00150 *>             TOLB = MAX(P,N)*norm(B)*MACHEPS.
00151 *>          The size of TOLA and TOLB may affect the size of backward
00152 *>          errors of the decomposition.
00153 *> \endverbatim
00154 *>
00155 *> \param[out] K
00156 *> \verbatim
00157 *>          K is INTEGER
00158 *> \endverbatim
00159 *>
00160 *> \param[out] L
00161 *> \verbatim
00162 *>          L is INTEGER
00163 *>
00164 *>          On exit, K and L specify the dimension of the subblocks
00165 *>          described in Purpose section.
00166 *>          K + L = effective numerical rank of (A**T,B**T)**T.
00167 *> \endverbatim
00168 *>
00169 *> \param[out] U
00170 *> \verbatim
00171 *>          U is DOUBLE PRECISION array, dimension (LDU,M)
00172 *>          If JOBU = 'U', U contains the orthogonal matrix U.
00173 *>          If JOBU = 'N', U is not referenced.
00174 *> \endverbatim
00175 *>
00176 *> \param[in] LDU
00177 *> \verbatim
00178 *>          LDU is INTEGER
00179 *>          The leading dimension of the array U. LDU >= max(1,M) if
00180 *>          JOBU = 'U'; LDU >= 1 otherwise.
00181 *> \endverbatim
00182 *>
00183 *> \param[out] V
00184 *> \verbatim
00185 *>          V is DOUBLE PRECISION array, dimension (LDV,P)
00186 *>          If JOBV = 'V', V contains the orthogonal matrix V.
00187 *>          If JOBV = 'N', V is not referenced.
00188 *> \endverbatim
00189 *>
00190 *> \param[in] LDV
00191 *> \verbatim
00192 *>          LDV is INTEGER
00193 *>          The leading dimension of the array V. LDV >= max(1,P) if
00194 *>          JOBV = 'V'; LDV >= 1 otherwise.
00195 *> \endverbatim
00196 *>
00197 *> \param[out] Q
00198 *> \verbatim
00199 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
00200 *>          If JOBQ = 'Q', Q contains the orthogonal matrix Q.
00201 *>          If JOBQ = 'N', Q is not referenced.
00202 *> \endverbatim
00203 *>
00204 *> \param[in] LDQ
00205 *> \verbatim
00206 *>          LDQ is INTEGER
00207 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
00208 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
00209 *> \endverbatim
00210 *>
00211 *> \param[out] IWORK
00212 *> \verbatim
00213 *>          IWORK is INTEGER array, dimension (N)
00214 *> \endverbatim
00215 *>
00216 *> \param[out] TAU
00217 *> \verbatim
00218 *>          TAU is DOUBLE PRECISION array, dimension (N)
00219 *> \endverbatim
00220 *>
00221 *> \param[out] WORK
00222 *> \verbatim
00223 *>          WORK is DOUBLE PRECISION array, dimension (max(3*N,M,P))
00224 *> \endverbatim
00225 *>
00226 *> \param[out] INFO
00227 *> \verbatim
00228 *>          INFO is INTEGER
00229 *>          = 0:  successful exit
00230 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00231 *> \endverbatim
00232 *
00233 *  Authors:
00234 *  ========
00235 *
00236 *> \author Univ. of Tennessee 
00237 *> \author Univ. of California Berkeley 
00238 *> \author Univ. of Colorado Denver 
00239 *> \author NAG Ltd. 
00240 *
00241 *> \date November 2011
00242 *
00243 *> \ingroup doubleOTHERcomputational
00244 *
00245 *> \par Further Details:
00246 *  =====================
00247 *>
00248 *>  The subroutine uses LAPACK subroutine DGEQPF for the QR factorization
00249 *>  with column pivoting to detect the effective numerical rank of the
00250 *>  a matrix. It may be replaced by a better rank determination strategy.
00251 *>
00252 *  =====================================================================
00253       SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
00254      $                   TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
00255      $                   IWORK, TAU, WORK, INFO )
00256 *
00257 *  -- LAPACK computational routine (version 3.4.0) --
00258 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00259 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00260 *     November 2011
00261 *
00262 *     .. Scalar Arguments ..
00263       CHARACTER          JOBQ, JOBU, JOBV
00264       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
00265       DOUBLE PRECISION   TOLA, TOLB
00266 *     ..
00267 *     .. Array Arguments ..
00268       INTEGER            IWORK( * )
00269       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00270      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
00271 *     ..
00272 *
00273 *  =====================================================================
00274 *
00275 *     .. Parameters ..
00276       DOUBLE PRECISION   ZERO, ONE
00277       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00278 *     ..
00279 *     .. Local Scalars ..
00280       LOGICAL            FORWRD, WANTQ, WANTU, WANTV
00281       INTEGER            I, J
00282 *     ..
00283 *     .. External Functions ..
00284       LOGICAL            LSAME
00285       EXTERNAL           LSAME
00286 *     ..
00287 *     .. External Subroutines ..
00288       EXTERNAL           DGEQPF, DGEQR2, DGERQ2, DLACPY, DLAPMT, DLASET,
00289      $                   DORG2R, DORM2R, DORMR2, XERBLA
00290 *     ..
00291 *     .. Intrinsic Functions ..
00292       INTRINSIC          ABS, MAX, MIN
00293 *     ..
00294 *     .. Executable Statements ..
00295 *
00296 *     Test the input parameters
00297 *
00298       WANTU = LSAME( JOBU, 'U' )
00299       WANTV = LSAME( JOBV, 'V' )
00300       WANTQ = LSAME( JOBQ, 'Q' )
00301       FORWRD = .TRUE.
00302 *
00303       INFO = 0
00304       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
00305          INFO = -1
00306       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00307          INFO = -2
00308       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
00309          INFO = -3
00310       ELSE IF( M.LT.0 ) THEN
00311          INFO = -4
00312       ELSE IF( P.LT.0 ) THEN
00313          INFO = -5
00314       ELSE IF( N.LT.0 ) THEN
00315          INFO = -6
00316       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00317          INFO = -8
00318       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
00319          INFO = -10
00320       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
00321          INFO = -16
00322       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
00323          INFO = -18
00324       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00325          INFO = -20
00326       END IF
00327       IF( INFO.NE.0 ) THEN
00328          CALL XERBLA( 'DGGSVP', -INFO )
00329          RETURN
00330       END IF
00331 *
00332 *     QR with column pivoting of B: B*P = V*( S11 S12 )
00333 *                                           (  0   0  )
00334 *
00335       DO 10 I = 1, N
00336          IWORK( I ) = 0
00337    10 CONTINUE
00338       CALL DGEQPF( P, N, B, LDB, IWORK, TAU, WORK, INFO )
00339 *
00340 *     Update A := A*P
00341 *
00342       CALL DLAPMT( FORWRD, M, N, A, LDA, IWORK )
00343 *
00344 *     Determine the effective rank of matrix B.
00345 *
00346       L = 0
00347       DO 20 I = 1, MIN( P, N )
00348          IF( ABS( B( I, I ) ).GT.TOLB )
00349      $      L = L + 1
00350    20 CONTINUE
00351 *
00352       IF( WANTV ) THEN
00353 *
00354 *        Copy the details of V, and form V.
00355 *
00356          CALL DLASET( 'Full', P, P, ZERO, ZERO, V, LDV )
00357          IF( P.GT.1 )
00358      $      CALL DLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
00359      $                   LDV )
00360          CALL DORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
00361       END IF
00362 *
00363 *     Clean up B
00364 *
00365       DO 40 J = 1, L - 1
00366          DO 30 I = J + 1, L
00367             B( I, J ) = ZERO
00368    30    CONTINUE
00369    40 CONTINUE
00370       IF( P.GT.L )
00371      $   CALL DLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB )
00372 *
00373       IF( WANTQ ) THEN
00374 *
00375 *        Set Q = I and Update Q := Q*P
00376 *
00377          CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00378          CALL DLAPMT( FORWRD, N, N, Q, LDQ, IWORK )
00379       END IF
00380 *
00381       IF( P.GE.L .AND. N.NE.L ) THEN
00382 *
00383 *        RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
00384 *
00385          CALL DGERQ2( L, N, B, LDB, TAU, WORK, INFO )
00386 *
00387 *        Update A := A*Z**T
00388 *
00389          CALL DORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A,
00390      $                LDA, WORK, INFO )
00391 *
00392          IF( WANTQ ) THEN
00393 *
00394 *           Update Q := Q*Z**T
00395 *
00396             CALL DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,
00397      $                   LDQ, WORK, INFO )
00398          END IF
00399 *
00400 *        Clean up B
00401 *
00402          CALL DLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB )
00403          DO 60 J = N - L + 1, N
00404             DO 50 I = J - N + L + 1, L
00405                B( I, J ) = ZERO
00406    50       CONTINUE
00407    60    CONTINUE
00408 *
00409       END IF
00410 *
00411 *     Let              N-L     L
00412 *                A = ( A11    A12 ) M,
00413 *
00414 *     then the following does the complete QR decomposition of A11:
00415 *
00416 *              A11 = U*(  0  T12 )*P1**T
00417 *                      (  0   0  )
00418 *
00419       DO 70 I = 1, N - L
00420          IWORK( I ) = 0
00421    70 CONTINUE
00422       CALL DGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, INFO )
00423 *
00424 *     Determine the effective rank of A11
00425 *
00426       K = 0
00427       DO 80 I = 1, MIN( M, N-L )
00428          IF( ABS( A( I, I ) ).GT.TOLA )
00429      $      K = K + 1
00430    80 CONTINUE
00431 *
00432 *     Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
00433 *
00434       CALL DORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA,
00435      $             TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
00436 *
00437       IF( WANTU ) THEN
00438 *
00439 *        Copy the details of U, and form U
00440 *
00441          CALL DLASET( 'Full', M, M, ZERO, ZERO, U, LDU )
00442          IF( M.GT.1 )
00443      $      CALL DLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
00444      $                   LDU )
00445          CALL DORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO )
00446       END IF
00447 *
00448       IF( WANTQ ) THEN
00449 *
00450 *        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1
00451 *
00452          CALL DLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK )
00453       END IF
00454 *
00455 *     Clean up A: set the strictly lower triangular part of
00456 *     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
00457 *
00458       DO 100 J = 1, K - 1
00459          DO 90 I = J + 1, K
00460             A( I, J ) = ZERO
00461    90    CONTINUE
00462   100 CONTINUE
00463       IF( M.GT.K )
00464      $   CALL DLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA )
00465 *
00466       IF( N-L.GT.K ) THEN
00467 *
00468 *        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
00469 *
00470          CALL DGERQ2( K, N-L, A, LDA, TAU, WORK, INFO )
00471 *
00472          IF( WANTQ ) THEN
00473 *
00474 *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
00475 *
00476             CALL DORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU,
00477      $                   Q, LDQ, WORK, INFO )
00478          END IF
00479 *
00480 *        Clean up A
00481 *
00482          CALL DLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA )
00483          DO 120 J = N - L - K + 1, N - L
00484             DO 110 I = J - N + L + K + 1, K
00485                A( I, J ) = ZERO
00486   110       CONTINUE
00487   120    CONTINUE
00488 *
00489       END IF
00490 *
00491       IF( M.GT.K ) THEN
00492 *
00493 *        QR factorization of A( K+1:M,N-L+1:N )
00494 *
00495          CALL DGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO )
00496 *
00497          IF( WANTU ) THEN
00498 *
00499 *           Update U(:,K+1:M) := U(:,K+1:M)*U1
00500 *
00501             CALL DORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ),
00502      $                   A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU,
00503      $                   WORK, INFO )
00504          END IF
00505 *
00506 *        Clean up
00507 *
00508          DO 140 J = N - L + 1, N
00509             DO 130 I = J - N + K + L + 1, M
00510                A( I, J ) = ZERO
00511   130       CONTINUE
00512   140    CONTINUE
00513 *
00514       END IF
00515 *
00516       RETURN
00517 *
00518 *     End of DGGSVP
00519 *
00520       END
 All Files Functions