LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine cgetsls ( character  TRANS,
integer  M,
integer  N,
integer  NRHS,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)
Purpose:

CGETSLS solves overdetermined or underdetermined complex linear systems involving an M-by-N matrix A, using a tall skinny QR or short wide LQ factorization of A. It is assumed that A has full rank.

The following options are provided:

  1. If TRANS = 'N' and m >= n: find the least squares solution of an overdetermined system, i.e., solve the least squares problem minimize || B - A*X ||.
  2. If TRANS = 'N' and m < n: find the minimum norm solution of an underdetermined system A * X = B.
  3. If TRANS = 'C' and m >= n: find the minimum norm solution of an undetermined system A**T * X = B.
  4. If TRANS = 'C' and m < n: find the least squares solution of an overdetermined system, i.e., solve the least squares problem minimize || B - A**T * X ||.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X.

Parameters
[in]TRANS
          TRANS is CHARACTER*1
          = 'N': the linear system involves A;
          = 'C': the linear system involves A**C.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of
          columns of the matrices B and X. NRHS >=0.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          A is overwritten by details of its QR or LQ
          factorization as returned by CGEQR or CGELQ.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in,out]B
          B is COMPLEX array, dimension (LDB,NRHS)
          On entry, the matrix B of right hand side vectors, stored
          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
          if TRANS = 'C'.
          On exit, if INFO = 0, B is overwritten by the solution
          vectors, stored columnwise:
          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
          squares solution vectors.
          if TRANS = 'N' and m < n, rows 1 to N of B contain the
          minimum norm solution vectors;
          if TRANS = 'C' and m >= n, rows 1 to M of B contain the
          minimum norm solution vectors;
          if TRANS = 'C' and m < n, rows 1 to M of B contain the
          least squares solution vectors.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= MAX(1,M,N).
[out]WORK
          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
          or optimal, if query was assumed) LWORK.
          See LWORK for details.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          If LWORK = -1 or -2, then a workspace query is assumed.
          If LWORK = -1, the routine calculates optimal size of WORK for the
          optimal performance and returns this value in WORK(1).
          If LWORK = -2, the routine calculates minimal size of WORK and 
          returns this value in WORK(1).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO =  i, the i-th diagonal element of the
                triangular factor of A is zero, so that A does not have
                full rank; the least squares solution could not be
                computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 162 of file cgetsls.f.

162 *
163 * -- LAPACK driver routine (version 3.7.0) --
164 * -- LAPACK is a software package provided by Univ. of Tennessee, --
165 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166 * December 2016
167 *
168 * .. Scalar Arguments ..
169  CHARACTER trans
170  INTEGER info, lda, ldb, lwork, m, n, nrhs
171 * ..
172 * .. Array Arguments ..
173  COMPLEX a( lda, * ), b( ldb, * ), work( * )
174 *
175 * ..
176 *
177 * =====================================================================
178 *
179 * .. Parameters ..
180  REAL zero, one
181  parameter ( zero = 0.0e0, one = 1.0e0 )
182  COMPLEX czero
183  parameter ( czero = ( 0.0e+0, 0.0e+0 ) )
184 * ..
185 * .. Local Scalars ..
186  LOGICAL lquery, tran
187  INTEGER i, iascl, ibscl, j, minmn, maxmn, brow,
188  $ scllen, mnk, tszo, tszm, lwo, lwm, lw1, lw2,
189  $ wsizeo, wsizem, info2
190  REAL anrm, bignum, bnrm, smlnum
191  COMPLEX tq( 5 ), workq
192 * ..
193 * .. External Functions ..
194  LOGICAL lsame
195  INTEGER ilaenv
196  REAL slamch, clange
197  EXTERNAL lsame, ilaenv, slabad, slamch, clange
198 * ..
199 * .. External Subroutines ..
200  EXTERNAL cgeqr, cgemqr, clascl, claset,
201  $ ctrtrs, xerbla, cgelq, cgemlq
202 * ..
203 * .. Intrinsic Functions ..
204  INTRINSIC REAL, max, min, int
205 * ..
206 * .. Executable Statements ..
207 *
208 * Test the input arguments.
209 *
210  info = 0
211  minmn = min( m, n )
212  maxmn = max( m, n )
213  mnk = max( minmn, nrhs )
214  tran = lsame( trans, 'C' )
215 *
216  lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
217  IF( .NOT.( lsame( trans, 'N' ) .OR.
218  $ lsame( trans, 'C' ) ) ) THEN
219  info = -1
220  ELSE IF( m.LT.0 ) THEN
221  info = -2
222  ELSE IF( n.LT.0 ) THEN
223  info = -3
224  ELSE IF( nrhs.LT.0 ) THEN
225  info = -4
226  ELSE IF( lda.LT.max( 1, m ) ) THEN
227  info = -6
228  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
229  info = -8
230  END IF
231 *
232  IF( info.EQ.0 ) THEN
233 *
234 * Determine the block size and minimum LWORK
235 *
236  IF( m.GE.n ) THEN
237  CALL cgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
238  tszo = int( tq( 1 ) )
239  lwo = int( workq )
240  CALL cgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
241  $ tszo, b, ldb, workq, -1, info2 )
242  lwo = max( lwo, int( workq ) )
243  CALL cgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
244  tszm = int( tq( 1 ) )
245  lwm = int( workq )
246  CALL cgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
247  $ tszm, b, ldb, workq, -1, info2 )
248  lwm = max( lwm, int( workq ) )
249  wsizeo = tszo + lwo
250  wsizem = tszm + lwm
251  ELSE
252  CALL cgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
253  tszo = int( tq( 1 ) )
254  lwo = int( workq )
255  CALL cgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
256  $ tszo, b, ldb, workq, -1, info2 )
257  lwo = max( lwo, int( workq ) )
258  CALL cgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
259  tszm = int( tq( 1 ) )
260  lwm = int( workq )
261  CALL cgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
262  $ tszo, b, ldb, workq, -1, info2 )
263  lwm = max( lwm, int( workq ) )
264  wsizeo = tszo + lwo
265  wsizem = tszm + lwm
266  END IF
267 *
268  IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
269  info = -10
270  END IF
271 *
272  END IF
273 *
274  IF( info.NE.0 ) THEN
275  CALL xerbla( 'CGETSLS', -info )
276  work( 1 ) = REAL( wsizeo )
277  RETURN
278  END IF
279  IF( lquery ) THEN
280  IF( lwork.EQ.-1 ) work( 1 ) = REAL( wsizeo )
281  IF( lwork.EQ.-2 ) work( 1 ) = REAL( wsizem )
282  RETURN
283  END IF
284  IF( lwork.LT.wsizeo ) THEN
285  lw1 = tszm
286  lw2 = lwm
287  ELSE
288  lw1 = tszo
289  lw2 = lwo
290  END IF
291 *
292 * Quick return if possible
293 *
294  IF( min( m, n, nrhs ).EQ.0 ) THEN
295  CALL claset( 'FULL', max( m, n ), nrhs, czero, czero,
296  $ b, ldb )
297  RETURN
298  END IF
299 *
300 * Get machine parameters
301 *
302  smlnum = slamch( 'S' ) / slamch( 'P' )
303  bignum = one / smlnum
304  CALL slabad( smlnum, bignum )
305 *
306 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
307 *
308  anrm = clange( 'M', m, n, a, lda, work )
309  iascl = 0
310  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
311 *
312 * Scale matrix norm up to SMLNUM
313 *
314  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
315  iascl = 1
316  ELSE IF( anrm.GT.bignum ) THEN
317 *
318 * Scale matrix norm down to BIGNUM
319 *
320  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
321  iascl = 2
322  ELSE IF( anrm.EQ.zero ) THEN
323 *
324 * Matrix all zero. Return zero solution.
325 *
326  CALL claset( 'F', maxmn, nrhs, czero, czero, b, ldb )
327  GO TO 50
328  END IF
329 *
330  brow = m
331  IF ( tran ) THEN
332  brow = n
333  END IF
334  bnrm = clange( 'M', brow, nrhs, b, ldb, work )
335  ibscl = 0
336  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
337 *
338 * Scale matrix norm up to SMLNUM
339 *
340  CALL clascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
341  $ info )
342  ibscl = 1
343  ELSE IF( bnrm.GT.bignum ) THEN
344 *
345 * Scale matrix norm down to BIGNUM
346 *
347  CALL clascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
348  $ info )
349  ibscl = 2
350  END IF
351 *
352  IF ( m.GE.n ) THEN
353 *
354 * compute QR factorization of A
355 *
356  CALL cgeqr( m, n, a, lda, work( lw2+1 ), lw1,
357  $ work( 1 ), lw2, info )
358  IF ( .NOT.tran ) THEN
359 *
360 * Least-Squares Problem min || A * X - B ||
361 *
362 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
363 *
364  CALL cgemqr( 'L' , 'C', m, nrhs, n, a, lda,
365  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
366  $ info )
367 *
368 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
369 *
370  CALL ctrtrs( 'U', 'N', 'N', n, nrhs,
371  $ a, lda, b, ldb, info )
372  IF( info.GT.0 ) THEN
373  RETURN
374  END IF
375  scllen = n
376  ELSE
377 *
378 * Overdetermined system of equations A**T * X = B
379 *
380 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
381 *
382  CALL ctrtrs( 'U', 'C', 'N', n, nrhs,
383  $ a, lda, b, ldb, info )
384 *
385  IF( info.GT.0 ) THEN
386  RETURN
387  END IF
388 *
389 * B(N+1:M,1:NRHS) = CZERO
390 *
391  DO 20 j = 1, nrhs
392  DO 10 i = n + 1, m
393  b( i, j ) = czero
394  10 CONTINUE
395  20 CONTINUE
396 *
397 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
398 *
399  CALL cgemqr( 'L', 'N', m, nrhs, n, a, lda,
400  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
401  $ info )
402 *
403  scllen = m
404 *
405  END IF
406 *
407  ELSE
408 *
409 * Compute LQ factorization of A
410 *
411  CALL cgelq( m, n, a, lda, work( lw2+1 ), lw1,
412  $ work( 1 ), lw2, info )
413 *
414 * workspace at least M, optimally M*NB.
415 *
416  IF( .NOT.tran ) THEN
417 *
418 * underdetermined system of equations A * X = B
419 *
420 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
421 *
422  CALL ctrtrs( 'L', 'N', 'N', m, nrhs,
423  $ a, lda, b, ldb, info )
424 *
425  IF( info.GT.0 ) THEN
426  RETURN
427  END IF
428 *
429 * B(M+1:N,1:NRHS) = 0
430 *
431  DO 40 j = 1, nrhs
432  DO 30 i = m + 1, n
433  b( i, j ) = czero
434  30 CONTINUE
435  40 CONTINUE
436 *
437 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
438 *
439  CALL cgemlq( 'L', 'C', n, nrhs, m, a, lda,
440  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
441  $ info )
442 *
443 * workspace at least NRHS, optimally NRHS*NB
444 *
445  scllen = n
446 *
447  ELSE
448 *
449 * overdetermined system min || A**T * X - B ||
450 *
451 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
452 *
453  CALL cgemlq( 'L', 'N', n, nrhs, m, a, lda,
454  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
455  $ info )
456 *
457 * workspace at least NRHS, optimally NRHS*NB
458 *
459 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
460 *
461  CALL ctrtrs( 'L', 'C', 'N', m, nrhs,
462  $ a, lda, b, ldb, info )
463 *
464  IF( info.GT.0 ) THEN
465  RETURN
466  END IF
467 *
468  scllen = m
469 *
470  END IF
471 *
472  END IF
473 *
474 * Undo scaling
475 *
476  IF( iascl.EQ.1 ) THEN
477  CALL clascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
478  $ info )
479  ELSE IF( iascl.EQ.2 ) THEN
480  CALL clascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
481  $ info )
482  END IF
483  IF( ibscl.EQ.1 ) THEN
484  CALL clascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
485  $ info )
486  ELSE IF( ibscl.EQ.2 ) THEN
487  CALL clascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
488  $ info )
489  END IF
490 *
491  50 CONTINUE
492  work( 1 ) = REAL( tszo + lwo )
493  RETURN
494 *
495 * End of ZGETSLS
496 *
cgemqr
subroutine cgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
Definition: cgemqr.f:171
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
cgelq
subroutine cgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
Definition: cgelq.f:161
ctrtrs
subroutine ctrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
CTRTRS
Definition: ctrtrs.f:142
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
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
clange
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
cgemlq
subroutine cgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
Definition: cgemlq.f:169
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cgeqr
subroutine cgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
Definition: cgeqr.f:162
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76

Here is the call graph for this function:

Here is the caller graph for this function: