Annotation of rpl/lapack/lapack/dggsvp3.f, revision 1.6

1.1       bertrand    1: *> \brief \b DGGSVP3
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.3       bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.1       bertrand    7: *
                      8: *> \htmlonly
1.3       bertrand    9: *> Download DGGSVP3 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvp3.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvp3.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvp3.f">
1.1       bertrand   15: *> [TXT]</a>
1.3       bertrand   16: *> \endhtmlonly
1.1       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DGGSVP3( 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 )
1.3       bertrand   24: *
1.1       bertrand   25: *       .. Scalar Arguments ..
                     26: *       CHARACTER          JOBQ, JOBU, JOBV
                     27: *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK
                     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: *       ..
1.3       bertrand   35: *
1.1       bertrand   36: *
                     37: *> \par Purpose:
                     38: *  =============
                     39: *>
                     40: *> \verbatim
                     41: *>
                     42: *> DGGSVP3 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
1.3       bertrand   60: *> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T.
1.1       bertrand   61: *>
                     62: *> This decomposition is the preprocessing step for computing the
                     63: *> Generalized Singular Value Decomposition (GSVD), see subroutine
                     64: *> DGGSVD3.
                     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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
                    140: *> \endverbatim
                    141: *>
                    142: *> \param[in] TOLB
                    143: *> \verbatim
                    144: *>          TOLB is DOUBLE PRECISION
                    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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
                    219: *> \endverbatim
                    220: *>
                    221: *> \param[out] WORK
                    222: *> \verbatim
                    223: *>          WORK is DOUBLE PRECISION 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: *
1.3       bertrand  248: *> \author Univ. of Tennessee
                    249: *> \author Univ. of California Berkeley
                    250: *> \author Univ. of Colorado Denver
                    251: *> \author NAG Ltd.
1.1       bertrand  252: *
                    253: *> \ingroup doubleOTHERcomputational
                    254: *
                    255: *> \par Further Details:
                    256: *  =====================
                    257: *>
                    258: *> \verbatim
                    259: *>
                    260: *>  The subroutine uses LAPACK subroutine DGEQP3 for the QR factorization
                    261: *>  with column pivoting to detect the effective numerical rank of the
                    262: *>  a matrix. It may be replaced by a better rank determination strategy.
                    263: *>
                    264: *>  DGGSVP3 replaces the deprecated subroutine DGGSVP.
                    265: *>
                    266: *> \endverbatim
                    267: *>
                    268: *  =====================================================================
                    269:       SUBROUTINE DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
                    270:      $                    TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
                    271:      $                    IWORK, TAU, WORK, LWORK, INFO )
                    272: *
1.6     ! bertrand  273: *  -- LAPACK computational routine --
1.1       bertrand  274: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    275: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    276: *
                    277:       IMPLICIT NONE
                    278: *
                    279: *     .. Scalar Arguments ..
                    280:       CHARACTER          JOBQ, JOBU, JOBV
                    281:       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
                    282:      $                   LWORK
                    283:       DOUBLE PRECISION   TOLA, TOLB
                    284: *     ..
                    285: *     .. Array Arguments ..
                    286:       INTEGER            IWORK( * )
                    287:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
                    288:      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
                    289: *     ..
                    290: *
                    291: *  =====================================================================
                    292: *
                    293: *     .. Parameters ..
                    294:       DOUBLE PRECISION   ZERO, ONE
                    295:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
                    296: *     ..
                    297: *     .. Local Scalars ..
                    298:       LOGICAL            FORWRD, WANTQ, WANTU, WANTV, LQUERY
                    299:       INTEGER            I, J, LWKOPT
                    300: *     ..
                    301: *     .. External Functions ..
                    302:       LOGICAL            LSAME
                    303:       EXTERNAL           LSAME
                    304: *     ..
                    305: *     .. External Subroutines ..
                    306:       EXTERNAL           DGEQP3, DGEQR2, DGERQ2, DLACPY, DLAPMT,
                    307:      $                   DLASET, DORG2R, DORM2R, DORMR2, XERBLA
                    308: *     ..
                    309: *     .. Intrinsic Functions ..
                    310:       INTRINSIC          ABS, MAX, MIN
                    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:       LQUERY = ( LWORK.EQ.-1 )
                    321:       LWKOPT = 1
                    322: *
                    323: *     Test the input arguments
                    324: *
                    325:       INFO = 0
                    326:       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
                    327:          INFO = -1
                    328:       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
                    329:          INFO = -2
                    330:       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
                    331:          INFO = -3
                    332:       ELSE IF( M.LT.0 ) THEN
                    333:          INFO = -4
                    334:       ELSE IF( P.LT.0 ) THEN
                    335:          INFO = -5
                    336:       ELSE IF( N.LT.0 ) THEN
                    337:          INFO = -6
                    338:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
                    339:          INFO = -8
                    340:       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
                    341:          INFO = -10
                    342:       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
                    343:          INFO = -16
                    344:       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
                    345:          INFO = -18
                    346:       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
                    347:          INFO = -20
                    348:       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
                    349:          INFO = -24
                    350:       END IF
                    351: *
                    352: *     Compute workspace
                    353: *
                    354:       IF( INFO.EQ.0 ) THEN
                    355:          CALL DGEQP3( P, N, B, LDB, IWORK, TAU, WORK, -1, INFO )
                    356:          LWKOPT = INT( WORK ( 1 ) )
                    357:          IF( WANTV ) THEN
                    358:             LWKOPT = MAX( LWKOPT, P )
                    359:          END IF
                    360:          LWKOPT = MAX( LWKOPT, MIN( N, P ) )
                    361:          LWKOPT = MAX( LWKOPT, M )
                    362:          IF( WANTQ ) THEN
                    363:             LWKOPT = MAX( LWKOPT, N )
                    364:          END IF
                    365:          CALL DGEQP3( M, N, A, LDA, IWORK, TAU, WORK, -1, INFO )
                    366:          LWKOPT = MAX( LWKOPT, INT( WORK ( 1 ) ) )
                    367:          LWKOPT = MAX( 1, LWKOPT )
                    368:          WORK( 1 ) = DBLE( LWKOPT )
                    369:       END IF
                    370: *
                    371:       IF( INFO.NE.0 ) THEN
                    372:          CALL XERBLA( 'DGGSVP3', -INFO )
                    373:          RETURN
                    374:       END IF
                    375:       IF( LQUERY ) THEN
                    376:          RETURN
                    377:       ENDIF
                    378: *
                    379: *     QR with column pivoting of B: B*P = V*( S11 S12 )
                    380: *                                           (  0   0  )
                    381: *
                    382:       DO 10 I = 1, N
                    383:          IWORK( I ) = 0
                    384:    10 CONTINUE
                    385:       CALL DGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, INFO )
                    386: *
                    387: *     Update A := A*P
                    388: *
                    389:       CALL DLAPMT( FORWRD, M, N, A, LDA, IWORK )
                    390: *
                    391: *     Determine the effective rank of matrix B.
                    392: *
                    393:       L = 0
                    394:       DO 20 I = 1, MIN( P, N )
                    395:          IF( ABS( B( I, I ) ).GT.TOLB )
                    396:      $      L = L + 1
                    397:    20 CONTINUE
                    398: *
                    399:       IF( WANTV ) THEN
                    400: *
                    401: *        Copy the details of V, and form V.
                    402: *
                    403:          CALL DLASET( 'Full', P, P, ZERO, ZERO, V, LDV )
                    404:          IF( P.GT.1 )
                    405:      $      CALL DLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
                    406:      $                   LDV )
                    407:          CALL DORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
                    408:       END IF
                    409: *
                    410: *     Clean up B
                    411: *
                    412:       DO 40 J = 1, L - 1
                    413:          DO 30 I = J + 1, L
                    414:             B( I, J ) = ZERO
                    415:    30    CONTINUE
                    416:    40 CONTINUE
                    417:       IF( P.GT.L )
                    418:      $   CALL DLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB )
                    419: *
                    420:       IF( WANTQ ) THEN
                    421: *
                    422: *        Set Q = I and Update Q := Q*P
                    423: *
                    424:          CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
                    425:          CALL DLAPMT( FORWRD, N, N, Q, LDQ, IWORK )
                    426:       END IF
                    427: *
                    428:       IF( P.GE.L .AND. N.NE.L ) THEN
                    429: *
                    430: *        RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
                    431: *
                    432:          CALL DGERQ2( L, N, B, LDB, TAU, WORK, INFO )
                    433: *
                    434: *        Update A := A*Z**T
                    435: *
                    436:          CALL DORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A,
                    437:      $                LDA, WORK, INFO )
                    438: *
                    439:          IF( WANTQ ) THEN
                    440: *
                    441: *           Update Q := Q*Z**T
                    442: *
                    443:             CALL DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,
                    444:      $                   LDQ, WORK, INFO )
                    445:          END IF
                    446: *
                    447: *        Clean up B
                    448: *
                    449:          CALL DLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB )
                    450:          DO 60 J = N - L + 1, N
                    451:             DO 50 I = J - N + L + 1, L
                    452:                B( I, J ) = ZERO
                    453:    50       CONTINUE
                    454:    60    CONTINUE
                    455: *
                    456:       END IF
                    457: *
                    458: *     Let              N-L     L
                    459: *                A = ( A11    A12 ) M,
                    460: *
                    461: *     then the following does the complete QR decomposition of A11:
                    462: *
                    463: *              A11 = U*(  0  T12 )*P1**T
                    464: *                      (  0   0  )
                    465: *
                    466:       DO 70 I = 1, N - L
                    467:          IWORK( I ) = 0
                    468:    70 CONTINUE
                    469:       CALL DGEQP3( M, N-L, A, LDA, IWORK, TAU, WORK, LWORK, INFO )
                    470: *
                    471: *     Determine the effective rank of A11
                    472: *
                    473:       K = 0
                    474:       DO 80 I = 1, MIN( M, N-L )
                    475:          IF( ABS( A( I, I ) ).GT.TOLA )
                    476:      $      K = K + 1
                    477:    80 CONTINUE
                    478: *
                    479: *     Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
                    480: *
                    481:       CALL DORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA,
                    482:      $             TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
                    483: *
                    484:       IF( WANTU ) THEN
                    485: *
                    486: *        Copy the details of U, and form U
                    487: *
                    488:          CALL DLASET( 'Full', M, M, ZERO, ZERO, U, LDU )
                    489:          IF( M.GT.1 )
                    490:      $      CALL DLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
                    491:      $                   LDU )
                    492:          CALL DORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO )
                    493:       END IF
                    494: *
                    495:       IF( WANTQ ) THEN
                    496: *
                    497: *        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1
                    498: *
                    499:          CALL DLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK )
                    500:       END IF
                    501: *
                    502: *     Clean up A: set the strictly lower triangular part of
                    503: *     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
                    504: *
                    505:       DO 100 J = 1, K - 1
                    506:          DO 90 I = J + 1, K
                    507:             A( I, J ) = ZERO
                    508:    90    CONTINUE
                    509:   100 CONTINUE
                    510:       IF( M.GT.K )
                    511:      $   CALL DLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA )
                    512: *
                    513:       IF( N-L.GT.K ) THEN
                    514: *
                    515: *        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
                    516: *
                    517:          CALL DGERQ2( K, N-L, A, LDA, TAU, WORK, INFO )
                    518: *
                    519:          IF( WANTQ ) THEN
                    520: *
                    521: *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
                    522: *
                    523:             CALL DORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU,
                    524:      $                   Q, LDQ, WORK, INFO )
                    525:          END IF
                    526: *
                    527: *        Clean up A
                    528: *
                    529:          CALL DLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA )
                    530:          DO 120 J = N - L - K + 1, N - L
                    531:             DO 110 I = J - N + L + K + 1, K
                    532:                A( I, J ) = ZERO
                    533:   110       CONTINUE
                    534:   120    CONTINUE
                    535: *
                    536:       END IF
                    537: *
                    538:       IF( M.GT.K ) THEN
                    539: *
                    540: *        QR factorization of A( K+1:M,N-L+1:N )
                    541: *
                    542:          CALL DGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO )
                    543: *
                    544:          IF( WANTU ) THEN
                    545: *
                    546: *           Update U(:,K+1:M) := U(:,K+1:M)*U1
                    547: *
                    548:             CALL DORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ),
                    549:      $                   A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU,
                    550:      $                   WORK, INFO )
                    551:          END IF
                    552: *
                    553: *        Clean up
                    554: *
                    555:          DO 140 J = N - L + 1, N
                    556:             DO 130 I = J - N + K + L + 1, M
                    557:                A( I, J ) = ZERO
                    558:   130       CONTINUE
                    559:   140    CONTINUE
                    560: *
                    561:       END IF
                    562: *
                    563:       WORK( 1 ) = DBLE( LWKOPT )
                    564:       RETURN
                    565: *
                    566: *     End of DGGSVP3
                    567: *
                    568:       END

CVSweb interface <joel.bertrand@systella.fr>