LAPACK 3.11.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
228 EXTERNAL lsame, ilaenv, slamch, clange
229* ..
230* .. External Subroutines ..
231 EXTERNAL cgelqt, cgeqrt, cgemlqt, cgemqrt, slabad,
233* ..
234* .. Intrinsic Functions ..
235 INTRINSIC real, 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 ) = real( 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 ) = real( 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 CALL slabad( smlnum, bignum )
315*
316* Scale A, B if max element outside range [SMLNUM,BIGNUM]
317*
318 anrm = clange( 'M', m, n, a, lda, rwork )
319 iascl = 0
320 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
321*
322* Scale matrix norm up to SMLNUM
323*
324 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
325 iascl = 1
326 ELSE IF( anrm.GT.bignum ) THEN
327*
328* Scale matrix norm down to BIGNUM
329*
330 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
331 iascl = 2
332 ELSE IF( anrm.EQ.zero ) THEN
333*
334* Matrix all zero. Return zero solution.
335*
336 CALL claset( 'Full', max( m, n ), nrhs, czero, czero, b, ldb )
337 work( 1 ) = real( lwopt )
338 RETURN
339 END IF
340*
341 brow = m
342 IF( tpsd )
343 $ brow = n
344 bnrm = clange( 'M', brow, nrhs, b, ldb, rwork )
345 ibscl = 0
346 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
347*
348* Scale matrix norm up to SMLNUM
349*
350 CALL clascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
351 $ info )
352 ibscl = 1
353 ELSE IF( bnrm.GT.bignum ) THEN
354*
355* Scale matrix norm down to BIGNUM
356*
357 CALL clascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
358 $ info )
359 ibscl = 2
360 END IF
361*
362 IF( m.GE.n ) THEN
363*
364* M > N:
365* Compute the blocked QR factorization of A,
366* using the compact WY representation of Q,
367* workspace at least N, optimally N*NB.
368*
369 CALL cgeqrt( m, n, nb, a, lda, work( 1 ), nb,
370 $ work( mn*nb+1 ), info )
371*
372 IF( .NOT.tpsd ) THEN
373*
374* M > N, A is not transposed:
375* Overdetermined system of equations,
376* least-squares problem, min || A * X - B ||.
377*
378* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
379* using the compact WY representation of Q,
380* workspace at least NRHS, optimally NRHS*NB.
381*
382 CALL cgemqrt( 'Left', 'Conjugate transpose', m, nrhs, n, nb,
383 $ a, lda, work( 1 ), nb, b, ldb,
384 $ work( mn*nb+1 ), info )
385*
386* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
387*
388 CALL ctrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
389 $ a, lda, b, ldb, info )
390*
391 IF( info.GT.0 ) THEN
392 RETURN
393 END IF
394*
395 scllen = n
396*
397 ELSE
398*
399* M > N, A is transposed:
400* Underdetermined system of equations,
401* minimum norm solution of A**T * X = B.
402*
403* Compute B := inv(R**T) * B in two row blocks of B.
404*
405* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
406*
407 CALL ctrtrs( 'Upper', 'Conjugate transpose', 'Non-unit',
408 $ n, nrhs, a, lda, b, ldb, info )
409*
410 IF( info.GT.0 ) THEN
411 RETURN
412 END IF
413*
414* Block 2: Zero out all rows below the N-th row in B:
415* B(N+1:M,1:NRHS) = ZERO
416*
417 DO j = 1, nrhs
418 DO i = n + 1, m
419 b( i, j ) = zero
420 END DO
421 END DO
422*
423* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
424* using the compact WY representation of Q,
425* workspace at least NRHS, optimally NRHS*NB.
426*
427 CALL cgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
428 $ a, lda, work( 1 ), nb, b, ldb,
429 $ work( mn*nb+1 ), info )
430*
431 scllen = m
432*
433 END IF
434*
435 ELSE
436*
437* M < N:
438* Compute the blocked LQ factorization of A,
439* using the compact WY representation of Q,
440* workspace at least M, optimally M*NB.
441*
442 CALL cgelqt( m, n, nb, a, lda, work( 1 ), nb,
443 $ work( mn*nb+1 ), info )
444*
445 IF( .NOT.tpsd ) THEN
446*
447* M < N, A is not transposed:
448* Underdetermined system of equations,
449* minimum norm solution of A * X = B.
450*
451* Compute B := inv(L) * B in two row blocks of B.
452*
453* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
454*
455 CALL ctrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
456 $ a, lda, b, ldb, info )
457*
458 IF( info.GT.0 ) THEN
459 RETURN
460 END IF
461*
462* Block 2: Zero out all rows below the M-th row in B:
463* B(M+1:N,1:NRHS) = ZERO
464*
465 DO j = 1, nrhs
466 DO i = m + 1, n
467 b( i, j ) = zero
468 END DO
469 END DO
470*
471* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
472* using the compact WY representation of Q,
473* workspace at least NRHS, optimally NRHS*NB.
474*
475 CALL cgemlqt( 'Left', 'Conjugate transpose', n, nrhs, m, nb,
476 $ a, lda, work( 1 ), nb, b, ldb,
477 $ work( mn*nb+1 ), info )
478*
479 scllen = n
480*
481 ELSE
482*
483* M < N, A is transposed:
484* Overdetermined system of equations,
485* least-squares problem, min || A**T * X - B ||.
486*
487* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
488* using the compact WY representation of Q,
489* workspace at least NRHS, optimally NRHS*NB.
490*
491 CALL cgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
492 $ a, lda, work( 1 ), nb, b, ldb,
493 $ work( mn*nb+1), info )
494*
495* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
496*
497 CALL ctrtrs( 'Lower', 'Conjugate transpose', 'Non-unit',
498 $ m, nrhs, a, lda, b, ldb, info )
499*
500 IF( info.GT.0 ) THEN
501 RETURN
502 END IF
503*
504 scllen = m
505*
506 END IF
507*
508 END IF
509*
510* Undo scaling
511*
512 IF( iascl.EQ.1 ) THEN
513 CALL clascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
514 $ info )
515 ELSE IF( iascl.EQ.2 ) THEN
516 CALL clascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
517 $ info )
518 END IF
519 IF( ibscl.EQ.1 ) THEN
520 CALL clascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
521 $ info )
522 ELSE IF( ibscl.EQ.2 ) THEN
523 CALL clascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
524 $ info )
525 END IF
526*
527 work( 1 ) = real( lwopt )
528*
529 RETURN
530*
531* End of CGELST
532*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 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
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
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 ctrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
CTRTRS
Definition: ctrtrs.f:140
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: