Return to zggsvp.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

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