Return to zggsvd.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

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