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

◆ cgels()

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

CGELS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
 CGELS solves overdetermined or underdetermined complex linear systems
 involving an M-by-N matrix A, or its conjugate-transpose, using a QR
 or 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 underdetermined system A**H * 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**H * 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.
            if M >= N, A is overwritten by details of its QR
                       factorization as returned by CGEQRF;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by CGELQF.
[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 the
          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.

Definition at line 180 of file cgels.f.

182*
183* -- LAPACK driver routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 CHARACTER TRANS
189 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
190* ..
191* .. Array Arguments ..
192 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
193* ..
194*
195* =====================================================================
196*
197* .. Parameters ..
198 REAL ZERO, ONE
199 parameter( zero = 0.0e+0, one = 1.0e+0 )
200 COMPLEX CZERO
201 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
202* ..
203* .. Local Scalars ..
204 LOGICAL LQUERY, TPSD
205 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
206 REAL ANRM, BIGNUM, BNRM, SMLNUM
207* ..
208* .. Local Arrays ..
209 REAL RWORK( 1 )
210* ..
211* .. External Functions ..
212 LOGICAL LSAME
213 INTEGER ILAENV
214 REAL CLANGE, SLAMCH
215 EXTERNAL lsame, ilaenv, clange, slamch
216* ..
217* .. External Subroutines ..
218 EXTERNAL cgelqf, cgeqrf, clascl, claset, ctrtrs, cunmlq,
220* ..
221* .. Intrinsic Functions ..
222 INTRINSIC max, min, real
223* ..
224* .. Executable Statements ..
225*
226* Test the input arguments.
227*
228 info = 0
229 mn = min( m, n )
230 lquery = ( lwork.EQ.-1 )
231 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'C' ) ) ) THEN
232 info = -1
233 ELSE IF( m.LT.0 ) THEN
234 info = -2
235 ELSE IF( n.LT.0 ) THEN
236 info = -3
237 ELSE IF( nrhs.LT.0 ) THEN
238 info = -4
239 ELSE IF( lda.LT.max( 1, m ) ) THEN
240 info = -6
241 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
242 info = -8
243 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND.
244 $ .NOT.lquery ) THEN
245 info = -10
246 END IF
247*
248* Figure out optimal block size
249*
250 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
251*
252 tpsd = .true.
253 IF( lsame( trans, 'N' ) )
254 $ tpsd = .false.
255*
256 IF( m.GE.n ) THEN
257 nb = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
258 IF( tpsd ) THEN
259 nb = max( nb, ilaenv( 1, 'CUNMQR', 'LN', m, nrhs, n,
260 $ -1 ) )
261 ELSE
262 nb = max( nb, ilaenv( 1, 'CUNMQR', 'LC', m, nrhs, n,
263 $ -1 ) )
264 END IF
265 ELSE
266 nb = ilaenv( 1, 'CGELQF', ' ', m, n, -1, -1 )
267 IF( tpsd ) THEN
268 nb = max( nb, ilaenv( 1, 'CUNMLQ', 'LC', n, nrhs, m,
269 $ -1 ) )
270 ELSE
271 nb = max( nb, ilaenv( 1, 'CUNMLQ', 'LN', n, nrhs, m,
272 $ -1 ) )
273 END IF
274 END IF
275*
276 wsize = max( 1, mn + max( mn, nrhs )*nb )
277 work( 1 ) = real( wsize )
278*
279 END IF
280*
281 IF( info.NE.0 ) THEN
282 CALL xerbla( 'CGELS ', -info )
283 RETURN
284 ELSE IF( lquery ) THEN
285 RETURN
286 END IF
287*
288* Quick return if possible
289*
290 IF( min( m, n, nrhs ).EQ.0 ) THEN
291 CALL claset( 'Full', max( m, n ), nrhs, czero, czero, b, ldb )
292 RETURN
293 END IF
294*
295* Get machine parameters
296*
297 smlnum = slamch( 'S' ) / slamch( 'P' )
298 bignum = one / smlnum
299 CALL slabad( smlnum, bignum )
300*
301* Scale A, B if max element outside range [SMLNUM,BIGNUM]
302*
303 anrm = clange( 'M', m, n, a, lda, rwork )
304 iascl = 0
305 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
306*
307* Scale matrix norm up to SMLNUM
308*
309 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
310 iascl = 1
311 ELSE IF( anrm.GT.bignum ) THEN
312*
313* Scale matrix norm down to BIGNUM
314*
315 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
316 iascl = 2
317 ELSE IF( anrm.EQ.zero ) THEN
318*
319* Matrix all zero. Return zero solution.
320*
321 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
322 GO TO 50
323 END IF
324*
325 brow = m
326 IF( tpsd )
327 $ brow = n
328 bnrm = clange( 'M', brow, nrhs, b, ldb, rwork )
329 ibscl = 0
330 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
331*
332* Scale matrix norm up to SMLNUM
333*
334 CALL clascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
335 $ info )
336 ibscl = 1
337 ELSE IF( bnrm.GT.bignum ) THEN
338*
339* Scale matrix norm down to BIGNUM
340*
341 CALL clascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
342 $ info )
343 ibscl = 2
344 END IF
345*
346 IF( m.GE.n ) THEN
347*
348* compute QR factorization of A
349*
350 CALL cgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
351 $ info )
352*
353* workspace at least N, optimally N*NB
354*
355 IF( .NOT.tpsd ) THEN
356*
357* Least-Squares Problem min || A * X - B ||
358*
359* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
360*
361 CALL cunmqr( 'Left', 'Conjugate transpose', m, nrhs, n, a,
362 $ lda, work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
363 $ info )
364*
365* workspace at least NRHS, optimally NRHS*NB
366*
367* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
368*
369 CALL ctrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
370 $ a, lda, b, ldb, info )
371*
372 IF( info.GT.0 ) THEN
373 RETURN
374 END IF
375*
376 scllen = n
377*
378 ELSE
379*
380* Underdetermined system of equations A**T * X = B
381*
382* B(1:N,1:NRHS) := inv(R**H) * B(1:N,1:NRHS)
383*
384 CALL ctrtrs( 'Upper', 'Conjugate transpose','Non-unit',
385 $ n, nrhs, 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) = ZERO
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 cunmqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
402 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
403 $ info )
404*
405* workspace at least NRHS, optimally NRHS*NB
406*
407 scllen = m
408*
409 END IF
410*
411 ELSE
412*
413* Compute LQ factorization of A
414*
415 CALL cgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
416 $ info )
417*
418* workspace at least M, optimally M*NB.
419*
420 IF( .NOT.tpsd ) THEN
421*
422* underdetermined system of equations A * X = B
423*
424* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
425*
426 CALL ctrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
427 $ a, lda, b, ldb, info )
428*
429 IF( info.GT.0 ) THEN
430 RETURN
431 END IF
432*
433* B(M+1:N,1:NRHS) = 0
434*
435 DO 40 j = 1, nrhs
436 DO 30 i = m + 1, n
437 b( i, j ) = czero
438 30 CONTINUE
439 40 CONTINUE
440*
441* B(1:N,1:NRHS) := Q(1:N,:)**H * B(1:M,1:NRHS)
442*
443 CALL cunmlq( 'Left', 'Conjugate transpose', n, nrhs, m, a,
444 $ lda, work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
445 $ info )
446*
447* workspace at least NRHS, optimally NRHS*NB
448*
449 scllen = n
450*
451 ELSE
452*
453* overdetermined system min || A**H * X - B ||
454*
455* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
456*
457 CALL cunmlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
458 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
459 $ info )
460*
461* workspace at least NRHS, optimally NRHS*NB
462*
463* B(1:M,1:NRHS) := inv(L**H) * B(1:M,1:NRHS)
464*
465 CALL ctrtrs( 'Lower', 'Conjugate transpose', 'Non-unit',
466 $ m, nrhs, a, lda, b, ldb, info )
467*
468 IF( info.GT.0 ) THEN
469 RETURN
470 END IF
471*
472 scllen = m
473*
474 END IF
475*
476 END IF
477*
478* Undo scaling
479*
480 IF( iascl.EQ.1 ) THEN
481 CALL clascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
482 $ info )
483 ELSE IF( iascl.EQ.2 ) THEN
484 CALL clascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
485 $ info )
486 END IF
487 IF( ibscl.EQ.1 ) THEN
488 CALL clascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
489 $ info )
490 ELSE IF( ibscl.EQ.2 ) THEN
491 CALL clascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
492 $ info )
493 END IF
494*
495 50 CONTINUE
496 work( 1 ) = real( wsize )
497*
498 RETURN
499*
500* End of CGELS
501*
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 cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:146
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.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
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 cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
Definition: cunmlq.f:168
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:168
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: