001:       SUBROUTINE SGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
002:      $                   LWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
014:      $                   WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SGGQRF computes a generalized QR factorization of an N-by-M matrix A
021: *  and an N-by-P matrix B:
022: *
023: *              A = Q*R,        B = Q*T*Z,
024: *
025: *  where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
026: *  matrix, and R and T assume one of the forms:
027: *
028: *  if N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N,
029: *                  (  0  ) N-M                         N   M-N
030: *                     M
031: *
032: *  where R11 is upper triangular, and
033: *
034: *  if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N-P,
035: *                   P-N  N                           ( T21 ) P
036: *                                                       P
037: *
038: *  where T12 or T21 is upper triangular.
039: *
040: *  In particular, if B is square and nonsingular, the GQR factorization
041: *  of A and B implicitly gives the QR factorization of inv(B)*A:
042: *
043: *               inv(B)*A = Z'*(inv(T)*R)
044: *
045: *  where inv(B) denotes the inverse of the matrix B, and Z' denotes the
046: *  transpose of the matrix Z.
047: *
048: *  Arguments
049: *  =========
050: *
051: *  N       (input) INTEGER
052: *          The number of rows of the matrices A and B. N >= 0.
053: *
054: *  M       (input) INTEGER
055: *          The number of columns of the matrix A.  M >= 0.
056: *
057: *  P       (input) INTEGER
058: *          The number of columns of the matrix B.  P >= 0.
059: *
060: *  A       (input/output) REAL array, dimension (LDA,M)
061: *          On entry, the N-by-M matrix A.
062: *          On exit, the elements on and above the diagonal of the array
063: *          contain the min(N,M)-by-M upper trapezoidal matrix R (R is
064: *          upper triangular if N >= M); the elements below the diagonal,
065: *          with the array TAUA, represent the orthogonal matrix Q as a
066: *          product of min(N,M) elementary reflectors (see Further
067: *          Details).
068: *
069: *  LDA     (input) INTEGER
070: *          The leading dimension of the array A. LDA >= max(1,N).
071: *
072: *  TAUA    (output) REAL array, dimension (min(N,M))
073: *          The scalar factors of the elementary reflectors which
074: *          represent the orthogonal matrix Q (see Further Details).
075: *
076: *  B       (input/output) REAL array, dimension (LDB,P)
077: *          On entry, the N-by-P matrix B.
078: *          On exit, if N <= P, the upper triangle of the subarray
079: *          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
080: *          if N > P, the elements on and above the (N-P)-th subdiagonal
081: *          contain the N-by-P upper trapezoidal matrix T; the remaining
082: *          elements, with the array TAUB, represent the orthogonal
083: *          matrix Z as a product of elementary reflectors (see Further
084: *          Details).
085: *
086: *  LDB     (input) INTEGER
087: *          The leading dimension of the array B. LDB >= max(1,N).
088: *
089: *  TAUB    (output) REAL array, dimension (min(N,P))
090: *          The scalar factors of the elementary reflectors which
091: *          represent the orthogonal matrix Z (see Further Details).
092: *
093: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
094: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
095: *
096: *  LWORK   (input) INTEGER
097: *          The dimension of the array WORK. LWORK >= max(1,N,M,P).
098: *          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
099: *          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 SORMQR.
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 real scalar, and v is a real 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 SORGQR.
128: *  To use Q to update another matrix, use LAPACK subroutine SORMQR.
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 real scalar, and v is a real 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 SORGRQ.
142: *  To use Z to update another matrix, use LAPACK subroutine SORMRQ.
143: *
144: *  =====================================================================
145: *
146: *     .. Local Scalars ..
147:       LOGICAL            LQUERY
148:       INTEGER            LOPT, LWKOPT, NB, NB1, NB2, NB3
149: *     ..
150: *     .. External Subroutines ..
151:       EXTERNAL           SGEQRF, SGERQF, SORMQR, XERBLA
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, 'SGEQRF', ' ', N, M, -1, -1 )
166:       NB2 = ILAENV( 1, 'SGERQF', ' ', N, P, -1, -1 )
167:       NB3 = ILAENV( 1, 'SORMQR', ' ', 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( 'SGGQRF', -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 SGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
195:       LOPT = WORK( 1 )
196: *
197: *     Update B := Q'*B.
198: *
199:       CALL SORMQR( 'Left', 'Transpose', N, P, MIN( N, M ), A, LDA, TAUA,
200:      $             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 SGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
206:       WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
207: *
208:       RETURN
209: *
210: *     End of SGGQRF
211: *
212:       END
213: