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

1.1     ! bertrand    1: *> \brief \b DGGSVP3
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at 
        !             6: *            http://www.netlib.org/lapack/explore-html/ 
        !             7: *
        !             8: *> \htmlonly
        !             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"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            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 )
        !            24: * 
        !            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: *       ..
        !            35: *  
        !            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
        !            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: *> 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: *
        !           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 doubleOTHERcomputational
        !           256: *
        !           257: *> \par Further Details:
        !           258: *  =====================
        !           259: *>
        !           260: *> \verbatim
        !           261: *>
        !           262: *>  The subroutine uses LAPACK subroutine DGEQP3 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: *>  DGGSVP3 replaces the deprecated subroutine DGGSVP.
        !           267: *>
        !           268: *> \endverbatim
        !           269: *>
        !           270: *  =====================================================================
        !           271:       SUBROUTINE DGGSVP3( 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.0) --
        !           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:       DOUBLE PRECISION   TOLA, TOLB
        !           287: *     ..
        !           288: *     .. Array Arguments ..
        !           289:       INTEGER            IWORK( * )
        !           290:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
        !           291:      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
        !           292: *     ..
        !           293: *
        !           294: *  =====================================================================
        !           295: *
        !           296: *     .. Parameters ..
        !           297:       DOUBLE PRECISION   ZERO, ONE
        !           298:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+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           DGEQP3, DGEQR2, DGERQ2, DLACPY, DLAPMT,
        !           310:      $                   DLASET, DORG2R, DORM2R, DORMR2, 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 DGEQP3( 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 DGEQP3( M, N, A, LDA, IWORK, TAU, WORK, -1, INFO )
        !           369:          LWKOPT = MAX( LWKOPT, INT( WORK ( 1 ) ) )
        !           370:          LWKOPT = MAX( 1, LWKOPT )
        !           371:          WORK( 1 ) = DBLE( LWKOPT )
        !           372:       END IF
        !           373: *
        !           374:       IF( INFO.NE.0 ) THEN
        !           375:          CALL XERBLA( 'DGGSVP3', -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 DGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, INFO )
        !           389: *
        !           390: *     Update A := A*P
        !           391: *
        !           392:       CALL DLAPMT( 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 DLASET( 'Full', P, P, ZERO, ZERO, V, LDV )
        !           407:          IF( P.GT.1 )
        !           408:      $      CALL DLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
        !           409:      $                   LDV )
        !           410:          CALL DORG2R( 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 DLASET( '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 DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
        !           428:          CALL DLAPMT( 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 DGERQ2( L, N, B, LDB, TAU, WORK, INFO )
        !           436: *
        !           437: *        Update A := A*Z**T
        !           438: *
        !           439:          CALL DORMR2( '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 DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,
        !           447:      $                   LDQ, WORK, INFO )
        !           448:          END IF
        !           449: *
        !           450: *        Clean up B
        !           451: *
        !           452:          CALL DLASET( '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 DGEQP3( 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 DORM2R( '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 DLASET( 'Full', M, M, ZERO, ZERO, U, LDU )
        !           492:          IF( M.GT.1 )
        !           493:      $      CALL DLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
        !           494:      $                   LDU )
        !           495:          CALL DORG2R( 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 DLAPMT( 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 DLASET( '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 DGERQ2( 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 DORMR2( '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 DLASET( '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 DGEQR2( 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 DORM2R( '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 ) = DBLE( LWKOPT )
        !           567:       RETURN
        !           568: *
        !           569: *     End of DGGSVP3
        !           570: *
        !           571:       END

CVSweb interface <joel.bertrand@systella.fr>