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*> \date November 2011
176*
177*> \ingroup doubleOTHERsolve
178*
179*  =====================================================================
180      SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
181     $                   INFO )
182*
183*  -- LAPACK driver routine (version 3.4.0) --
184*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
185*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*     November 2011
187*
188*     .. Scalar Arguments ..
189      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
190*     ..
191*     .. Array Arguments ..
192      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( * ), D( * ),
193     $                   WORK( * ), X( * )
194*     ..
195*
196*  =====================================================================
197*
198*     .. Parameters ..
199      DOUBLE PRECISION   ONE
200      PARAMETER          ( ONE = 1.0D+0 )
201*     ..
202*     .. Local Scalars ..
203      LOGICAL            LQUERY
204      INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
205     $                   NB4, NR
206*     ..
207*     .. External Subroutines ..
208      EXTERNAL           DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ,
209     $                   DTRMV, DTRTRS, XERBLA
210*     ..
211*     .. External Functions ..
212      INTEGER            ILAENV
213      EXTERNAL           ILAENV
214*     ..
215*     .. Intrinsic Functions ..
216      INTRINSIC          INT, MAX, MIN
217*     ..
218*     .. Executable Statements ..
219*
220*     Test the input parameters
221*
222      INFO = 0
223      MN = MIN( M, N )
224      LQUERY = ( LWORK.EQ.-1 )
225      IF( M.LT.0 ) THEN
226         INFO = -1
227      ELSE IF( N.LT.0 ) THEN
228         INFO = -2
229      ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
230         INFO = -3
231      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
232         INFO = -5
233      ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
234         INFO = -7
235      END IF
236*
237*     Calculate workspace
238*
239      IF( INFO.EQ.0) THEN
240         IF( N.EQ.0 ) THEN
241            LWKMIN = 1
242            LWKOPT = 1
243         ELSE
244            NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
245            NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
246            NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 )
247            NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, P, -1 )
248            NB = MAX( NB1, NB2, NB3, NB4 )
249            LWKMIN = M + N + P
250            LWKOPT = P + MN + MAX( M, N )*NB
251         END IF
252         WORK( 1 ) = LWKOPT
253*
254         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
255            INFO = -12
256         END IF
257      END IF
258*
259      IF( INFO.NE.0 ) THEN
260         CALL XERBLA( 'DGGLSE', -INFO )
261         RETURN
262      ELSE IF( LQUERY ) THEN
263         RETURN
264      END IF
265*
266*     Quick return if possible
267*
268      IF( N.EQ.0 )
269     $   RETURN
270*
271*     Compute the GRQ factorization of matrices B and A:
272*
273*            B*Q**T = (  0  T12 ) P   Z**T*A*Q**T = ( R11 R12 ) N-P
274*                        N-P  P                     (  0  R22 ) M+P-N
275*                                                      N-P  P
276*
277*     where T12 and R11 are upper triangular, and Q and Z are
278*     orthogonal.
279*
280      CALL DGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
281     $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
282      LOPT = WORK( P+MN+1 )
283*
284*     Update c = Z**T *c = ( c1 ) N-P
285*                          ( c2 ) M+P-N
286*
287      CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
288     $             C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
289      LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
290*
291*     Solve T12*x2 = d for x2
292*
293      IF( P.GT.0 ) THEN
294         CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
295     $                B( 1, N-P+1 ), LDB, D, P, INFO )
296*
297         IF( INFO.GT.0 ) THEN
298            INFO = 1
299            RETURN
300         END IF
301*
302*        Put the solution in X
303*
304         CALL DCOPY( P, D, 1, X( N-P+1 ), 1 )
305*
306*        Update c1
307*
308         CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
309     $               D, 1, ONE, C, 1 )
310      END IF
311*
312*     Solve R11*x1 = c1 for x1
313*
314      IF( N.GT.P ) THEN
315         CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
316     $                A, LDA, C, N-P, INFO )
317*
318         IF( INFO.GT.0 ) THEN
319            INFO = 2
320            RETURN
321         END IF
322*
323*        Put the solutions in X
324*
325         CALL DCOPY( N-P, C, 1, X, 1 )
326      END IF
327*
328*     Compute the residual vector:
329*
330      IF( M.LT.N ) THEN
331         NR = M + P - N
332         IF( NR.GT.0 )
333     $      CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
334     $                  LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
335      ELSE
336         NR = P
337      END IF
338      IF( NR.GT.0 ) THEN
339         CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR,
340     $               A( N-P+1, N-P+1 ), LDA, D, 1 )
341         CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
342      END IF
343*
344*     Backward transformation x = Q**T*x
345*
346      CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
347     $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
348      WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
349*
350      RETURN
351*
352*     End of DGGLSE
353*
354      END
355