SRC/dggsvd3.f

   1 *> \brief <b> DGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGGSVD3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvd3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvd3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvd3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
  22 *                           LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
  23 *                           LWORK, IWORK, 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 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IWORK( * )
  31 *       DOUBLE PRECISION   A( LDA, * ), ALPHA( * ), B( LDB, * ),
  32 *      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
  33 *      $                   V( LDV, * ), WORK( * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> DGGSVD3 computes the generalized singular value decomposition (GSVD)
  43 *> of an M-by-N real matrix A and P-by-N real matrix B:
  44 *>
  45 *>       U**T*A*Q = D1*( 0 R ),    V**T*B*Q = D2*( 0 R )
  46 *>
  47 *> where U, V and Q are orthogonal matrices.
  48 *> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T,
  49 *> then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
  50 *> D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
  51 *> following structures, respectively:
  52 *>
  53 *> If M-K-L >= 0,
  54 *>
  55 *>                     K  L
  56 *>        D1 =     K ( I  0 )
  57 *>                 L ( 0  C )
  58 *>             M-K-L ( 0  0 )
  59 *>
  60 *>                   K  L
  61 *>        D2 =   L ( 0  S )
  62 *>             P-L ( 0  0 )
  63 *>
  64 *>                 N-K-L  K    L
  65 *>   ( 0 R ) = K (  0   R11  R12 )
  66 *>             L (  0    0   R22 )
  67 *>
  68 *> where
  69 *>
  70 *>   C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
  71 *>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
  72 *>   C**2 + S**2 = I.
  73 *>
  74 *>   R is stored in A(1:K+L,N-K-L+1:N) on exit.
  75 *>
  76 *> If M-K-L < 0,
  77 *>
  78 *>                   K M-K K+L-M
  79 *>        D1 =   K ( I  0    0   )
  80 *>             M-K ( 0  C    0   )
  81 *>
  82 *>                     K M-K K+L-M
  83 *>        D2 =   M-K ( 0  S    0  )
  84 *>             K+L-M ( 0  0    I  )
  85 *>               P-L ( 0  0    0  )
  86 *>
  87 *>                    N-K-L  K   M-K  K+L-M
  88 *>   ( 0 R ) =     K ( 0    R11  R12  R13  )
  89 *>               M-K ( 0     0   R22  R23  )
  90 *>             K+L-M ( 0     0    0   R33  )
  91 *>
  92 *> where
  93 *>
  94 *>   C = diag( ALPHA(K+1), ... , ALPHA(M) ),
  95 *>   S = diag( BETA(K+1),  ... , BETA(M) ),
  96 *>   C**2 + S**2 = I.
  97 *>
  98 *>   (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
  99 *>   ( 0  R22 R23 )
 100 *>   in B(M-K+1:L,N+M-K-L+1:N) on exit.
 101 *>
 102 *> The routine computes C, S, R, and optionally the orthogonal
 103 *> transformation matrices U, V and Q.
 104 *>
 105 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
 106 *> A and B implicitly gives the SVD of A*inv(B):
 107 *>                      A*inv(B) = U*(D1*inv(D2))*V**T.
 108 *> If ( A**T,B**T)**T  has orthonormal columns, then the GSVD of A and B is
 109 *> also equal to the CS decomposition of A and B. Furthermore, the GSVD
 110 *> can be used to derive the solution of the eigenvalue problem:
 111 *>                      A**T*A x = lambda* B**T*B x.
 112 *> In some literature, the GSVD of A and B is presented in the form
 113 *>                  U**T*A*X = ( 0 D1 ),   V**T*B*X = ( 0 D2 )
 114 *> where U and V are orthogonal and X is nonsingular, D1 and D2 are
 115 *> ``diagonal''.  The former GSVD form can be converted to the latter
 116 *> form by taking the nonsingular matrix X as
 117 *>
 118 *>                      X = Q*( I   0    )
 119 *>                            ( 0 inv(R) ).
 120 *> \endverbatim
 121 *
 122 *  Arguments:
 123 *  ==========
 124 *
 125 *> \param[in] JOBU
 126 *> \verbatim
 127 *>          JOBU is CHARACTER*1
 128 *>          = 'U':  Orthogonal matrix U is computed;
 129 *>          = 'N':  U is not computed.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] JOBV
 133 *> \verbatim
 134 *>          JOBV is CHARACTER*1
 135 *>          = 'V':  Orthogonal matrix V is computed;
 136 *>          = 'N':  V is not computed.
 137 *> \endverbatim
 138 *>
 139 *> \param[in] JOBQ
 140 *> \verbatim
 141 *>          JOBQ is CHARACTER*1
 142 *>          = 'Q':  Orthogonal matrix Q is computed;
 143 *>          = 'N':  Q is not computed.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] M
 147 *> \verbatim
 148 *>          M is INTEGER
 149 *>          The number of rows of the matrix A.  M >= 0.
 150 *> \endverbatim
 151 *>
 152 *> \param[in] N
 153 *> \verbatim
 154 *>          N is INTEGER
 155 *>          The number of columns of the matrices A and B.  N >= 0.
 156 *> \endverbatim
 157 *>
 158 *> \param[in] P
 159 *> \verbatim
 160 *>          P is INTEGER
 161 *>          The number of rows of the matrix B.  P >= 0.
 162 *> \endverbatim
 163 *>
 164 *> \param[out] K
 165 *> \verbatim
 166 *>          K is INTEGER
 167 *> \endverbatim
 168 *>
 169 *> \param[out] L
 170 *> \verbatim
 171 *>          L is INTEGER
 172 *>
 173 *>          On exit, K and L specify the dimension of the subblocks
 174 *>          described in Purpose.
 175 *>          K + L = effective numerical rank of (A**T,B**T)**T.
 176 *> \endverbatim
 177 *>
 178 *> \param[in,out] A
 179 *> \verbatim
 180 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 181 *>          On entry, the M-by-N matrix A.
 182 *>          On exit, A contains the triangular matrix R, or part of R.
 183 *>          See Purpose for details.
 184 *> \endverbatim
 185 *>
 186 *> \param[in] LDA
 187 *> \verbatim
 188 *>          LDA is INTEGER
 189 *>          The leading dimension of the array A. LDA >= max(1,M).
 190 *> \endverbatim
 191 *>
 192 *> \param[in,out] B
 193 *> \verbatim
 194 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
 195 *>          On entry, the P-by-N matrix B.
 196 *>          On exit, B contains the triangular matrix R if M-K-L < 0.
 197 *>          See Purpose for details.
 198 *> \endverbatim
 199 *>
 200 *> \param[in] LDB
 201 *> \verbatim
 202 *>          LDB is INTEGER
 203 *>          The leading dimension of the array B. LDB >= max(1,P).
 204 *> \endverbatim
 205 *>
 206 *> \param[out] ALPHA
 207 *> \verbatim
 208 *>          ALPHA is DOUBLE PRECISION array, dimension (N)
 209 *> \endverbatim
 210 *>
 211 *> \param[out] BETA
 212 *> \verbatim
 213 *>          BETA is DOUBLE PRECISION array, dimension (N)
 214 *>
 215 *>          On exit, ALPHA and BETA contain the generalized singular
 216 *>          value pairs of A and B;
 217 *>            ALPHA(1:K) = 1,
 218 *>            BETA(1:K)  = 0,
 219 *>          and if M-K-L >= 0,
 220 *>            ALPHA(K+1:K+L) = C,
 221 *>            BETA(K+1:K+L)  = S,
 222 *>          or if M-K-L < 0,
 223 *>            ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
 224 *>            BETA(K+1:M) =S, BETA(M+1:K+L) =1
 225 *>          and
 226 *>            ALPHA(K+L+1:N) = 0
 227 *>            BETA(K+L+1:N)  = 0
 228 *> \endverbatim
 229 *>
 230 *> \param[out] U
 231 *> \verbatim
 232 *>          U is DOUBLE PRECISION array, dimension (LDU,M)
 233 *>          If JOBU = 'U', U contains the M-by-M orthogonal matrix U.
 234 *>          If JOBU = 'N', U is not referenced.
 235 *> \endverbatim
 236 *>
 237 *> \param[in] LDU
 238 *> \verbatim
 239 *>          LDU is INTEGER
 240 *>          The leading dimension of the array U. LDU >= max(1,M) if
 241 *>          JOBU = 'U'; LDU >= 1 otherwise.
 242 *> \endverbatim
 243 *>
 244 *> \param[out] V
 245 *> \verbatim
 246 *>          V is DOUBLE PRECISION array, dimension (LDV,P)
 247 *>          If JOBV = 'V', V contains the P-by-P orthogonal matrix V.
 248 *>          If JOBV = 'N', V is not referenced.
 249 *> \endverbatim
 250 *>
 251 *> \param[in] LDV
 252 *> \verbatim
 253 *>          LDV is INTEGER
 254 *>          The leading dimension of the array V. LDV >= max(1,P) if
 255 *>          JOBV = 'V'; LDV >= 1 otherwise.
 256 *> \endverbatim
 257 *>
 258 *> \param[out] Q
 259 *> \verbatim
 260 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
 261 *>          If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.
 262 *>          If JOBQ = 'N', Q is not referenced.
 263 *> \endverbatim
 264 *>
 265 *> \param[in] LDQ
 266 *> \verbatim
 267 *>          LDQ is INTEGER
 268 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
 269 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
 270 *> \endverbatim
 271 *>
 272 *> \param[out] WORK
 273 *> \verbatim
 274 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 275 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 276 *> \endverbatim
 277 *>
 278 *> \param[in] LWORK
 279 *> \verbatim
 280 *>          LWORK is INTEGER
 281 *>          The dimension of the array WORK.
 282 *>
 283 *>          If LWORK = -1, then a workspace query is assumed; the routine
 284 *>          only calculates the optimal size of the WORK array, returns
 285 *>          this value as the first entry of the WORK array, and no error
 286 *>          message related to LWORK is issued by XERBLA.
 287 *> \endverbatim
 288 *>
 289 *> \param[out] IWORK
 290 *> \verbatim
 291 *>          IWORK is INTEGER array, dimension (N)
 292 *>          On exit, IWORK stores the sorting information. More
 293 *>          precisely, the following loop will sort ALPHA
 294 *>             for I = K+1, min(M,K+L)
 295 *>                 swap ALPHA(I) and ALPHA(IWORK(I))
 296 *>             endfor
 297 *>          such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
 298 *> \endverbatim
 299 *>
 300 *> \param[out] INFO
 301 *> \verbatim
 302 *>          INFO is INTEGER
 303 *>          = 0:  successful exit.
 304 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 305 *>          > 0:  if INFO = 1, the Jacobi-type procedure failed to
 306 *>                converge.  For further details, see subroutine DTGSJA.
 307 *> \endverbatim
 308 *
 309 *> \par Internal Parameters:
 310 *  =========================
 311 *>
 312 *> \verbatim
 313 *>  TOLA    DOUBLE PRECISION
 314 *>  TOLB    DOUBLE PRECISION
 315 *>          TOLA and TOLB are the thresholds to determine the effective
 316 *>          rank of (A**T,B**T)**T. Generally, they are set to
 317 *>                   TOLA = MAX(M,N)*norm(A)*MACHEPS,
 318 *>                   TOLB = MAX(P,N)*norm(B)*MACHEPS.
 319 *>          The size of TOLA and TOLB may affect the size of backward
 320 *>          errors of the decomposition.
 321 *> \endverbatim
 322 *
 323 *  Authors:
 324 *  ========
 325 *
 326 *> \author Univ. of Tennessee
 327 *> \author Univ. of California Berkeley
 328 *> \author Univ. of Colorado Denver
 329 *> \author NAG Ltd.
 330 *
 331 *> \date August 2015
 332 *
 333 *> \ingroup doubleGEsing
 334 *
 335 *> \par Contributors:
 336 *  ==================
 337 *>
 338 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 339 *>     California at Berkeley, USA
 340 *>
 341 *
 342 *> \par Further Details:
 343 *  =====================
 344 *>
 345 *>  DGGSVD3 replaces the deprecated subroutine DGGSVD.
 346 *>
 347 *  =====================================================================
 348       SUBROUTINE DGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
 349      $                    LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
 350      $                    WORK, LWORK, IWORK, INFO )
 351 *
 352 *  -- LAPACK driver routine (version 3.6.0) --
 353 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 354 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 355 *     August 2015
 356 *
 357 *     .. Scalar Arguments ..
 358       CHARACTER          JOBQ, JOBU, JOBV
 359       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
 360      $                   LWORK
 361 *     ..
 362 *     .. Array Arguments ..
 363       INTEGER            IWORK( * )
 364       DOUBLE PRECISION   A( LDA, * ), ALPHA( * ), B( LDB, * ),
 365      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
 366      $                   V( LDV, * ), WORK( * )
 367 *     ..
 368 *
 369 *  =====================================================================
 370 *
 371 *     .. Local Scalars ..
 372       LOGICAL            WANTQ, WANTU, WANTV, LQUERY
 373       INTEGER            I, IBND, ISUB, J, NCYCLE, LWKOPT
 374       DOUBLE PRECISION   ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
 375 *     ..
 376 *     .. External Functions ..
 377       LOGICAL            LSAME
 378       DOUBLE PRECISION   DLAMCH, DLANGE
 379       EXTERNAL           LSAME, DLAMCH, DLANGE
 380 *     ..
 381 *     .. External Subroutines ..
 382       EXTERNAL           DCOPY, DGGSVP3, DTGSJA, XERBLA
 383 *     ..
 384 *     .. Intrinsic Functions ..
 385       INTRINSIC          MAX, MIN
 386 *     ..
 387 *     .. Executable Statements ..
 388 *
 389 *     Decode and test the input parameters
 390 *
 391       WANTU = LSAME( JOBU, 'U' )
 392       WANTV = LSAME( JOBV, 'V' )
 393       WANTQ = LSAME( JOBQ, 'Q' )
 394       LQUERY = ( LWORK.EQ.-1 )
 395       LWKOPT = 1
 396 *
 397 *     Test the input arguments
 398 *
 399       INFO = 0
 400       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
 401          INFO = -1
 402       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 403          INFO = -2
 404       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
 405          INFO = -3
 406       ELSE IF( M.LT.0 ) THEN
 407          INFO = -4
 408       ELSE IF( N.LT.0 ) THEN
 409          INFO = -5
 410       ELSE IF( P.LT.0 ) THEN
 411          INFO = -6
 412       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 413          INFO = -10
 414       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
 415          INFO = -12
 416       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
 417          INFO = -16
 418       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
 419          INFO = -18
 420       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
 421          INFO = -20
 422       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
 423          INFO = -24
 424       END IF
 425 *
 426 *     Compute workspace
 427 *
 428       IF( INFO.EQ.0 ) THEN
 429          CALL DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
 430      $                 TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK,
 431      $                 WORK, -1, INFO )
 432          LWKOPT = N + INT( WORK( 1 ) )
 433          LWKOPT = MAX( 2*N, LWKOPT )
 434          LWKOPT = MAX( 1, LWKOPT )
 435          WORK( 1 ) = DBLE( LWKOPT )
 436       END IF
 437 *
 438       IF( INFO.NE.0 ) THEN
 439          CALL XERBLA( 'DGGSVD3', -INFO )
 440          RETURN
 441       END IF
 442       IF( LQUERY ) THEN
 443          RETURN
 444       ENDIF
 445 *
 446 *     Compute the Frobenius norm of matrices A and B
 447 *
 448       ANORM = DLANGE( '1', M, N, A, LDA, WORK )
 449       BNORM = DLANGE( '1', P, N, B, LDB, WORK )
 450 *
 451 *     Get machine precision and set up threshold for determining
 452 *     the effective numerical rank of the matrices A and B.
 453 *
 454       ULP = DLAMCH( 'Precision' )
 455       UNFL = DLAMCH( 'Safe Minimum' )
 456       TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP
 457       TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP
 458 *
 459 *     Preprocessing
 460 *
 461       CALL DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
 462      $              TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK,
 463      $              WORK( N+1 ), LWORK-N, INFO )
 464 *
 465 *     Compute the GSVD of two upper "triangular" matrices
 466 *
 467       CALL DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
 468      $             TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
 469      $             WORK, NCYCLE, INFO )
 470 *
 471 *     Sort the singular values and store the pivot indices in IWORK
 472 *     Copy ALPHA to WORK, then sort ALPHA in WORK
 473 *
 474       CALL DCOPY( N, ALPHA, 1, WORK, 1 )
 475       IBND = MIN( L, M-K )
 476       DO 20 I = 1, IBND
 477 *
 478 *        Scan for largest ALPHA(K+I)
 479 *
 480          ISUB = I
 481          SMAX = WORK( K+I )
 482          DO 10 J = I + 1, IBND
 483             TEMP = WORK( K+J )
 484             IF( TEMP.GT.SMAX ) THEN
 485                ISUB = J
 486                SMAX = TEMP
 487             END IF
 488    10    CONTINUE
 489          IF( ISUB.NE.I ) THEN
 490             WORK( K+ISUB ) = WORK( K+I )
 491             WORK( K+I ) = SMAX
 492             IWORK( K+I ) = K + ISUB
 493          ELSE
 494             IWORK( K+I ) = K + I
 495          END IF
 496    20 CONTINUE
 497 *
 498       WORK( 1 ) = DBLE( LWKOPT )
 499       RETURN
 500 *
 501 *     End of DGGSVD3
 502 *
 503       END