LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
cgetsls.f
Go to the documentation of this file.
1 *> \brief \b CGETSLS
2 *
3 * Definition:
4 * ===========
5 *
6 * SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
7 * $ WORK, LWORK, INFO )
8 *
9 * .. Scalar Arguments ..
10 * CHARACTER TRANS
11 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
12 * ..
13 * .. Array Arguments ..
14 * COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
15 * ..
16 *
17 *
18 *> \par Purpose:
19 * =============
20 *>
21 *> \verbatim
22 *>
23 *> CGETSLS solves overdetermined or underdetermined complex linear systems
24 *> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
25 *> factorization of A. It is assumed that A has full rank.
26 *>
27 *>
28 *>
29 *> The following options are provided:
30 *>
31 *> 1. If TRANS = 'N' and m >= n: find the least squares solution of
32 *> an overdetermined system, i.e., solve the least squares problem
33 *> minimize || B - A*X ||.
34 *>
35 *> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
36 *> an underdetermined system A * X = B.
37 *>
38 *> 3. If TRANS = 'C' and m >= n: find the minimum norm solution of
39 *> an undetermined system A**T * X = B.
40 *>
41 *> 4. If TRANS = 'C' and m < n: find the least squares solution of
42 *> an overdetermined system, i.e., solve the least squares problem
43 *> minimize || B - A**T * X ||.
44 *>
45 *> Several right hand side vectors b and solution vectors x can be
46 *> handled in a single call; they are stored as the columns of the
47 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
48 *> matrix X.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[in] TRANS
55 *> \verbatim
56 *> TRANS is CHARACTER*1
57 *> = 'N': the linear system involves A;
58 *> = 'C': the linear system involves A**H.
59 *> \endverbatim
60 *>
61 *> \param[in] M
62 *> \verbatim
63 *> M is INTEGER
64 *> The number of rows of the matrix A. M >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] N
68 *> \verbatim
69 *> N is INTEGER
70 *> The number of columns of the matrix A. N >= 0.
71 *> \endverbatim
72 *>
73 *> \param[in] NRHS
74 *> \verbatim
75 *> NRHS is INTEGER
76 *> The number of right hand sides, i.e., the number of
77 *> columns of the matrices B and X. NRHS >=0.
78 *> \endverbatim
79 *>
80 *> \param[in,out] A
81 *> \verbatim
82 *> A is COMPLEX array, dimension (LDA,N)
83 *> On entry, the M-by-N matrix A.
84 *> On exit,
85 *> A is overwritten by details of its QR or LQ
86 *> factorization as returned by CGEQR or CGELQ.
87 *> \endverbatim
88 *>
89 *> \param[in] LDA
90 *> \verbatim
91 *> LDA is INTEGER
92 *> The leading dimension of the array A. LDA >= max(1,M).
93 *> \endverbatim
94 *>
95 *> \param[in,out] B
96 *> \verbatim
97 *> B is COMPLEX array, dimension (LDB,NRHS)
98 *> On entry, the matrix B of right hand side vectors, stored
99 *> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
100 *> if TRANS = 'C'.
101 *> On exit, if INFO = 0, B is overwritten by the solution
102 *> vectors, stored columnwise:
103 *> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
104 *> squares solution vectors.
105 *> if TRANS = 'N' and m < n, rows 1 to N of B contain the
106 *> minimum norm solution vectors;
107 *> if TRANS = 'C' and m >= n, rows 1 to M of B contain the
108 *> minimum norm solution vectors;
109 *> if TRANS = 'C' and m < n, rows 1 to M of B contain the
110 *> least squares solution vectors.
111 *> \endverbatim
112 *>
113 *> \param[in] LDB
114 *> \verbatim
115 *> LDB is INTEGER
116 *> The leading dimension of the array B. LDB >= MAX(1,M,N).
117 *> \endverbatim
118 *>
119 *> \param[out] WORK
120 *> \verbatim
121 *> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
122 *> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
123 *> or optimal, if query was assumed) LWORK.
124 *> See LWORK for details.
125 *> \endverbatim
126 *>
127 *> \param[in] LWORK
128 *> \verbatim
129 *> LWORK is INTEGER
130 *> The dimension of the array WORK.
131 *> If LWORK = -1 or -2, then a workspace query is assumed.
132 *> If LWORK = -1, the routine calculates optimal size of WORK for the
133 *> optimal performance and returns this value in WORK(1).
134 *> If LWORK = -2, the routine calculates minimal size of WORK and
135 *> returns this value in WORK(1).
136 *> \endverbatim
137 *>
138 *> \param[out] INFO
139 *> \verbatim
140 *> INFO is INTEGER
141 *> = 0: successful exit
142 *> < 0: if INFO = -i, the i-th argument had an illegal value
143 *> > 0: if INFO = i, the i-th diagonal element of the
144 *> triangular factor of A is zero, so that A does not have
145 *> full rank; the least squares solution could not be
146 *> computed.
147 *> \endverbatim
148 *
149 * Authors:
150 * ========
151 *
152 *> \author Univ. of Tennessee
153 *> \author Univ. of California Berkeley
154 *> \author Univ. of Colorado Denver
155 *> \author NAG Ltd.
156 *
157 *> \date June 2017
158 *
159 *> \ingroup complexGEsolve
160 *
161 * =====================================================================
162  SUBROUTINE cgetsls( TRANS, M, N, NRHS, A, LDA, B, LDB,
163  $ WORK, LWORK, INFO )
164 *
165 * -- LAPACK driver routine (version 3.7.1) --
166 * -- LAPACK is a software package provided by Univ. of Tennessee, --
167 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168 * June 2017
169 *
170 * .. Scalar Arguments ..
171  CHARACTER TRANS
172  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
173 * ..
174 * .. Array Arguments ..
175  COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
176 *
177 * ..
178 *
179 * =====================================================================
180 *
181 * .. Parameters ..
182  REAL ZERO, ONE
183  parameter( zero = 0.0e0, one = 1.0e0 )
184  COMPLEX CZERO
185  parameter( czero = ( 0.0e+0, 0.0e+0 ) )
186 * ..
187 * .. Local Scalars ..
188  LOGICAL LQUERY, TRAN
189  INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
190  $ scllen, mnk, tszo, tszm, lwo, lwm, lw1, lw2,
191  $ wsizeo, wsizem, info2
192  REAL ANRM, BIGNUM, BNRM, SMLNUM, DUM( 1 )
193  COMPLEX TQ( 5 ), WORKQ( 1 )
194 * ..
195 * .. External Functions ..
196  LOGICAL LSAME
197  INTEGER ILAENV
198  REAL SLAMCH, CLANGE
199  EXTERNAL lsame, ilaenv, slabad, slamch, clange
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL cgeqr, cgemqr, clascl, claset,
203  $ ctrtrs, xerbla, cgelq, cgemlq
204 * ..
205 * .. Intrinsic Functions ..
206  INTRINSIC real, max, min, int
207 * ..
208 * .. Executable Statements ..
209 *
210 * Test the input arguments.
211 *
212  info = 0
213  minmn = min( m, n )
214  maxmn = max( m, n )
215  mnk = max( minmn, nrhs )
216  tran = lsame( trans, 'C' )
217 *
218  lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
219  IF( .NOT.( lsame( trans, 'N' ) .OR.
220  $ lsame( trans, 'C' ) ) ) THEN
221  info = -1
222  ELSE IF( m.LT.0 ) THEN
223  info = -2
224  ELSE IF( n.LT.0 ) THEN
225  info = -3
226  ELSE IF( nrhs.LT.0 ) THEN
227  info = -4
228  ELSE IF( lda.LT.max( 1, m ) ) THEN
229  info = -6
230  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
231  info = -8
232  END IF
233 *
234  IF( info.EQ.0 ) THEN
235 *
236 * Determine the block size and minimum LWORK
237 *
238  IF( m.GE.n ) THEN
239  CALL cgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
240  tszo = int( tq( 1 ) )
241  lwo = int( workq( 1 ) )
242  CALL cgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
243  $ tszo, b, ldb, workq, -1, info2 )
244  lwo = max( lwo, int( workq( 1 ) ) )
245  CALL cgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
246  tszm = int( tq( 1 ) )
247  lwm = int( workq( 1 ) )
248  CALL cgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
249  $ tszm, b, ldb, workq, -1, info2 )
250  lwm = max( lwm, int( workq( 1 ) ) )
251  wsizeo = tszo + lwo
252  wsizem = tszm + lwm
253  ELSE
254  CALL cgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
255  tszo = int( tq( 1 ) )
256  lwo = int( workq( 1 ) )
257  CALL cgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
258  $ tszo, b, ldb, workq, -1, info2 )
259  lwo = max( lwo, int( workq( 1 ) ) )
260  CALL cgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
261  tszm = int( tq( 1 ) )
262  lwm = int( workq( 1 ) )
263  CALL cgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
264  $ tszo, b, ldb, workq, -1, info2 )
265  lwm = max( lwm, int( workq( 1 ) ) )
266  wsizeo = tszo + lwo
267  wsizem = tszm + lwm
268  END IF
269 *
270  IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
271  info = -10
272  END IF
273 *
274  END IF
275 *
276  IF( info.NE.0 ) THEN
277  CALL xerbla( 'CGETSLS', -info )
278  work( 1 ) = real( wsizeo )
279  RETURN
280  END IF
281  IF( lquery ) THEN
282  IF( lwork.EQ.-1 ) work( 1 ) = real( wsizeo )
283  IF( lwork.EQ.-2 ) work( 1 ) = real( wsizem )
284  RETURN
285  END IF
286  IF( lwork.LT.wsizeo ) THEN
287  lw1 = tszm
288  lw2 = lwm
289  ELSE
290  lw1 = tszo
291  lw2 = lwo
292  END IF
293 *
294 * Quick return if possible
295 *
296  IF( min( m, n, nrhs ).EQ.0 ) THEN
297  CALL claset( 'FULL', max( m, n ), nrhs, czero, czero,
298  $ b, ldb )
299  RETURN
300  END IF
301 *
302 * Get machine parameters
303 *
304  smlnum = slamch( 'S' ) / slamch( 'P' )
305  bignum = one / smlnum
306  CALL slabad( smlnum, bignum )
307 *
308 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
309 *
310  anrm = clange( 'M', m, n, a, lda, dum )
311  iascl = 0
312  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
313 *
314 * Scale matrix norm up to SMLNUM
315 *
316  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
317  iascl = 1
318  ELSE IF( anrm.GT.bignum ) THEN
319 *
320 * Scale matrix norm down to BIGNUM
321 *
322  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
323  iascl = 2
324  ELSE IF( anrm.EQ.zero ) THEN
325 *
326 * Matrix all zero. Return zero solution.
327 *
328  CALL claset( 'F', maxmn, nrhs, czero, czero, b, ldb )
329  GO TO 50
330  END IF
331 *
332  brow = m
333  IF ( tran ) THEN
334  brow = n
335  END IF
336  bnrm = clange( 'M', brow, nrhs, b, ldb, dum )
337  ibscl = 0
338  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
339 *
340 * Scale matrix norm up to SMLNUM
341 *
342  CALL clascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
343  $ info )
344  ibscl = 1
345  ELSE IF( bnrm.GT.bignum ) THEN
346 *
347 * Scale matrix norm down to BIGNUM
348 *
349  CALL clascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
350  $ info )
351  ibscl = 2
352  END IF
353 *
354  IF ( m.GE.n ) THEN
355 *
356 * compute QR factorization of A
357 *
358  CALL cgeqr( m, n, a, lda, work( lw2+1 ), lw1,
359  $ work( 1 ), lw2, info )
360  IF ( .NOT.tran ) THEN
361 *
362 * Least-Squares Problem min || A * X - B ||
363 *
364 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
365 *
366  CALL cgemqr( 'L' , 'C', m, nrhs, n, a, lda,
367  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
368  $ info )
369 *
370 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
371 *
372  CALL ctrtrs( 'U', 'N', 'N', n, nrhs,
373  $ a, lda, b, ldb, info )
374  IF( info.GT.0 ) THEN
375  RETURN
376  END IF
377  scllen = n
378  ELSE
379 *
380 * Overdetermined system of equations A**T * X = B
381 *
382 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
383 *
384  CALL ctrtrs( 'U', 'C', 'N', n, nrhs,
385  $ a, lda, b, ldb, info )
386 *
387  IF( info.GT.0 ) THEN
388  RETURN
389  END IF
390 *
391 * B(N+1:M,1:NRHS) = CZERO
392 *
393  DO 20 j = 1, nrhs
394  DO 10 i = n + 1, m
395  b( i, j ) = czero
396  10 CONTINUE
397  20 CONTINUE
398 *
399 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
400 *
401  CALL cgemqr( 'L', 'N', m, nrhs, n, a, lda,
402  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
403  $ info )
404 *
405  scllen = m
406 *
407  END IF
408 *
409  ELSE
410 *
411 * Compute LQ factorization of A
412 *
413  CALL cgelq( m, n, a, lda, work( lw2+1 ), lw1,
414  $ work( 1 ), lw2, info )
415 *
416 * workspace at least M, optimally M*NB.
417 *
418  IF( .NOT.tran ) THEN
419 *
420 * underdetermined system of equations A * X = B
421 *
422 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
423 *
424  CALL ctrtrs( 'L', 'N', 'N', m, nrhs,
425  $ a, lda, b, ldb, info )
426 *
427  IF( info.GT.0 ) THEN
428  RETURN
429  END IF
430 *
431 * B(M+1:N,1:NRHS) = 0
432 *
433  DO 40 j = 1, nrhs
434  DO 30 i = m + 1, n
435  b( i, j ) = czero
436  30 CONTINUE
437  40 CONTINUE
438 *
439 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
440 *
441  CALL cgemlq( 'L', 'C', n, nrhs, m, a, lda,
442  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
443  $ info )
444 *
445 * workspace at least NRHS, optimally NRHS*NB
446 *
447  scllen = n
448 *
449  ELSE
450 *
451 * overdetermined system min || A**T * X - B ||
452 *
453 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
454 *
455  CALL cgemlq( 'L', 'N', n, nrhs, m, a, lda,
456  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
457  $ info )
458 *
459 * workspace at least NRHS, optimally NRHS*NB
460 *
461 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
462 *
463  CALL ctrtrs( 'L', 'C', 'N', m, nrhs,
464  $ a, lda, b, ldb, info )
465 *
466  IF( info.GT.0 ) THEN
467  RETURN
468  END IF
469 *
470  scllen = m
471 *
472  END IF
473 *
474  END IF
475 *
476 * Undo scaling
477 *
478  IF( iascl.EQ.1 ) THEN
479  CALL clascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
480  $ info )
481  ELSE IF( iascl.EQ.2 ) THEN
482  CALL clascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
483  $ info )
484  END IF
485  IF( ibscl.EQ.1 ) THEN
486  CALL clascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
487  $ info )
488  ELSE IF( ibscl.EQ.2 ) THEN
489  CALL clascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
490  $ info )
491  END IF
492 *
493  50 CONTINUE
494  work( 1 ) = real( tszo + lwo )
495  RETURN
496 *
497 * End of ZGETSLS
498 *
499  END
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
ctrtrs
subroutine ctrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
CTRTRS
Definition: ctrtrs.f:142
cgelq
subroutine cgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
CGELQ
Definition: cgelq.f:172
cgeqr
subroutine cgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
CGEQR
Definition: cgeqr.f:174
cgemlq
subroutine cgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
CGEMLQ
Definition: cgemlq.f:170
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
cgemqr
subroutine cgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
CGEMQR
Definition: cgemqr.f:172
clascl
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145
cgetsls
subroutine cgetsls(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
CGETSLS
Definition: cgetsls.f:164