[BACK]Return to dggsvp.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dggsvp.f between versions 1.6 and 1.19

version 1.6, 2010/08/13 21:03:46 version 1.19, 2023/08/07 08:38:51
Line 1 Line 1
   *> \brief \b DGGSVP
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DGGSVP + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvp.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvp.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvp.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
   *                          TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
   *                          IWORK, TAU, WORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          JOBQ, JOBU, JOBV
   *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
   *       DOUBLE PRECISION   TOLA, TOLB
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IWORK( * )
   *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
   *      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> This routine is deprecated and has been replaced by routine DGGSVP3.
   *>
   *> DGGSVP computes orthogonal matrices U, V and Q such that
   *>
   *>                    N-K-L  K    L
   *>  U**T*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
   *>                 L ( 0     0   A23 )
   *>             M-K-L ( 0     0    0  )
   *>
   *>                  N-K-L  K    L
   *>         =     K ( 0    A12  A13 )  if M-K-L < 0;
   *>             M-K ( 0     0   A23 )
   *>
   *>                  N-K-L  K    L
   *>  V**T*B*Q =   L ( 0     0   B13 )
   *>             P-L ( 0     0    0  )
   *>
   *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
   *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
   *> otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
   *> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T.
   *>
   *> This decomposition is the preprocessing step for computing the
   *> Generalized Singular Value Decomposition (GSVD), see subroutine
   *> DGGSVD.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] JOBU
   *> \verbatim
   *>          JOBU is CHARACTER*1
   *>          = 'U':  Orthogonal matrix U is computed;
   *>          = 'N':  U is not computed.
   *> \endverbatim
   *>
   *> \param[in] JOBV
   *> \verbatim
   *>          JOBV is CHARACTER*1
   *>          = 'V':  Orthogonal matrix V is computed;
   *>          = 'N':  V is not computed.
   *> \endverbatim
   *>
   *> \param[in] JOBQ
   *> \verbatim
   *>          JOBQ is CHARACTER*1
   *>          = 'Q':  Orthogonal matrix Q is computed;
   *>          = 'N':  Q is not computed.
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The number of rows of the matrix A.  M >= 0.
   *> \endverbatim
   *>
   *> \param[in] P
   *> \verbatim
   *>          P is INTEGER
   *>          The number of rows of the matrix B.  P >= 0.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The number of columns of the matrices A and B.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] A
   *> \verbatim
   *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   *>          On entry, the M-by-N matrix A.
   *>          On exit, A contains the triangular (or trapezoidal) matrix
   *>          described in the Purpose section.
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A. LDA >= max(1,M).
   *> \endverbatim
   *>
   *> \param[in,out] B
   *> \verbatim
   *>          B is DOUBLE PRECISION array, dimension (LDB,N)
   *>          On entry, the P-by-N matrix B.
   *>          On exit, B contains the triangular matrix described in
   *>          the Purpose section.
   *> \endverbatim
   *>
   *> \param[in] LDB
   *> \verbatim
   *>          LDB is INTEGER
   *>          The leading dimension of the array B. LDB >= max(1,P).
   *> \endverbatim
   *>
   *> \param[in] TOLA
   *> \verbatim
   *>          TOLA is DOUBLE PRECISION
   *> \endverbatim
   *>
   *> \param[in] TOLB
   *> \verbatim
   *>          TOLB is DOUBLE PRECISION
   *>
   *>          TOLA and TOLB are the thresholds to determine the effective
   *>          numerical rank of matrix B and a subblock of A. Generally,
   *>          they are set to
   *>             TOLA = MAX(M,N)*norm(A)*MACHEPS,
   *>             TOLB = MAX(P,N)*norm(B)*MACHEPS.
   *>          The size of TOLA and TOLB may affect the size of backward
   *>          errors of the decomposition.
   *> \endverbatim
   *>
   *> \param[out] K
   *> \verbatim
   *>          K is INTEGER
   *> \endverbatim
   *>
   *> \param[out] L
   *> \verbatim
   *>          L is INTEGER
   *>
   *>          On exit, K and L specify the dimension of the subblocks
   *>          described in Purpose section.
   *>          K + L = effective numerical rank of (A**T,B**T)**T.
   *> \endverbatim
   *>
   *> \param[out] U
   *> \verbatim
   *>          U is DOUBLE PRECISION array, dimension (LDU,M)
   *>          If JOBU = 'U', U contains the orthogonal matrix U.
   *>          If JOBU = 'N', U is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDU
   *> \verbatim
   *>          LDU is INTEGER
   *>          The leading dimension of the array U. LDU >= max(1,M) if
   *>          JOBU = 'U'; LDU >= 1 otherwise.
   *> \endverbatim
   *>
   *> \param[out] V
   *> \verbatim
   *>          V is DOUBLE PRECISION array, dimension (LDV,P)
   *>          If JOBV = 'V', V contains the orthogonal matrix V.
   *>          If JOBV = 'N', V is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDV
   *> \verbatim
   *>          LDV is INTEGER
   *>          The leading dimension of the array V. LDV >= max(1,P) if
   *>          JOBV = 'V'; LDV >= 1 otherwise.
   *> \endverbatim
   *>
   *> \param[out] Q
   *> \verbatim
   *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
   *>          If JOBQ = 'Q', Q contains the orthogonal matrix Q.
   *>          If JOBQ = 'N', Q is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDQ
   *> \verbatim
   *>          LDQ is INTEGER
   *>          The leading dimension of the array Q. LDQ >= max(1,N) if
   *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (N)
   *> \endverbatim
   *>
   *> \param[out] TAU
   *> \verbatim
   *>          TAU is DOUBLE PRECISION array, dimension (N)
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (max(3*N,M,P))
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \ingroup doubleOTHERcomputational
   *
   *> \par Further Details:
   *  =====================
   *>
   *>  The subroutine uses LAPACK subroutine DGEQPF for the QR factorization
   *>  with column pivoting to detect the effective numerical rank of the
   *>  a matrix. It may be replaced by a better rank determination strategy.
   *>
   *  =====================================================================
       SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,        SUBROUTINE DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
      $                   TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,       $                   TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
      $                   IWORK, TAU, WORK, INFO )       $                   IWORK, TAU, WORK, INFO )
 *  *
 *  -- LAPACK routine (version 3.2) --  *  -- LAPACK computational routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2006  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBQ, JOBU, JOBV        CHARACTER          JOBQ, JOBU, JOBV
Line 18 Line 269
      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )       $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DGGSVP computes orthogonal matrices U, V and Q such that  
 *  
 *                   N-K-L  K    L  
 *   U'*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;  
 *                L ( 0     0   A23 )  
 *            M-K-L ( 0     0    0  )  
 *  
 *                   N-K-L  K    L  
 *          =     K ( 0    A12  A13 )  if M-K-L < 0;  
 *              M-K ( 0     0   A23 )  
 *  
 *                 N-K-L  K    L  
 *   V'*B*Q =   L ( 0     0   B13 )  
 *            P-L ( 0     0    0  )  
 *  
 *  where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular  
 *  upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,  
 *  otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective  
 *  numerical rank of the (M+P)-by-N matrix (A',B')'.  Z' denotes the  
 *  transpose of Z.  
 *  
 *  This decomposition is the preprocessing step for computing the  
 *  Generalized Singular Value Decomposition (GSVD), see subroutine  
 *  DGGSVD.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  JOBU    (input) CHARACTER*1  
 *          = 'U':  Orthogonal matrix U is computed;  
 *          = 'N':  U is not computed.  
 *  
 *  JOBV    (input) CHARACTER*1  
 *          = 'V':  Orthogonal matrix V is computed;  
 *          = 'N':  V is not computed.  
 *  
 *  JOBQ    (input) CHARACTER*1  
 *          = 'Q':  Orthogonal matrix Q is computed;  
 *          = 'N':  Q is not computed.  
 *  
 *  M       (input) INTEGER  
 *          The number of rows of the matrix A.  M >= 0.  
 *  
 *  P       (input) INTEGER  
 *          The number of rows of the matrix B.  P >= 0.  
 *  
 *  N       (input) INTEGER  
 *          The number of columns of the matrices A and B.  N >= 0.  
 *  
 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)  
 *          On entry, the M-by-N matrix A.  
 *          On exit, A contains the triangular (or trapezoidal) matrix  
 *          described in the Purpose section.  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A. LDA >= max(1,M).  
 *  
 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)  
 *          On entry, the P-by-N matrix B.  
 *          On exit, B contains the triangular matrix described in  
 *          the Purpose section.  
 *  
 *  LDB     (input) INTEGER  
 *          The leading dimension of the array B. LDB >= max(1,P).  
 *  
 *  TOLA    (input) DOUBLE PRECISION  
 *  TOLB    (input) DOUBLE PRECISION  
 *          TOLA and TOLB are the thresholds to determine the effective  
 *          numerical rank of matrix B and a subblock of A. Generally,  
 *          they are set to  
 *             TOLA = MAX(M,N)*norm(A)*MAZHEPS,  
 *             TOLB = MAX(P,N)*norm(B)*MAZHEPS.  
 *          The size of TOLA and TOLB may affect the size of backward  
 *          errors of the decomposition.  
 *  
 *  K       (output) INTEGER  
 *  L       (output) INTEGER  
 *          On exit, K and L specify the dimension of the subblocks  
 *          described in Purpose.  
 *          K + L = effective numerical rank of (A',B')'.  
 *  
 *  U       (output) DOUBLE PRECISION array, dimension (LDU,M)  
 *          If JOBU = 'U', U contains the orthogonal matrix U.  
 *          If JOBU = 'N', U is not referenced.  
 *  
 *  LDU     (input) INTEGER  
 *          The leading dimension of the array U. LDU >= max(1,M) if  
 *          JOBU = 'U'; LDU >= 1 otherwise.  
 *  
 *  V       (output) DOUBLE PRECISION array, dimension (LDV,P)  
 *          If JOBV = 'V', V contains the orthogonal matrix V.  
 *          If JOBV = 'N', V is not referenced.  
 *  
 *  LDV     (input) INTEGER  
 *          The leading dimension of the array V. LDV >= max(1,P) if  
 *          JOBV = 'V'; LDV >= 1 otherwise.  
 *  
 *  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)  
 *          If JOBQ = 'Q', Q contains the orthogonal matrix Q.  
 *          If JOBQ = 'N', Q is not referenced.  
 *  
 *  LDQ     (input) INTEGER  
 *          The leading dimension of the array Q. LDQ >= max(1,N) if  
 *          JOBQ = 'Q'; LDQ >= 1 otherwise.  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (N)  
 *  
 *  TAU     (workspace) DOUBLE PRECISION array, dimension (N)  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P))  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value.  
 *  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  The subroutine uses LAPACK subroutine DGEQPF for the QR factorization  
 *  with column pivoting to detect the effective numerical rank of the  
 *  a matrix. It may be replaced by a better rank determination strategy.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 258 Line 383
 *  *
          CALL DGERQ2( L, N, B, LDB, TAU, WORK, INFO )           CALL DGERQ2( L, N, B, LDB, TAU, WORK, INFO )
 *  *
 *        Update A := A*Z'  *        Update A := A*Z**T
 *  *
          CALL DORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A,           CALL DORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A,
      $                LDA, WORK, INFO )       $                LDA, WORK, INFO )
 *  *
          IF( WANTQ ) THEN           IF( WANTQ ) THEN
 *  *
 *           Update Q := Q*Z'  *           Update Q := Q*Z**T
 *  *
             CALL DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,              CALL DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,
      $                   LDQ, WORK, INFO )       $                   LDQ, WORK, INFO )
Line 287 Line 412
 *  *
 *     then the following does the complete QR decomposition of A11:  *     then the following does the complete QR decomposition of A11:
 *  *
 *              A11 = U*(  0  T12 )*P1'  *              A11 = U*(  0  T12 )*P1**T
 *                      (  0   0  )  *                      (  0   0  )
 *  *
       DO 70 I = 1, N - L        DO 70 I = 1, N - L
Line 303 Line 428
      $      K = K + 1       $      K = K + 1
    80 CONTINUE     80 CONTINUE
 *  *
 *     Update A12 := U'*A12, where A12 = A( 1:M, N-L+1:N )  *     Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
 *  *
       CALL DORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA,        CALL DORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA,
      $             TAU, A( 1, N-L+1 ), LDA, WORK, INFO )       $             TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
Line 345 Line 470
 *  *
          IF( WANTQ ) THEN           IF( WANTQ ) THEN
 *  *
 *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1'  *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
 *  *
             CALL DORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU,              CALL DORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU,
      $                   Q, LDQ, WORK, INFO )       $                   Q, LDQ, WORK, INFO )
Removed from v.1.6  
changed lines
  Added in v.1.19

CVSweb interface <joel.bertrand@systella.fr>