LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cgelst()

subroutine cgelst ( 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 
)

CGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.

Download CGELST + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CGELST solves overdetermined or underdetermined real linear systems
 involving an M-by-N matrix A, or its conjugate-transpose, using a QR
 or LQ factorization of A with compact WY representation of Q.
 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 underdetermined 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**H.
[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,
            if M >= N, A is overwritten by details of its QR
                       factorization as returned by CGEQRT;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by CGELQT.
[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; the residual sum of squares for the
          solution in each column is given by the sum of squares of
          modulus of elements N+1 to M in that column;
          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; the residual sum of squares
          for the solution in each column is given by the sum of
          squares of the modulus of elements M+1 to N in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= MAX(1,M,N).
[out]WORK
          WORK is COMPLEX array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= max( 1, MN + max( MN, NRHS ) ).
          For optimal performance,
          LWORK >= max( 1, (MN + max( MN, NRHS ))*NB ).
          where MN = min(M,N) and NB is the optimum block size.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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.
Contributors:
  November 2022,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

Definition at line 192 of file cgelst.f.

194*
195* -- LAPACK driver routine --
196* -- LAPACK is a software package provided by Univ. of Tennessee, --
197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198*
199* .. Scalar Arguments ..
200 CHARACTER TRANS
201 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
202* ..
203* .. Array Arguments ..
204 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 REAL ZERO, ONE
211 parameter( zero = 0.0e+0, one = 1.0e+0 )
212 COMPLEX CZERO
213 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
214* ..
215* .. Local Scalars ..
216 LOGICAL LQUERY, TPSD
217 INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
218 $ NB, NBMIN, SCLLEN
219 REAL ANRM, BIGNUM, BNRM, SMLNUM
220* ..
221* .. Local Arrays ..
222 REAL RWORK( 1 )
223* ..
224* .. External Functions ..
225 LOGICAL LSAME
226 INTEGER ILAENV
227 REAL SLAMCH, CLANGE, SROUNDUP_LWORK
229* ..
230* .. External Subroutines ..
231 EXTERNAL cgelqt, cgeqrt, cgemlqt, cgemqrt,
233* ..
234* .. Intrinsic Functions ..
235 INTRINSIC max, min
236* ..
237* .. Executable Statements ..
238*
239* Test the input arguments.
240*
241 info = 0
242 mn = min( m, n )
243 lquery = ( lwork.EQ.-1 )
244 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'C' ) ) ) THEN
245 info = -1
246 ELSE IF( m.LT.0 ) THEN
247 info = -2
248 ELSE IF( n.LT.0 ) THEN
249 info = -3
250 ELSE IF( nrhs.LT.0 ) THEN
251 info = -4
252 ELSE IF( lda.LT.max( 1, m ) ) THEN
253 info = -6
254 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
255 info = -8
256 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND. .NOT.lquery )
257 $ THEN
258 info = -10
259 END IF
260*
261* Figure out optimal block size and optimal workspace size
262*
263 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
264*
265 tpsd = .true.
266 IF( lsame( trans, 'N' ) )
267 $ tpsd = .false.
268*
269 nb = ilaenv( 1, 'CGELST', ' ', m, n, -1, -1 )
270*
271 mnnrhs = max( mn, nrhs )
272 lwopt = max( 1, (mn+mnnrhs)*nb )
273 work( 1 ) = sroundup_lwork( lwopt )
274*
275 END IF
276*
277 IF( info.NE.0 ) THEN
278 CALL xerbla( 'CGELST ', -info )
279 RETURN
280 ELSE IF( lquery ) THEN
281 RETURN
282 END IF
283*
284* Quick return if possible
285*
286 IF( min( m, n, nrhs ).EQ.0 ) THEN
287 CALL claset( 'Full', max( m, n ), nrhs, czero, czero, b, ldb )
288 work( 1 ) = sroundup_lwork( lwopt )
289 RETURN
290 END IF
291*
292* *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
293*
294 IF( nb.GT.mn ) nb = mn
295*
296* Determine the block size from the supplied LWORK
297* ( at this stage we know that LWORK >= (minimum required workspace,
298* but it may be less than optimal)
299*
300 nb = min( nb, lwork/( mn + mnnrhs ) )
301*
302* The minimum value of NB, when blocked code is used
303*
304 nbmin = max( 2, ilaenv( 2, 'CGELST', ' ', m, n, -1, -1 ) )
305*
306 IF( nb.LT.nbmin ) THEN
307 nb = 1
308 END IF
309*
310* Get machine parameters
311*
312 smlnum = slamch( 'S' ) / slamch( 'P' )
313 bignum = one / smlnum
314*
315* Scale A, B if max element outside range [SMLNUM,BIGNUM]
316*
317 anrm = clange( 'M', m, n, a, lda, rwork )
318 iascl = 0
319 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
320*
321* Scale matrix norm up to SMLNUM
322*
323 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
324 iascl = 1
325 ELSE IF( anrm.GT.bignum ) THEN
326*
327* Scale matrix norm down to BIGNUM
328*
329 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
330 iascl = 2
331 ELSE IF( anrm.EQ.zero ) THEN
332*
333* Matrix all zero. Return zero solution.
334*
335 CALL claset( 'Full', max( m, n ), nrhs, czero, czero, b, ldb )
336 work( 1 ) = sroundup_lwork( lwopt )
337 RETURN
338 END IF
339*
340 brow = m
341 IF( tpsd )
342 $ brow = n
343 bnrm = clange( 'M', brow, nrhs, b, ldb, rwork )
344 ibscl = 0
345 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
346*
347* Scale matrix norm up to SMLNUM
348*
349 CALL clascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
350 $ info )
351 ibscl = 1
352 ELSE IF( bnrm.GT.bignum ) THEN
353*
354* Scale matrix norm down to BIGNUM
355*
356 CALL clascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
357 $ info )
358 ibscl = 2
359 END IF
360*
361 IF( m.GE.n ) THEN
362*
363* M > N:
364* Compute the blocked QR factorization of A,
365* using the compact WY representation of Q,
366* workspace at least N, optimally N*NB.
367*
368 CALL cgeqrt( m, n, nb, a, lda, work( 1 ), nb,
369 $ work( mn*nb+1 ), info )
370*
371 IF( .NOT.tpsd ) THEN
372*
373* M > N, A is not transposed:
374* Overdetermined system of equations,
375* least-squares problem, min || A * X - B ||.
376*
377* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
378* using the compact WY representation of Q,
379* workspace at least NRHS, optimally NRHS*NB.
380*
381 CALL cgemqrt( 'Left', 'Conjugate transpose', m, nrhs, n, nb,
382 $ a, lda, work( 1 ), nb, b, ldb,
383 $ work( mn*nb+1 ), info )
384*
385* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
386*
387 CALL ctrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
388 $ a, lda, b, ldb, info )
389*
390 IF( info.GT.0 ) THEN
391 RETURN
392 END IF
393*
394 scllen = n
395*
396 ELSE
397*
398* M > N, A is transposed:
399* Underdetermined system of equations,
400* minimum norm solution of A**T * X = B.
401*
402* Compute B := inv(R**T) * B in two row blocks of B.
403*
404* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
405*
406 CALL ctrtrs( 'Upper', 'Conjugate transpose', 'Non-unit',
407 $ n, nrhs, a, lda, b, ldb, info )
408*
409 IF( info.GT.0 ) THEN
410 RETURN
411 END IF
412*
413* Block 2: Zero out all rows below the N-th row in B:
414* B(N+1:M,1:NRHS) = ZERO
415*
416 DO j = 1, nrhs
417 DO i = n + 1, m
418 b( i, j ) = zero
419 END DO
420 END DO
421*
422* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
423* using the compact WY representation of Q,
424* workspace at least NRHS, optimally NRHS*NB.
425*
426 CALL cgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
427 $ a, lda, work( 1 ), nb, b, ldb,
428 $ work( mn*nb+1 ), info )
429*
430 scllen = m
431*
432 END IF
433*
434 ELSE
435*
436* M < N:
437* Compute the blocked LQ factorization of A,
438* using the compact WY representation of Q,
439* workspace at least M, optimally M*NB.
440*
441 CALL cgelqt( m, n, nb, a, lda, work( 1 ), nb,
442 $ work( mn*nb+1 ), info )
443*
444 IF( .NOT.tpsd ) THEN
445*
446* M < N, A is not transposed:
447* Underdetermined system of equations,
448* minimum norm solution of A * X = B.
449*
450* Compute B := inv(L) * B in two row blocks of B.
451*
452* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
453*
454 CALL ctrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
455 $ a, lda, b, ldb, info )
456*
457 IF( info.GT.0 ) THEN
458 RETURN
459 END IF
460*
461* Block 2: Zero out all rows below the M-th row in B:
462* B(M+1:N,1:NRHS) = ZERO
463*
464 DO j = 1, nrhs
465 DO i = m + 1, n
466 b( i, j ) = zero
467 END DO
468 END DO
469*
470* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
471* using the compact WY representation of Q,
472* workspace at least NRHS, optimally NRHS*NB.
473*
474 CALL cgemlqt( 'Left', 'Conjugate transpose', n, nrhs, m, nb,
475 $ a, lda, work( 1 ), nb, b, ldb,
476 $ work( mn*nb+1 ), info )
477*
478 scllen = n
479*
480 ELSE
481*
482* M < N, A is transposed:
483* Overdetermined system of equations,
484* least-squares problem, min || A**T * X - B ||.
485*
486* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
487* using the compact WY representation of Q,
488* workspace at least NRHS, optimally NRHS*NB.
489*
490 CALL cgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
491 $ a, lda, work( 1 ), nb, b, ldb,
492 $ work( mn*nb+1), info )
493*
494* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
495*
496 CALL ctrtrs( 'Lower', 'Conjugate transpose', 'Non-unit',
497 $ m, nrhs, a, lda, b, ldb, info )
498*
499 IF( info.GT.0 ) THEN
500 RETURN
501 END IF
502*
503 scllen = m
504*
505 END IF
506*
507 END IF
508*
509* Undo scaling
510*
511 IF( iascl.EQ.1 ) THEN
512 CALL clascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
513 $ info )
514 ELSE IF( iascl.EQ.2 ) THEN
515 CALL clascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
516 $ info )
517 END IF
518 IF( ibscl.EQ.1 ) THEN
519 CALL clascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
520 $ info )
521 ELSE IF( ibscl.EQ.2 ) THEN
522 CALL clascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
523 $ info )
524 END IF
525*
526 work( 1 ) = sroundup_lwork( lwopt )
527*
528 RETURN
529*
530* End of CGELST
531*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgelqt(m, n, mb, a, lda, t, ldt, work, info)
CGELQT
Definition cgelqt.f:124
subroutine cgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
CGEMLQT
Definition cgemlqt.f:153
subroutine cgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
CGEMQRT
Definition cgemqrt.f:168
subroutine cgeqrt(m, n, nb, a, lda, t, ldt, work, info)
CGEQRT
Definition cgeqrt.f:141
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:115
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:143
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:106
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ctrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
CTRTRS
Definition ctrtrs.f:140
Here is the call graph for this function:
Here is the caller graph for this function: