SRC/DEPRECATED/dggsvp.f

   1 *> \brief \b DGGSVP
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGGSVP + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvp.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvp.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvp.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
  22 *                          TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
  23 *                          IWORK, TAU, WORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBQ, JOBU, JOBV
  27 *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
  28 *       DOUBLE PRECISION   TOLA, TOLB
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IWORK( * )
  32 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  33 *      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> This routine is deprecated and has been replaced by routine DGGSVP3.
  43 *>
  44 *> DGGSVP computes orthogonal matrices U, V and Q such that
  45 *>
  46 *>                    N-K-L  K    L
  47 *>  U**T*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
  48 *>                 L ( 0     0   A23 )
  49 *>             M-K-L ( 0     0    0  )
  50 *>
  51 *>                  N-K-L  K    L
  52 *>         =     K ( 0    A12  A13 )  if M-K-L < 0;
  53 *>             M-K ( 0     0   A23 )
  54 *>
  55 *>                  N-K-L  K    L
  56 *>  V**T*B*Q =   L ( 0     0   B13 )
  57 *>             P-L ( 0     0    0  )
  58 *>
  59 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
  60 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
  61 *> otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
  62 *> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T.
  63 *>
  64 *> This decomposition is the preprocessing step for computing the
  65 *> Generalized Singular Value Decomposition (GSVD), see subroutine
  66 *> DGGSVD.
  67 *> \endverbatim
  68 *
  69 *  Arguments:
  70 *  ==========
  71 *
  72 *> \param[in] JOBU
  73 *> \verbatim
  74 *>          JOBU is CHARACTER*1
  75 *>          = 'U':  Orthogonal matrix U is computed;
  76 *>          = 'N':  U is not computed.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] JOBV
  80 *> \verbatim
  81 *>          JOBV is CHARACTER*1
  82 *>          = 'V':  Orthogonal matrix V is computed;
  83 *>          = 'N':  V is not computed.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] JOBQ
  87 *> \verbatim
  88 *>          JOBQ is CHARACTER*1
  89 *>          = 'Q':  Orthogonal matrix Q is computed;
  90 *>          = 'N':  Q is not computed.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] M
  94 *> \verbatim
  95 *>          M is INTEGER
  96 *>          The number of rows of the matrix A.  M >= 0.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] P
 100 *> \verbatim
 101 *>          P is INTEGER
 102 *>          The number of rows of the matrix B.  P >= 0.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] N
 106 *> \verbatim
 107 *>          N is INTEGER
 108 *>          The number of columns of the matrices A and B.  N >= 0.
 109 *> \endverbatim
 110 *>
 111 *> \param[in,out] A
 112 *> \verbatim
 113 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 114 *>          On entry, the M-by-N matrix A.
 115 *>          On exit, A contains the triangular (or trapezoidal) matrix
 116 *>          described in the Purpose section.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] LDA
 120 *> \verbatim
 121 *>          LDA is INTEGER
 122 *>          The leading dimension of the array A. LDA >= max(1,M).
 123 *> \endverbatim
 124 *>
 125 *> \param[in,out] B
 126 *> \verbatim
 127 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
 128 *>          On entry, the P-by-N matrix B.
 129 *>          On exit, B contains the triangular matrix described in
 130 *>          the Purpose section.
 131 *> \endverbatim
 132 *>
 133 *> \param[in] LDB
 134 *> \verbatim
 135 *>          LDB is INTEGER
 136 *>          The leading dimension of the array B. LDB >= max(1,P).
 137 *> \endverbatim
 138 *>
 139 *> \param[in] TOLA
 140 *> \verbatim
 141 *>          TOLA is DOUBLE PRECISION
 142 *> \endverbatim
 143 *>
 144 *> \param[in] TOLB
 145 *> \verbatim
 146 *>          TOLB is DOUBLE PRECISION
 147 *>
 148 *>          TOLA and TOLB are the thresholds to determine the effective
 149 *>          numerical rank of matrix B and a subblock of A. Generally,
 150 *>          they are set to
 151 *>             TOLA = MAX(M,N)*norm(A)*MACHEPS,
 152 *>             TOLB = MAX(P,N)*norm(B)*MACHEPS.
 153 *>          The size of TOLA and TOLB may affect the size of backward
 154 *>          errors of the decomposition.
 155 *> \endverbatim
 156 *>
 157 *> \param[out] K
 158 *> \verbatim
 159 *>          K is INTEGER
 160 *> \endverbatim
 161 *>
 162 *> \param[out] L
 163 *> \verbatim
 164 *>          L is INTEGER
 165 *>
 166 *>          On exit, K and L specify the dimension of the subblocks
 167 *>          described in Purpose section.
 168 *>          K + L = effective numerical rank of (A**T,B**T)**T.
 169 *> \endverbatim
 170 *>
 171 *> \param[out] U
 172 *> \verbatim
 173 *>          U is DOUBLE PRECISION array, dimension (LDU,M)
 174 *>          If JOBU = 'U', U contains the orthogonal matrix U.
 175 *>          If JOBU = 'N', U is not referenced.
 176 *> \endverbatim
 177 *>
 178 *> \param[in] LDU
 179 *> \verbatim
 180 *>          LDU is INTEGER
 181 *>          The leading dimension of the array U. LDU >= max(1,M) if
 182 *>          JOBU = 'U'; LDU >= 1 otherwise.
 183 *> \endverbatim
 184 *>
 185 *> \param[out] V
 186 *> \verbatim
 187 *>          V is DOUBLE PRECISION array, dimension (LDV,P)
 188 *>          If JOBV = 'V', V contains the orthogonal matrix V.
 189 *>          If JOBV = 'N', V is not referenced.
 190 *> \endverbatim
 191 *>
 192 *> \param[in] LDV
 193 *> \verbatim
 194 *>          LDV is INTEGER
 195 *>          The leading dimension of the array V. LDV >= max(1,P) if
 196 *>          JOBV = 'V'; LDV >= 1 otherwise.
 197 *> \endverbatim
 198 *>
 199 *> \param[out] Q
 200 *> \verbatim
 201 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
 202 *>          If JOBQ = 'Q', Q contains the orthogonal matrix Q.
 203 *>          If JOBQ = 'N', Q is not referenced.
 204 *> \endverbatim
 205 *>
 206 *> \param[in] LDQ
 207 *> \verbatim
 208 *>          LDQ is INTEGER
 209 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
 210 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
 211 *> \endverbatim
 212 *>
 213 *> \param[out] IWORK
 214 *> \verbatim
 215 *>          IWORK is INTEGER array, dimension (N)
 216 *> \endverbatim
 217 *>
 218 *> \param[out] TAU
 219 *> \verbatim
 220 *>          TAU is DOUBLE PRECISION array, dimension (N)
 221 *> \endverbatim
 222 *>
 223 *> \param[out] WORK
 224 *> \verbatim
 225 *>          WORK is DOUBLE PRECISION array, dimension (max(3*N,M,P))
 226 *> \endverbatim
 227 *>
 228 *> \param[out] INFO
 229 *> \verbatim
 230 *>          INFO is INTEGER
 231 *>          = 0:  successful exit
 232 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 233 *> \endverbatim
 234 *
 235 *  Authors:
 236 *  ========
 237 *
 238 *> \author Univ. of Tennessee
 239 *> \author Univ. of California Berkeley
 240 *> \author Univ. of Colorado Denver
 241 *> \author NAG Ltd.
 242 *
 243 *> \date November 2011
 244 *
 245 *> \ingroup doubleOTHERcomputational
 246 *
 247 *> \par Further Details:
 248 *  =====================
 249 *>
 250 *>  The subroutine uses LAPACK subroutine DGEQPF for the QR factorization
 251 *>  with column pivoting to detect the effective numerical rank of the
 252 *>  a matrix. It may be replaced by a better rank determination strategy.
 253 *>
 254 *  =====================================================================
 255       SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
 256      $                   TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
 257      $                   IWORK, TAU, WORK, INFO )
 258 *
 259 *  -- LAPACK computational routine (version 3.4.0) --
 260 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 261 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 262 *     November 2011
 263 *
 264 *     .. Scalar Arguments ..
 265       CHARACTER          JOBQ, JOBU, JOBV
 266       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
 267       DOUBLE PRECISION   TOLA, TOLB
 268 *     ..
 269 *     .. Array Arguments ..
 270       INTEGER            IWORK( * )
 271       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 272      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
 273 *     ..
 274 *
 275 *  =====================================================================
 276 *
 277 *     .. Parameters ..
 278       DOUBLE PRECISION   ZERO, ONE
 279       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 280 *     ..
 281 *     .. Local Scalars ..
 282       LOGICAL            FORWRD, WANTQ, WANTU, WANTV
 283       INTEGER            I, J
 284 *     ..
 285 *     .. External Functions ..
 286       LOGICAL            LSAME
 287       EXTERNAL           LSAME
 288 *     ..
 289 *     .. External Subroutines ..
 290       EXTERNAL           DGEQPF, DGEQR2, DGERQ2, DLACPY, DLAPMT, DLASET,
 291      $                   DORG2R, DORM2R, DORMR2, XERBLA
 292 *     ..
 293 *     .. Intrinsic Functions ..
 294       INTRINSIC          ABS, MAX, MIN
 295 *     ..
 296 *     .. Executable Statements ..
 297 *
 298 *     Test the input parameters
 299 *
 300       WANTU = LSAME( JOBU, 'U' )
 301       WANTV = LSAME( JOBV, 'V' )
 302       WANTQ = LSAME( JOBQ, 'Q' )
 303       FORWRD = .TRUE.
 304 *
 305       INFO = 0
 306       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
 307          INFO = -1
 308       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 309          INFO = -2
 310       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
 311          INFO = -3
 312       ELSE IF( M.LT.0 ) THEN
 313          INFO = -4
 314       ELSE IF( P.LT.0 ) THEN
 315          INFO = -5
 316       ELSE IF( N.LT.0 ) THEN
 317          INFO = -6
 318       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 319          INFO = -8
 320       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
 321          INFO = -10
 322       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
 323          INFO = -16
 324       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
 325          INFO = -18
 326       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
 327          INFO = -20
 328       END IF
 329       IF( INFO.NE.0 ) THEN
 330          CALL XERBLA( 'DGGSVP', -INFO )
 331          RETURN
 332       END IF
 333 *
 334 *     QR with column pivoting of B: B*P = V*( S11 S12 )
 335 *                                           (  0   0  )
 336 *
 337       DO 10 I = 1, N
 338          IWORK( I ) = 0
 339    10 CONTINUE
 340       CALL DGEQPF( P, N, B, LDB, IWORK, TAU, WORK, INFO )
 341 *
 342 *     Update A := A*P
 343 *
 344       CALL DLAPMT( FORWRD, M, N, A, LDA, IWORK )
 345 *
 346 *     Determine the effective rank of matrix B.
 347 *
 348       L = 0
 349       DO 20 I = 1, MIN( P, N )
 350          IF( ABS( B( I, I ) ).GT.TOLB )
 351      $      L = L + 1
 352    20 CONTINUE
 353 *
 354       IF( WANTV ) THEN
 355 *
 356 *        Copy the details of V, and form V.
 357 *
 358          CALL DLASET( 'Full', P, P, ZERO, ZERO, V, LDV )
 359          IF( P.GT.1 )
 360      $      CALL DLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
 361      $                   LDV )
 362          CALL DORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
 363       END IF
 364 *
 365 *     Clean up B
 366 *
 367       DO 40 J = 1, L - 1
 368          DO 30 I = J + 1, L
 369             B( I, J ) = ZERO
 370    30    CONTINUE
 371    40 CONTINUE
 372       IF( P.GT.L )
 373      $   CALL DLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB )
 374 *
 375       IF( WANTQ ) THEN
 376 *
 377 *        Set Q = I and Update Q := Q*P
 378 *
 379          CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
 380          CALL DLAPMT( FORWRD, N, N, Q, LDQ, IWORK )
 381       END IF
 382 *
 383       IF( P.GE.L .AND. N.NE.L ) THEN
 384 *
 385 *        RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
 386 *
 387          CALL DGERQ2( L, N, B, LDB, TAU, WORK, INFO )
 388 *
 389 *        Update A := A*Z**T
 390 *
 391          CALL DORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A,
 392      $                LDA, WORK, INFO )
 393 *
 394          IF( WANTQ ) THEN
 395 *
 396 *           Update Q := Q*Z**T
 397 *
 398             CALL DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,
 399      $                   LDQ, WORK, INFO )
 400          END IF
 401 *
 402 *        Clean up B
 403 *
 404          CALL DLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB )
 405          DO 60 J = N - L + 1, N
 406             DO 50 I = J - N + L + 1, L
 407                B( I, J ) = ZERO
 408    50       CONTINUE
 409    60    CONTINUE
 410 *
 411       END IF
 412 *
 413 *     Let              N-L     L
 414 *                A = ( A11    A12 ) M,
 415 *
 416 *     then the following does the complete QR decomposition of A11:
 417 *
 418 *              A11 = U*(  0  T12 )*P1**T
 419 *                      (  0   0  )
 420 *
 421       DO 70 I = 1, N - L
 422          IWORK( I ) = 0
 423    70 CONTINUE
 424       CALL DGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, INFO )
 425 *
 426 *     Determine the effective rank of A11
 427 *
 428       K = 0
 429       DO 80 I = 1, MIN( M, N-L )
 430          IF( ABS( A( I, I ) ).GT.TOLA )
 431      $      K = K + 1
 432    80 CONTINUE
 433 *
 434 *     Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
 435 *
 436       CALL DORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA,
 437      $             TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
 438 *
 439       IF( WANTU ) THEN
 440 *
 441 *        Copy the details of U, and form U
 442 *
 443          CALL DLASET( 'Full', M, M, ZERO, ZERO, U, LDU )
 444          IF( M.GT.1 )
 445      $      CALL DLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
 446      $                   LDU )
 447          CALL DORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO )
 448       END IF
 449 *
 450       IF( WANTQ ) THEN
 451 *
 452 *        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1
 453 *
 454          CALL DLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK )
 455       END IF
 456 *
 457 *     Clean up A: set the strictly lower triangular part of
 458 *     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
 459 *
 460       DO 100 J = 1, K - 1
 461          DO 90 I = J + 1, K
 462             A( I, J ) = ZERO
 463    90    CONTINUE
 464   100 CONTINUE
 465       IF( M.GT.K )
 466      $   CALL DLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA )
 467 *
 468       IF( N-L.GT.K ) THEN
 469 *
 470 *        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
 471 *
 472          CALL DGERQ2( K, N-L, A, LDA, TAU, WORK, INFO )
 473 *
 474          IF( WANTQ ) THEN
 475 *
 476 *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
 477 *
 478             CALL DORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU,
 479      $                   Q, LDQ, WORK, INFO )
 480          END IF
 481 *
 482 *        Clean up A
 483 *
 484          CALL DLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA )
 485          DO 120 J = N - L - K + 1, N - L
 486             DO 110 I = J - N + L + K + 1, K
 487                A( I, J ) = ZERO
 488   110       CONTINUE
 489   120    CONTINUE
 490 *
 491       END IF
 492 *
 493       IF( M.GT.K ) THEN
 494 *
 495 *        QR factorization of A( K+1:M,N-L+1:N )
 496 *
 497          CALL DGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO )
 498 *
 499          IF( WANTU ) THEN
 500 *
 501 *           Update U(:,K+1:M) := U(:,K+1:M)*U1
 502 *
 503             CALL DORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ),
 504      $                   A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU,
 505      $                   WORK, INFO )
 506          END IF
 507 *
 508 *        Clean up
 509 *
 510          DO 140 J = N - L + 1, N
 511             DO 130 I = J - N + K + L + 1, M
 512                A( I, J ) = ZERO
 513   130       CONTINUE
 514   140    CONTINUE
 515 *
 516       END IF
 517 *
 518       RETURN
 519 *
 520 *     End of DGGSVP
 521 *
 522       END