SRC/sggsvp3.f

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