Return to zggqrf.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
    2:      $                   LWORK, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
   11: *     ..
   12: *     .. Array Arguments ..
   13:       COMPLEX*16         A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
   14:      $                   WORK( * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  ZGGQRF computes a generalized QR factorization of an N-by-M matrix A
   21: *  and an N-by-P matrix B:
   22: *
   23: *              A = Q*R,        B = Q*T*Z,
   24: *
   25: *  where Q is an N-by-N unitary matrix, Z is a P-by-P unitary matrix,
   26: *  and R and T assume one of the forms:
   27: *
   28: *  if N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N,
   29: *                  (  0  ) N-M                         N   M-N
   30: *                     M
   31: *
   32: *  where R11 is upper triangular, and
   33: *
   34: *  if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N-P,
   35: *                   P-N  N                           ( T21 ) P
   36: *                                                       P
   37: *
   38: *  where T12 or T21 is upper triangular.
   39: *
   40: *  In particular, if B is square and nonsingular, the GQR factorization
   41: *  of A and B implicitly gives the QR factorization of inv(B)*A:
   42: *
   43: *               inv(B)*A = Z'*(inv(T)*R)
   44: *
   45: *  where inv(B) denotes the inverse of the matrix B, and Z' denotes the
   46: *  conjugate transpose of matrix Z.
   47: *
   48: *  Arguments
   49: *  =========
   50: *
   51: *  N       (input) INTEGER
   52: *          The number of rows of the matrices A and B. N >= 0.
   53: *
   54: *  M       (input) INTEGER
   55: *          The number of columns of the matrix A.  M >= 0.
   56: *
   57: *  P       (input) INTEGER
   58: *          The number of columns of the matrix B.  P >= 0.
   59: *
   60: *  A       (input/output) COMPLEX*16 array, dimension (LDA,M)
   61: *          On entry, the N-by-M matrix A.
   62: *          On exit, the elements on and above the diagonal of the array
   63: *          contain the min(N,M)-by-M upper trapezoidal matrix R (R is
   64: *          upper triangular if N >= M); the elements below the diagonal,
   65: *          with the array TAUA, represent the unitary matrix Q as a
   66: *          product of min(N,M) elementary reflectors (see Further
   67: *          Details).
   68: *
   69: *  LDA     (input) INTEGER
   70: *          The leading dimension of the array A. LDA >= max(1,N).
   71: *
   72: *  TAUA    (output) COMPLEX*16 array, dimension (min(N,M))
   73: *          The scalar factors of the elementary reflectors which
   74: *          represent the unitary matrix Q (see Further Details).
   75: *
   76: *  B       (input/output) COMPLEX*16 array, dimension (LDB,P)
   77: *          On entry, the N-by-P matrix B.
   78: *          On exit, if N <= P, the upper triangle of the subarray
   79: *          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
   80: *          if N > P, the elements on and above the (N-P)-th subdiagonal
   81: *          contain the N-by-P upper trapezoidal matrix T; the remaining
   82: *          elements, with the array TAUB, represent the unitary
   83: *          matrix Z as a product of elementary reflectors (see Further
   84: *          Details).
   85: *
   86: *  LDB     (input) INTEGER
   87: *          The leading dimension of the array B. LDB >= max(1,N).
   88: *
   89: *  TAUB    (output) COMPLEX*16 array, dimension (min(N,P))
   90: *          The scalar factors of the elementary reflectors which
   91: *          represent the unitary matrix Z (see Further Details).
   92: *
   93: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
   94: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   95: *
   96: *  LWORK   (input) INTEGER
   97: *          The dimension of the array WORK. LWORK >= max(1,N,M,P).
   98: *          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
   99: *          where NB1 is the optimal blocksize for the QR factorization
  100: *          of an N-by-M matrix, NB2 is the optimal blocksize for the
  101: *          RQ factorization of an N-by-P matrix, and NB3 is the optimal
  102: *          blocksize for a call of ZUNMQR.
  103: *
  104: *          If LWORK = -1, then a workspace query is assumed; the routine
  105: *          only calculates the optimal size of the WORK array, returns
  106: *          this value as the first entry of the WORK array, and no error
  107: *          message related to LWORK is issued by XERBLA.
  108: *
  109: *  INFO    (output) INTEGER
  110: *           = 0:  successful exit
  111: *           < 0:  if INFO = -i, the i-th argument had an illegal value.
  112: *
  113: *  Further Details
  114: *  ===============
  115: *
  116: *  The matrix Q is represented as a product of elementary reflectors
  117: *
  118: *     Q = H(1) H(2) . . . H(k), where k = min(n,m).
  119: *
  120: *  Each H(i) has the form
  121: *
  122: *     H(i) = I - taua * v * v'
  123: *
  124: *  where taua is a complex scalar, and v is a complex vector with
  125: *  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
  126: *  and taua in TAUA(i).
  127: *  To form Q explicitly, use LAPACK subroutine ZUNGQR.
  128: *  To use Q to update another matrix, use LAPACK subroutine ZUNMQR.
  129: *
  130: *  The matrix Z is represented as a product of elementary reflectors
  131: *
  132: *     Z = H(1) H(2) . . . H(k), where k = min(n,p).
  133: *
  134: *  Each H(i) has the form
  135: *
  136: *     H(i) = I - taub * v * v'
  137: *
  138: *  where taub is a complex scalar, and v is a complex vector with
  139: *  v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
  140: *  B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
  141: *  To form Z explicitly, use LAPACK subroutine ZUNGRQ.
  142: *  To use Z to update another matrix, use LAPACK subroutine ZUNMRQ.
  143: *
  144: *  =====================================================================
  145: *
  146: *     .. Local Scalars ..
  147:       LOGICAL            LQUERY
  148:       INTEGER            LOPT, LWKOPT, NB, NB1, NB2, NB3
  149: *     ..
  150: *     .. External Subroutines ..
  151:       EXTERNAL           XERBLA, ZGEQRF, ZGERQF, ZUNMQR
  152: *     ..
  153: *     .. External Functions ..
  154:       INTEGER            ILAENV
  155:       EXTERNAL           ILAENV
  156: *     ..
  157: *     .. Intrinsic Functions ..
  158:       INTRINSIC          INT, MAX, MIN
  159: *     ..
  160: *     .. Executable Statements ..
  161: *
  162: *     Test the input parameters
  163: *
  164:       INFO = 0
  165:       NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 )
  166:       NB2 = ILAENV( 1, 'ZGERQF', ' ', N, P, -1, -1 )
  167:       NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 )
  168:       NB = MAX( NB1, NB2, NB3 )
  169:       LWKOPT = MAX( N, M, P )*NB
  170:       WORK( 1 ) = LWKOPT
  171:       LQUERY = ( LWORK.EQ.-1 )
  172:       IF( N.LT.0 ) THEN
  173:          INFO = -1
  174:       ELSE IF( M.LT.0 ) THEN
  175:          INFO = -2
  176:       ELSE IF( P.LT.0 ) THEN
  177:          INFO = -3
  178:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  179:          INFO = -5
  180:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  181:          INFO = -8
  182:       ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN
  183:          INFO = -11
  184:       END IF
  185:       IF( INFO.NE.0 ) THEN
  186:          CALL XERBLA( 'ZGGQRF', -INFO )
  187:          RETURN
  188:       ELSE IF( LQUERY ) THEN
  189:          RETURN
  190:       END IF
  191: *
  192: *     QR factorization of N-by-M matrix A: A = Q*R
  193: *
  194:       CALL ZGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
  195:       LOPT = WORK( 1 )
  196: *
  197: *     Update B := Q'*B.
  198: *
  199:       CALL ZUNMQR( 'Left', 'Conjugate Transpose', N, P, MIN( N, M ), A,
  200:      $             LDA, TAUA, B, LDB, WORK, LWORK, INFO )
  201:       LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
  202: *
  203: *     RQ factorization of N-by-P matrix B: B = T*Z.
  204: *
  205:       CALL ZGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
  206:       WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
  207: *
  208:       RETURN
  209: *
  210: *     End of ZGGQRF
  211: *
  212:       END