1*> \brief <b> DGGLSE solves overdetermined or underdetermined systems for OTHER matrices</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGGLSE + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgglse.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgglse.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgglse.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
22*                          INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
26*       ..
27*       .. Array Arguments ..
28*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( * ), D( * ),
29*      $                   WORK( * ), X( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DGGLSE solves the linear equality-constrained least squares (LSE)
39*> problem:
40*>
41*>         minimize || c - A*x ||_2   subject to   B*x = d
42*>
43*> where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
44*> M-vector, and d is a given P-vector. It is assumed that
45*> P <= N <= M+P, and
46*>
47*>          rank(B) = P and  rank( (A) ) = N.
48*>                               ( (B) )
49*>
50*> These conditions ensure that the LSE problem has a unique solution,
51*> which is obtained using a generalized RQ factorization of the
52*> matrices (B, A) given by
53*>
54*>    B = (0 R)*Q,   A = Z*T*Q.
55*> \endverbatim
56*
57*  Arguments:
58*  ==========
59*
60*> \param[in] M
61*> \verbatim
62*>          M is INTEGER
63*>          The number of rows of the matrix A.  M >= 0.
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*>          N is INTEGER
69*>          The number of columns of the matrices A and B. N >= 0.
70*> \endverbatim
71*>
72*> \param[in] P
73*> \verbatim
74*>          P is INTEGER
75*>          The number of rows of the matrix B. 0 <= P <= N <= M+P.
76*> \endverbatim
77*>
78*> \param[in,out] A
79*> \verbatim
80*>          A is DOUBLE PRECISION array, dimension (LDA,N)
81*>          On entry, the M-by-N matrix A.
82*>          On exit, the elements on and above the diagonal of the array
83*>          contain the min(M,N)-by-N upper trapezoidal matrix T.
84*> \endverbatim
85*>
86*> \param[in] LDA
87*> \verbatim
88*>          LDA is INTEGER
89*>          The leading dimension of the array A. LDA >= max(1,M).
90*> \endverbatim
91*>
92*> \param[in,out] B
93*> \verbatim
94*>          B is DOUBLE PRECISION array, dimension (LDB,N)
95*>          On entry, the P-by-N matrix B.
96*>          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
97*>          contains the P-by-P upper triangular matrix R.
98*> \endverbatim
99*>
100*> \param[in] LDB
101*> \verbatim
102*>          LDB is INTEGER
103*>          The leading dimension of the array B. LDB >= max(1,P).
104*> \endverbatim
105*>
106*> \param[in,out] C
107*> \verbatim
108*>          C is DOUBLE PRECISION array, dimension (M)
109*>          On entry, C contains the right hand side vector for the
110*>          least squares part of the LSE problem.
111*>          On exit, the residual sum of squares for the solution
112*>          is given by the sum of squares of elements N-P+1 to M of
113*>          vector C.
114*> \endverbatim
115*>
116*> \param[in,out] D
117*> \verbatim
118*>          D is DOUBLE PRECISION array, dimension (P)
119*>          On entry, D contains the right hand side vector for the
120*>          constrained equation.
121*>          On exit, D is destroyed.
122*> \endverbatim
123*>
124*> \param[out] X
125*> \verbatim
126*>          X is DOUBLE PRECISION array, dimension (N)
127*>          On exit, X is the solution of the LSE problem.
128*> \endverbatim
129*>
130*> \param[out] WORK
131*> \verbatim
132*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
133*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
134*> \endverbatim
135*>
136*> \param[in] LWORK
137*> \verbatim
138*>          LWORK is INTEGER
139*>          The dimension of the array WORK. LWORK >= max(1,M+N+P).
140*>          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
141*>          where NB is an upper bound for the optimal blocksizes for
142*>          DGEQRF, SGERQF, DORMQR and SORMRQ.
143*>
144*>          If LWORK = -1, then a workspace query is assumed; the routine
145*>          only calculates the optimal size of the WORK array, returns
146*>          this value as the first entry of the WORK array, and no error
147*>          message related to LWORK is issued by XERBLA.
148*> \endverbatim
149*>
150*> \param[out] INFO
151*> \verbatim
152*>          INFO is INTEGER
153*>          = 0:  successful exit.
154*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
155*>          = 1:  the upper triangular factor R associated with B in the
156*>                generalized RQ factorization of the pair (B, A) is
157*>                singular, so that rank(B) < P; the least squares
158*>                solution could not be computed.
159*>          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
160*>                T associated with A in the generalized RQ factorization
161*>                of the pair (B, A) is singular, so that
162*>                rank( (A) ) < N; the least squares solution could not
163*>                    ( (B) )
164*>                be computed.
165*> \endverbatim
166*
167*  Authors:
168*  ========
169*
170*> \author Univ. of Tennessee
171*> \author Univ. of California Berkeley
172*> \author Univ. of Colorado Denver
173*> \author NAG Ltd.
174*
175*> \ingroup doubleOTHERsolve
176*
177*  =====================================================================
178      SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
179     $                   INFO )
180*
181*  -- LAPACK driver routine --
182*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
183*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184*
185*     .. Scalar Arguments ..
186      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
187*     ..
188*     .. Array Arguments ..
189      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190     $                   WORK( * ), X( * )
191*     ..
192*
193*  =====================================================================
194*
195*     .. Parameters ..
196      DOUBLE PRECISION   ONE
197      PARAMETER          ( ONE = 1.0D+0 )
198*     ..
199*     .. Local Scalars ..
200      LOGICAL            LQUERY
201      INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
202     $                   NB4, NR
203*     ..
204*     .. External Subroutines ..
205      EXTERNAL           DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ,
206     $                   DTRMV, DTRTRS, XERBLA
207*     ..
208*     .. External Functions ..
209      INTEGER            ILAENV
210      EXTERNAL           ILAENV
211*     ..
212*     .. Intrinsic Functions ..
213      INTRINSIC          INT, MAX, MIN
214*     ..
215*     .. Executable Statements ..
216*
217*     Test the input parameters
218*
219      INFO = 0
220      MN = MIN( M, N )
221      LQUERY = ( LWORK.EQ.-1 )
222      IF( M.LT.0 ) THEN
223         INFO = -1
224      ELSE IF( N.LT.0 ) THEN
225         INFO = -2
226      ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
227         INFO = -3
228      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
229         INFO = -5
230      ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
231         INFO = -7
232      END IF
233*
234*     Calculate workspace
235*
236      IF( INFO.EQ.0) THEN
237         IF( N.EQ.0 ) THEN
238            LWKMIN = 1
239            LWKOPT = 1
240         ELSE
241            NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
242            NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
243            NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 )
244            NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, P, -1 )
245            NB = MAX( NB1, NB2, NB3, NB4 )
246            LWKMIN = M + N + P
247            LWKOPT = P + MN + MAX( M, N )*NB
248         END IF
249         WORK( 1 ) = LWKOPT
250*
251         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
252            INFO = -12
253         END IF
254      END IF
255*
256      IF( INFO.NE.0 ) THEN
257         CALL XERBLA( 'DGGLSE', -INFO )
258         RETURN
259      ELSE IF( LQUERY ) THEN
260         RETURN
261      END IF
262*
263*     Quick return if possible
264*
265      IF( N.EQ.0 )
266     $   RETURN
267*
268*     Compute the GRQ factorization of matrices B and A:
269*
270*            B*Q**T = (  0  T12 ) P   Z**T*A*Q**T = ( R11 R12 ) N-P
271*                        N-P  P                     (  0  R22 ) M+P-N
272*                                                      N-P  P
273*
274*     where T12 and R11 are upper triangular, and Q and Z are
275*     orthogonal.
276*
277      CALL DGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
278     $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
279      LOPT = WORK( P+MN+1 )
280*
281*     Update c = Z**T *c = ( c1 ) N-P
282*                          ( c2 ) M+P-N
283*
284      CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
285     $             C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
286      LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
287*
288*     Solve T12*x2 = d for x2
289*
290      IF( P.GT.0 ) THEN
291         CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
292     $                B( 1, N-P+1 ), LDB, D, P, INFO )
293*
294         IF( INFO.GT.0 ) THEN
295            INFO = 1
296            RETURN
297         END IF
298*
299*        Put the solution in X
300*
301         CALL DCOPY( P, D, 1, X( N-P+1 ), 1 )
302*
303*        Update c1
304*
305         CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
306     $               D, 1, ONE, C, 1 )
307      END IF
308*
309*     Solve R11*x1 = c1 for x1
310*
311      IF( N.GT.P ) THEN
312         CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
313     $                A, LDA, C, N-P, INFO )
314*
315         IF( INFO.GT.0 ) THEN
316            INFO = 2
317            RETURN
318         END IF
319*
320*        Put the solutions in X
321*
322         CALL DCOPY( N-P, C, 1, X, 1 )
323      END IF
324*
325*     Compute the residual vector:
326*
327      IF( M.LT.N ) THEN
328         NR = M + P - N
329         IF( NR.GT.0 )
330     $      CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
331     $                  LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
332      ELSE
333         NR = P
334      END IF
335      IF( NR.GT.0 ) THEN
336         CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR,
337     $               A( N-P+1, N-P+1 ), LDA, D, 1 )
338         CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
339      END IF
340*
341*     Backward transformation x = Q**T*x
342*
343      CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
344     $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
345      WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
346*
347      RETURN
348*
349*     End of DGGLSE
350*
351      END
352