001:       SUBROUTINE DGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
002:      $                   LWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
010: *     ..
011: *     .. Array Arguments ..
012:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
013:      $                   WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DGGQRF computes a generalized QR factorization of an N-by-M matrix A
020: *  and an N-by-P matrix B:
021: *
022: *              A = Q*R,        B = Q*T*Z,
023: *
024: *  where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
025: *  matrix, and R and T assume one of the forms:
026: *
027: *  if N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N,
028: *                  (  0  ) N-M                         N   M-N
029: *                     M
030: *
031: *  where R11 is upper triangular, and
032: *
033: *  if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N-P,
034: *                   P-N  N                           ( T21 ) P
035: *                                                       P
036: *
037: *  where T12 or T21 is upper triangular.
038: *
039: *  In particular, if B is square and nonsingular, the GQR factorization
040: *  of A and B implicitly gives the QR factorization of inv(B)*A:
041: *
042: *               inv(B)*A = Z'*(inv(T)*R)
043: *
044: *  where inv(B) denotes the inverse of the matrix B, and Z' denotes the
045: *  transpose of the matrix Z.
046: *
047: *  Arguments
048: *  =========
049: *
050: *  N       (input) INTEGER
051: *          The number of rows of the matrices A and B. N >= 0.
052: *
053: *  M       (input) INTEGER
054: *          The number of columns of the matrix A.  M >= 0.
055: *
056: *  P       (input) INTEGER
057: *          The number of columns of the matrix B.  P >= 0.
058: *
059: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,M)
060: *          On entry, the N-by-M matrix A.
061: *          On exit, the elements on and above the diagonal of the array
062: *          contain the min(N,M)-by-M upper trapezoidal matrix R (R is
063: *          upper triangular if N >= M); the elements below the diagonal,
064: *          with the array TAUA, represent the orthogonal matrix Q as a
065: *          product of min(N,M) elementary reflectors (see Further
066: *          Details).
067: *
068: *  LDA     (input) INTEGER
069: *          The leading dimension of the array A. LDA >= max(1,N).
070: *
071: *  TAUA    (output) DOUBLE PRECISION array, dimension (min(N,M))
072: *          The scalar factors of the elementary reflectors which
073: *          represent the orthogonal matrix Q (see Further Details).
074: *
075: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,P)
076: *          On entry, the N-by-P matrix B.
077: *          On exit, if N <= P, the upper triangle of the subarray
078: *          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
079: *          if N > P, the elements on and above the (N-P)-th subdiagonal
080: *          contain the N-by-P upper trapezoidal matrix T; the remaining
081: *          elements, with the array TAUB, represent the orthogonal
082: *          matrix Z as a product of elementary reflectors (see Further
083: *          Details).
084: *
085: *  LDB     (input) INTEGER
086: *          The leading dimension of the array B. LDB >= max(1,N).
087: *
088: *  TAUB    (output) DOUBLE PRECISION array, dimension (min(N,P))
089: *          The scalar factors of the elementary reflectors which
090: *          represent the orthogonal matrix Z (see Further Details).
091: *
092: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
093: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
094: *
095: *  LWORK   (input) INTEGER
096: *          The dimension of the array WORK. LWORK >= max(1,N,M,P).
097: *          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
098: *          where NB1 is the optimal blocksize for the QR factorization
099: *          of an N-by-M matrix, NB2 is the optimal blocksize for the
100: *          RQ factorization of an N-by-P matrix, and NB3 is the optimal
101: *          blocksize for a call of DORMQR.
102: *
103: *          If LWORK = -1, then a workspace query is assumed; the routine
104: *          only calculates the optimal size of the WORK array, returns
105: *          this value as the first entry of the WORK array, and no error
106: *          message related to LWORK is issued by XERBLA.
107: *
108: *  INFO    (output) INTEGER
109: *          = 0:  successful exit
110: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
111: *
112: *  Further Details
113: *  ===============
114: *
115: *  The matrix Q is represented as a product of elementary reflectors
116: *
117: *     Q = H(1) H(2) . . . H(k), where k = min(n,m).
118: *
119: *  Each H(i) has the form
120: *
121: *     H(i) = I - taua * v * v'
122: *
123: *  where taua is a real scalar, and v is a real vector with
124: *  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
125: *  and taua in TAUA(i).
126: *  To form Q explicitly, use LAPACK subroutine DORGQR.
127: *  To use Q to update another matrix, use LAPACK subroutine DORMQR.
128: *
129: *  The matrix Z is represented as a product of elementary reflectors
130: *
131: *     Z = H(1) H(2) . . . H(k), where k = min(n,p).
132: *
133: *  Each H(i) has the form
134: *
135: *     H(i) = I - taub * v * v'
136: *
137: *  where taub is a real scalar, and v is a real vector with
138: *  v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
139: *  B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
140: *  To form Z explicitly, use LAPACK subroutine DORGRQ.
141: *  To use Z to update another matrix, use LAPACK subroutine DORMRQ.
142: *
143: *  =====================================================================
144: *
145: *     .. Local Scalars ..
146:       LOGICAL            LQUERY
147:       INTEGER            LOPT, LWKOPT, NB, NB1, NB2, NB3
148: *     ..
149: *     .. External Subroutines ..
150:       EXTERNAL           DGEQRF, DGERQF, DORMQR, XERBLA
151: *     ..
152: *     .. External Functions ..
153:       INTEGER            ILAENV
154:       EXTERNAL           ILAENV
155: *     ..
156: *     .. Intrinsic Functions ..
157:       INTRINSIC          INT, MAX, MIN
158: *     ..
159: *     .. Executable Statements ..
160: *
161: *     Test the input parameters
162: *
163:       INFO = 0
164:       NB1 = ILAENV( 1, 'DGEQRF', ' ', N, M, -1, -1 )
165:       NB2 = ILAENV( 1, 'DGERQF', ' ', N, P, -1, -1 )
166:       NB3 = ILAENV( 1, 'DORMQR', ' ', N, M, P, -1 )
167:       NB = MAX( NB1, NB2, NB3 )
168:       LWKOPT = MAX( N, M, P )*NB
169:       WORK( 1 ) = LWKOPT
170:       LQUERY = ( LWORK.EQ.-1 )
171:       IF( N.LT.0 ) THEN
172:          INFO = -1
173:       ELSE IF( M.LT.0 ) THEN
174:          INFO = -2
175:       ELSE IF( P.LT.0 ) THEN
176:          INFO = -3
177:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
178:          INFO = -5
179:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
180:          INFO = -8
181:       ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN
182:          INFO = -11
183:       END IF
184:       IF( INFO.NE.0 ) THEN
185:          CALL XERBLA( 'DGGQRF', -INFO )
186:          RETURN
187:       ELSE IF( LQUERY ) THEN
188:          RETURN
189:       END IF
190: *
191: *     QR factorization of N-by-M matrix A: A = Q*R
192: *
193:       CALL DGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
194:       LOPT = WORK( 1 )
195: *
196: *     Update B := Q'*B.
197: *
198:       CALL DORMQR( 'Left', 'Transpose', N, P, MIN( N, M ), A, LDA, TAUA,
199:      $             B, LDB, WORK, LWORK, INFO )
200:       LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
201: *
202: *     RQ factorization of N-by-P matrix B: B = T*Z.
203: *
204:       CALL DGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
205:       WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
206: *
207:       RETURN
208: *
209: *     End of DGGQRF
210: *
211:       END
212: