OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
complex
LAPACK » General Matrices » Linear Solve

Functions

subroutine cgels (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
  CGELS solves overdetermined or underdetermined systems for GE matrices
subroutine cgelsd (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info)
  CGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
subroutine cgelss (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, info)
  CGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine cgelsy (m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, rwork, info)
  CGELSY solves overdetermined or underdetermined systems for GE matrices
subroutine cgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
  CGESV computes the solution to system of linear equations A * X = B for GE matrices (simple driver)
subroutine cgesvx (fact, trans, n, nrhs, a, lda, af, ldaf, ipiv, equed, r, c, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
  CGESVX computes the solution to system of linear equations A * X = B for GE matrices
subroutine cgesvxx (fact, trans, n, nrhs, a, lda, af, ldaf, ipiv, equed, r, c, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, rwork, info)
  CGESVXX computes the solution to system of linear equations A * X = B for GE matrices
subroutine cgetsls (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
 CGETSLS
subroutine cgelsx (m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, rwork, info)
  CGELSX solves overdetermined or underdetermined systems for GE matrices

Detailed Description

This is the group of complex solve driver functions for GE matrices

Function Documentation

◆ 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,
219 $ cunmqr, slabad, xerbla
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*
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
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:115
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
cgelqf
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
Definition cgelqf.f:143
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:143
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:106
cunmqr
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
cunmlq
subroutine cunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMLQ
Definition cunmlq.f:168
ctrtrs
subroutine ctrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
CTRTRS
Definition ctrtrs.f:140
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ cgelsd()

subroutine cgelsd ( integer m,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) s,
real rcond,
integer rank,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer info )

CGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices

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

Purpose:
!>
!> CGELSD computes the minimum-norm solution to a real linear least
!> squares problem:
!>     minimize 2-norm(| b - A*x |)
!> using the singular value decomposition (SVD) of A. A is an M-by-N
!> matrix which may be rank-deficient.
!>
!> 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.
!>
!> The problem is solved in three steps:
!> (1) Reduce the coefficient matrix A to bidiagonal form with
!>     Householder transformations, reducing the original problem
!>     into a  (BLS)
!> (2) Solve the BLS using a divide and conquer approach.
!> (3) Apply back all the Householder transformations to solve
!>     the original least squares problem.
!>
!> The effective rank of A is determined by treating as zero those
!> singular values which are less than RCOND times the largest singular
!> value.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[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 has been destroyed.
!>
[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 M-by-NRHS right hand side matrix B.
!>          On exit, B is overwritten by the N-by-NRHS solution matrix X.
!>          If m >= n and RANK = n, the residual sum-of-squares for
!>          the solution in the i-th column is given by the sum of
!>          squares of the modulus of elements n+1:m in that column.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,M,N).
!>
[out]S
!>          S is REAL array, dimension (min(M,N))
!>          The singular values of A in decreasing order.
!>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
!>
[in]RCOND
!>          RCOND is REAL
!>          RCOND is used to determine the effective rank of A.
!>          Singular values S(i) <= RCOND*S(1) are treated as zero.
!>          If RCOND < 0, machine precision is used instead.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the number of singular values
!>          which are greater than RCOND*S(1).
!>
[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 must be at least 1.
!>          The exact minimum amount of workspace needed depends on M,
!>          N and NRHS. As long as LWORK is at least
!>              2 * N + N * NRHS
!>          if M is greater than or equal to N or
!>              2 * M + M * NRHS
!>          if M is less than N, the code will execute correctly.
!>          For good performance, LWORK should generally be larger.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the array WORK and the
!>          minimum sizes of the arrays RWORK and IWORK, and returns
!>          these values as the first entries of the WORK, RWORK and
!>          IWORK arrays, and no error message related to LWORK is issued
!>          by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (MAX(1,LRWORK))
!>          LRWORK >=
!>             10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
!>             MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
!>          if M is greater than or equal to N or
!>             10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
!>             MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
!>          if M is less than N, the code will execute correctly.
!>          SMLSIZ is returned by ILAENV and is equal to the maximum
!>          size of the subproblems at the bottom of the computation
!>          tree (usually about 25), and
!>             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
!>          On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN),
!>          where MINMN = MIN( M,N ).
!>          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  the algorithm for computing the SVD failed to converge;
!>                if INFO = i, i off-diagonal elements of an intermediate
!>                bidiagonal form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 223 of file cgelsd.f.

225*
226* -- LAPACK driver routine --
227* -- LAPACK is a software package provided by Univ. of Tennessee, --
228* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
229*
230* .. Scalar Arguments ..
231 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
232 REAL RCOND
233* ..
234* .. Array Arguments ..
235 INTEGER IWORK( * )
236 REAL RWORK( * ), S( * )
237 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
238* ..
239*
240* =====================================================================
241*
242* .. Parameters ..
243 REAL ZERO, ONE, TWO
244 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
245 COMPLEX CZERO
246 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
247* ..
248* .. Local Scalars ..
249 LOGICAL LQUERY
250 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
251 $ LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
252 $ MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
253 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
254* ..
255* .. External Subroutines ..
256 EXTERNAL cgebrd, cgelqf, cgeqrf, clacpy,
257 $ clalsd, clascl, claset, cunmbr,
258 $ cunmlq, cunmqr, slabad, slascl,
259 $ slaset, xerbla
260* ..
261* .. External Functions ..
262 INTEGER ILAENV
263 REAL CLANGE, SLAMCH
264 EXTERNAL clange, slamch, ilaenv
265* ..
266* .. Intrinsic Functions ..
267 INTRINSIC int, log, max, min, real
268* ..
269* .. Executable Statements ..
270*
271* Test the input arguments.
272*
273 info = 0
274 minmn = min( m, n )
275 maxmn = max( m, n )
276 lquery = ( lwork.EQ.-1 )
277 IF( m.LT.0 ) THEN
278 info = -1
279 ELSE IF( n.LT.0 ) THEN
280 info = -2
281 ELSE IF( nrhs.LT.0 ) THEN
282 info = -3
283 ELSE IF( lda.LT.max( 1, m ) ) THEN
284 info = -5
285 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
286 info = -7
287 END IF
288*
289* Compute workspace.
290* (Note: Comments in the code beginning "Workspace:" describe the
291* minimal amount of workspace needed at that point in the code,
292* as well as the preferred amount for good performance.
293* NB refers to the optimal block size for the immediately
294* following subroutine, as returned by ILAENV.)
295*
296 IF( info.EQ.0 ) THEN
297 minwrk = 1
298 maxwrk = 1
299 liwork = 1
300 lrwork = 1
301 IF( minmn.GT.0 ) THEN
302 smlsiz = ilaenv( 9, 'CGELSD', ' ', 0, 0, 0, 0 )
303 mnthr = ilaenv( 6, 'CGELSD', ' ', m, n, nrhs, -1 )
304 nlvl = max( int( log( real( minmn ) / real( smlsiz + 1 ) ) /
305 $ log( two ) ) + 1, 0 )
306 liwork = 3*minmn*nlvl + 11*minmn
307 mm = m
308 IF( m.GE.n .AND. m.GE.mnthr ) THEN
309*
310* Path 1a - overdetermined, with many more rows than
311* columns.
312*
313 mm = n
314 maxwrk = max( maxwrk, n*ilaenv( 1, 'CGEQRF', ' ', m, n,
315 $ -1, -1 ) )
316 maxwrk = max( maxwrk, nrhs*ilaenv( 1, 'CUNMQR', 'LC', m,
317 $ nrhs, n, -1 ) )
318 END IF
319 IF( m.GE.n ) THEN
320*
321* Path 1 - overdetermined or exactly determined.
322*
323 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs +
324 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
325 maxwrk = max( maxwrk, 2*n + ( mm + n )*ilaenv( 1,
326 $ 'CGEBRD', ' ', mm, n, -1, -1 ) )
327 maxwrk = max( maxwrk, 2*n + nrhs*ilaenv( 1, 'CUNMBR',
328 $ 'QLC', mm, nrhs, n, -1 ) )
329 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
330 $ 'CUNMBR', 'PLN', n, nrhs, n, -1 ) )
331 maxwrk = max( maxwrk, 2*n + n*nrhs )
332 minwrk = max( 2*n + mm, 2*n + n*nrhs )
333 END IF
334 IF( n.GT.m ) THEN
335 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs +
336 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
337 IF( n.GE.mnthr ) THEN
338*
339* Path 2a - underdetermined, with many more columns
340* than rows.
341*
342 maxwrk = m + m*ilaenv( 1, 'CGELQF', ' ', m, n, -1,
343 $ -1 )
344 maxwrk = max( maxwrk, m*m + 4*m + 2*m*ilaenv( 1,
345 $ 'CGEBRD', ' ', m, m, -1, -1 ) )
346 maxwrk = max( maxwrk, m*m + 4*m + nrhs*ilaenv( 1,
347 $ 'CUNMBR', 'QLC', m, nrhs, m, -1 ) )
348 maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*ilaenv( 1,
349 $ 'CUNMLQ', 'LC', n, nrhs, m, -1 ) )
350 IF( nrhs.GT.1 ) THEN
351 maxwrk = max( maxwrk, m*m + m + m*nrhs )
352 ELSE
353 maxwrk = max( maxwrk, m*m + 2*m )
354 END IF
355 maxwrk = max( maxwrk, m*m + 4*m + m*nrhs )
356! XXX: Ensure the Path 2a case below is triggered. The workspace
357! calculation should use queries for all routines eventually.
358 maxwrk = max( maxwrk,
359 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
360 ELSE
361*
362* Path 2 - underdetermined.
363*
364 maxwrk = 2*m + ( n + m )*ilaenv( 1, 'CGEBRD', ' ', m,
365 $ n, -1, -1 )
366 maxwrk = max( maxwrk, 2*m + nrhs*ilaenv( 1, 'CUNMBR',
367 $ 'QLC', m, nrhs, m, -1 ) )
368 maxwrk = max( maxwrk, 2*m + m*ilaenv( 1, 'CUNMBR',
369 $ 'PLN', n, nrhs, m, -1 ) )
370 maxwrk = max( maxwrk, 2*m + m*nrhs )
371 END IF
372 minwrk = max( 2*m + n, 2*m + m*nrhs )
373 END IF
374 END IF
375 minwrk = min( minwrk, maxwrk )
376 work( 1 ) = maxwrk
377 iwork( 1 ) = liwork
378 rwork( 1 ) = lrwork
379*
380 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
381 info = -12
382 END IF
383 END IF
384*
385 IF( info.NE.0 ) THEN
386 CALL xerbla( 'CGELSD', -info )
387 RETURN
388 ELSE IF( lquery ) THEN
389 RETURN
390 END IF
391*
392* Quick return if possible.
393*
394 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
395 rank = 0
396 RETURN
397 END IF
398*
399* Get machine parameters.
400*
401 eps = slamch( 'P' )
402 sfmin = slamch( 'S' )
403 smlnum = sfmin / eps
404 bignum = one / smlnum
405 CALL slabad( smlnum, bignum )
406*
407* Scale A if max entry outside range [SMLNUM,BIGNUM].
408*
409 anrm = clange( 'M', m, n, a, lda, rwork )
410 iascl = 0
411 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
412*
413* Scale matrix norm up to SMLNUM
414*
415 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
416 iascl = 1
417 ELSE IF( anrm.GT.bignum ) THEN
418*
419* Scale matrix norm down to BIGNUM.
420*
421 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
422 iascl = 2
423 ELSE IF( anrm.EQ.zero ) THEN
424*
425* Matrix all zero. Return zero solution.
426*
427 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
428 CALL slaset( 'F', minmn, 1, zero, zero, s, 1 )
429 rank = 0
430 GO TO 10
431 END IF
432*
433* Scale B if max entry outside range [SMLNUM,BIGNUM].
434*
435 bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
436 ibscl = 0
437 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
438*
439* Scale matrix norm up to SMLNUM.
440*
441 CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
442 ibscl = 1
443 ELSE IF( bnrm.GT.bignum ) THEN
444*
445* Scale matrix norm down to BIGNUM.
446*
447 CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
448 ibscl = 2
449 END IF
450*
451* If M < N make sure B(M+1:N,:) = 0
452*
453 IF( m.LT.n )
454 $ CALL claset( 'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
455*
456* Overdetermined case.
457*
458 IF( m.GE.n ) THEN
459*
460* Path 1 - overdetermined or exactly determined.
461*
462 mm = m
463 IF( m.GE.mnthr ) THEN
464*
465* Path 1a - overdetermined, with many more rows than columns
466*
467 mm = n
468 itau = 1
469 nwork = itau + n
470*
471* Compute A=Q*R.
472* (RWorkspace: need N)
473* (CWorkspace: need N, prefer N*NB)
474*
475 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
476 $ lwork-nwork+1, info )
477*
478* Multiply B by transpose(Q).
479* (RWorkspace: need N)
480* (CWorkspace: need NRHS, prefer NRHS*NB)
481*
482 CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, work( itau ), b,
483 $ ldb, work( nwork ), lwork-nwork+1, info )
484*
485* Zero out below R.
486*
487 IF( n.GT.1 ) THEN
488 CALL claset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
489 $ lda )
490 END IF
491 END IF
492*
493 itauq = 1
494 itaup = itauq + n
495 nwork = itaup + n
496 ie = 1
497 nrwork = ie + n
498*
499* Bidiagonalize R in A.
500* (RWorkspace: need N)
501* (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
502*
503 CALL cgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
504 $ work( itaup ), work( nwork ), lwork-nwork+1,
505 $ info )
506*
507* Multiply B by transpose of left bidiagonalizing vectors of R.
508* (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
509*
510 CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, work( itauq ),
511 $ b, ldb, work( nwork ), lwork-nwork+1, info )
512*
513* Solve the bidiagonal least squares problem.
514*
515 CALL clalsd( 'U', smlsiz, n, nrhs, s, rwork( ie ), b, ldb,
516 $ rcond, rank, work( nwork ), rwork( nrwork ),
517 $ iwork, info )
518 IF( info.NE.0 ) THEN
519 GO TO 10
520 END IF
521*
522* Multiply B by right bidiagonalizing vectors of R.
523*
524 CALL cunmbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
525 $ b, ldb, work( nwork ), lwork-nwork+1, info )
526*
527 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
528 $ max( m, 2*m-4, nrhs, n-3*m ) ) THEN
529*
530* Path 2a - underdetermined, with many more columns than rows
531* and sufficient workspace for an efficient algorithm.
532*
533 ldwork = m
534 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
535 $ m*lda+m+m*nrhs ) )ldwork = lda
536 itau = 1
537 nwork = m + 1
538*
539* Compute A=L*Q.
540* (CWorkspace: need 2*M, prefer M+M*NB)
541*
542 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
543 $ lwork-nwork+1, info )
544 il = nwork
545*
546* Copy L to WORK(IL), zeroing out above its diagonal.
547*
548 CALL clacpy( 'L', m, m, a, lda, work( il ), ldwork )
549 CALL claset( 'U', m-1, m-1, czero, czero, work( il+ldwork ),
550 $ ldwork )
551 itauq = il + ldwork*m
552 itaup = itauq + m
553 nwork = itaup + m
554 ie = 1
555 nrwork = ie + m
556*
557* Bidiagonalize L in WORK(IL).
558* (RWorkspace: need M)
559* (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
560*
561 CALL cgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
562 $ work( itauq ), work( itaup ), work( nwork ),
563 $ lwork-nwork+1, info )
564*
565* Multiply B by transpose of left bidiagonalizing vectors of L.
566* (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
567*
568 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, work( il ), ldwork,
569 $ work( itauq ), b, ldb, work( nwork ),
570 $ lwork-nwork+1, info )
571*
572* Solve the bidiagonal least squares problem.
573*
574 CALL clalsd( 'U', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
575 $ rcond, rank, work( nwork ), rwork( nrwork ),
576 $ iwork, info )
577 IF( info.NE.0 ) THEN
578 GO TO 10
579 END IF
580*
581* Multiply B by right bidiagonalizing vectors of L.
582*
583 CALL cunmbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
584 $ work( itaup ), b, ldb, work( nwork ),
585 $ lwork-nwork+1, info )
586*
587* Zero out below first M rows of B.
588*
589 CALL claset( 'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
590 nwork = itau + m
591*
592* Multiply transpose(Q) by B.
593* (CWorkspace: need NRHS, prefer NRHS*NB)
594*
595 CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, work( itau ), b,
596 $ ldb, work( nwork ), lwork-nwork+1, info )
597*
598 ELSE
599*
600* Path 2 - remaining underdetermined cases.
601*
602 itauq = 1
603 itaup = itauq + m
604 nwork = itaup + m
605 ie = 1
606 nrwork = ie + m
607*
608* Bidiagonalize A.
609* (RWorkspace: need M)
610* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
611*
612 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
613 $ work( itaup ), work( nwork ), lwork-nwork+1,
614 $ info )
615*
616* Multiply B by transpose of left bidiagonalizing vectors.
617* (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
618*
619 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda, work( itauq ),
620 $ b, ldb, work( nwork ), lwork-nwork+1, info )
621*
622* Solve the bidiagonal least squares problem.
623*
624 CALL clalsd( 'L', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
625 $ rcond, rank, work( nwork ), rwork( nrwork ),
626 $ iwork, info )
627 IF( info.NE.0 ) THEN
628 GO TO 10
629 END IF
630*
631* Multiply B by right bidiagonalizing vectors of A.
632*
633 CALL cunmbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
634 $ b, ldb, work( nwork ), lwork-nwork+1, info )
635*
636 END IF
637*
638* Undo scaling.
639*
640 IF( iascl.EQ.1 ) THEN
641 CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
642 CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
643 $ info )
644 ELSE IF( iascl.EQ.2 ) THEN
645 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
646 CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
647 $ info )
648 END IF
649 IF( ibscl.EQ.1 ) THEN
650 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
651 ELSE IF( ibscl.EQ.2 ) THEN
652 CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
653 END IF
654*
655 10 CONTINUE
656 work( 1 ) = maxwrk
657 iwork( 1 ) = liwork
658 rwork( 1 ) = lrwork
659 RETURN
660*
661* End of CGELSD
662*
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:143
cgebrd
subroutine cgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
CGEBRD
Definition cgebrd.f:206
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
clalsd
subroutine clalsd(uplo, smlsiz, n, nrhs, d, e, b, ldb, rcond, rank, work, rwork, iwork, info)
CLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition clalsd.f:186
cunmbr
subroutine cunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMBR
Definition cunmbr.f:197

◆ cgelss()

subroutine cgelss ( integer m,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) s,
real rcond,
integer rank,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

CGELSS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> CGELSS computes the minimum norm solution to a complex linear
!> least squares problem:
!>
!> Minimize 2-norm(| b - A*x |).
!>
!> using the singular value decomposition (SVD) of A. A is an M-by-N
!> matrix which may be rank-deficient.
!>
!> 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.
!>
!> The effective rank of A is determined by treating as zero those
!> singular values which are less than RCOND times the largest singular
!> value.
!>
Parameters
[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, the first min(m,n) rows of A are overwritten with
!>          its right singular vectors, stored rowwise.
!>
[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 M-by-NRHS right hand side matrix B.
!>          On exit, B is overwritten by the N-by-NRHS solution matrix X.
!>          If m >= n and RANK = n, the residual sum-of-squares for
!>          the solution in the i-th column is given by the sum of
!>          squares of the modulus of elements n+1:m in that column.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,M,N).
!>
[out]S
!>          S is REAL array, dimension (min(M,N))
!>          The singular values of A in decreasing order.
!>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
!>
[in]RCOND
!>          RCOND is REAL
!>          RCOND is used to determine the effective rank of A.
!>          Singular values S(i) <= RCOND*S(1) are treated as zero.
!>          If RCOND < 0, machine precision is used instead.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the number of singular values
!>          which are greater than RCOND*S(1).
!>
[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 >= 1, and also:
!>          LWORK >=  2*min(M,N) + max(M,N,NRHS)
!>          For good performance, LWORK should generally be larger.
!>
!>          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]RWORK
!>          RWORK is REAL array, dimension (5*min(M,N))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  the algorithm for computing the SVD failed to converge;
!>                if INFO = i, i off-diagonal elements of an intermediate
!>                bidiagonal form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 176 of file cgelss.f.

178*
179* -- LAPACK driver routine --
180* -- LAPACK is a software package provided by Univ. of Tennessee, --
181* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182*
183* .. Scalar Arguments ..
184 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
185 REAL RCOND
186* ..
187* .. Array Arguments ..
188 REAL RWORK( * ), S( * )
189 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 REAL ZERO, ONE
196 parameter( zero = 0.0e+0, one = 1.0e+0 )
197 COMPLEX CZERO, CONE
198 parameter( czero = ( 0.0e+0, 0.0e+0 ),
199 $ cone = ( 1.0e+0, 0.0e+0 ) )
200* ..
201* .. Local Scalars ..
202 LOGICAL LQUERY
203 INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
204 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
205 $ MAXWRK, MINMN, MINWRK, MM, MNTHR
206 INTEGER LWORK_CGEQRF, LWORK_CUNMQR, LWORK_CGEBRD,
207 $ LWORK_CUNMBR, LWORK_CUNGBR, LWORK_CUNMLQ,
208 $ LWORK_CGELQF
209 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
210* ..
211* .. Local Arrays ..
212 COMPLEX DUM( 1 )
213* ..
214* .. External Subroutines ..
215 EXTERNAL cbdsqr, ccopy, cgebrd, cgelqf, cgemm, cgemv,
216 $ cgeqrf, clacpy, clascl, claset, csrscl, cungbr,
217 $ cunmbr, cunmlq, cunmqr, slabad, slascl, slaset,
218 $ xerbla
219* ..
220* .. External Functions ..
221 INTEGER ILAENV
222 REAL CLANGE, SLAMCH
223 EXTERNAL ilaenv, clange, slamch
224* ..
225* .. Intrinsic Functions ..
226 INTRINSIC max, min
227* ..
228* .. Executable Statements ..
229*
230* Test the input arguments
231*
232 info = 0
233 minmn = min( m, n )
234 maxmn = max( m, n )
235 lquery = ( lwork.EQ.-1 )
236 IF( m.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( nrhs.LT.0 ) THEN
241 info = -3
242 ELSE IF( lda.LT.max( 1, m ) ) THEN
243 info = -5
244 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
245 info = -7
246 END IF
247*
248* Compute workspace
249* (Note: Comments in the code beginning "Workspace:" describe the
250* minimal amount of workspace needed at that point in the code,
251* as well as the preferred amount for good performance.
252* CWorkspace refers to complex workspace, and RWorkspace refers
253* to real workspace. NB refers to the optimal block size for the
254* immediately following subroutine, as returned by ILAENV.)
255*
256 IF( info.EQ.0 ) THEN
257 minwrk = 1
258 maxwrk = 1
259 IF( minmn.GT.0 ) THEN
260 mm = m
261 mnthr = ilaenv( 6, 'CGELSS', ' ', m, n, nrhs, -1 )
262 IF( m.GE.n .AND. m.GE.mnthr ) THEN
263*
264* Path 1a - overdetermined, with many more rows than
265* columns
266*
267* Compute space needed for CGEQRF
268 CALL cgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
269 lwork_cgeqrf = real( dum(1) )
270* Compute space needed for CUNMQR
271 CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
272 $ ldb, dum(1), -1, info )
273 lwork_cunmqr = real( dum(1) )
274 mm = n
275 maxwrk = max( maxwrk, n + n*ilaenv( 1, 'CGEQRF', ' ', m,
276 $ n, -1, -1 ) )
277 maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'CUNMQR', 'LC',
278 $ m, nrhs, n, -1 ) )
279 END IF
280 IF( m.GE.n ) THEN
281*
282* Path 1 - overdetermined or exactly determined
283*
284* Compute space needed for CGEBRD
285 CALL cgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
286 $ -1, info )
287 lwork_cgebrd = real( dum(1) )
288* Compute space needed for CUNMBR
289 CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
290 $ b, ldb, dum(1), -1, info )
291 lwork_cunmbr = real( dum(1) )
292* Compute space needed for CUNGBR
293 CALL cungbr( 'P', n, n, n, a, lda, dum(1),
294 $ dum(1), -1, info )
295 lwork_cungbr = real( dum(1) )
296* Compute total workspace needed
297 maxwrk = max( maxwrk, 2*n + lwork_cgebrd )
298 maxwrk = max( maxwrk, 2*n + lwork_cunmbr )
299 maxwrk = max( maxwrk, 2*n + lwork_cungbr )
300 maxwrk = max( maxwrk, n*nrhs )
301 minwrk = 2*n + max( nrhs, m )
302 END IF
303 IF( n.GT.m ) THEN
304 minwrk = 2*m + max( nrhs, n )
305 IF( n.GE.mnthr ) THEN
306*
307* Path 2a - underdetermined, with many more columns
308* than rows
309*
310* Compute space needed for CGELQF
311 CALL cgelqf( m, n, a, lda, dum(1), dum(1),
312 $ -1, info )
313 lwork_cgelqf = real( dum(1) )
314* Compute space needed for CGEBRD
315 CALL cgebrd( m, m, a, lda, s, s, dum(1), dum(1),
316 $ dum(1), -1, info )
317 lwork_cgebrd = real( dum(1) )
318* Compute space needed for CUNMBR
319 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
320 $ dum(1), b, ldb, dum(1), -1, info )
321 lwork_cunmbr = real( dum(1) )
322* Compute space needed for CUNGBR
323 CALL cungbr( 'P', m, m, m, a, lda, dum(1),
324 $ dum(1), -1, info )
325 lwork_cungbr = real( dum(1) )
326* Compute space needed for CUNMLQ
327 CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
328 $ b, ldb, dum(1), -1, info )
329 lwork_cunmlq = real( dum(1) )
330* Compute total workspace needed
331 maxwrk = m + lwork_cgelqf
332 maxwrk = max( maxwrk, 3*m + m*m + lwork_cgebrd )
333 maxwrk = max( maxwrk, 3*m + m*m + lwork_cunmbr )
334 maxwrk = max( maxwrk, 3*m + m*m + lwork_cungbr )
335 IF( nrhs.GT.1 ) THEN
336 maxwrk = max( maxwrk, m*m + m + m*nrhs )
337 ELSE
338 maxwrk = max( maxwrk, m*m + 2*m )
339 END IF
340 maxwrk = max( maxwrk, m + lwork_cunmlq )
341 ELSE
342*
343* Path 2 - underdetermined
344*
345* Compute space needed for CGEBRD
346 CALL cgebrd( m, n, a, lda, s, s, dum(1), dum(1),
347 $ dum(1), -1, info )
348 lwork_cgebrd = real( dum(1) )
349* Compute space needed for CUNMBR
350 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
351 $ dum(1), b, ldb, dum(1), -1, info )
352 lwork_cunmbr = real( dum(1) )
353* Compute space needed for CUNGBR
354 CALL cungbr( 'P', m, n, m, a, lda, dum(1),
355 $ dum(1), -1, info )
356 lwork_cungbr = real( dum(1) )
357 maxwrk = 2*m + lwork_cgebrd
358 maxwrk = max( maxwrk, 2*m + lwork_cunmbr )
359 maxwrk = max( maxwrk, 2*m + lwork_cungbr )
360 maxwrk = max( maxwrk, n*nrhs )
361 END IF
362 END IF
363 maxwrk = max( minwrk, maxwrk )
364 END IF
365 work( 1 ) = maxwrk
366*
367 IF( lwork.LT.minwrk .AND. .NOT.lquery )
368 $ info = -12
369 END IF
370*
371 IF( info.NE.0 ) THEN
372 CALL xerbla( 'CGELSS', -info )
373 RETURN
374 ELSE IF( lquery ) THEN
375 RETURN
376 END IF
377*
378* Quick return if possible
379*
380 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
381 rank = 0
382 RETURN
383 END IF
384*
385* Get machine parameters
386*
387 eps = slamch( 'P' )
388 sfmin = slamch( 'S' )
389 smlnum = sfmin / eps
390 bignum = one / smlnum
391 CALL slabad( smlnum, bignum )
392*
393* Scale A if max element outside range [SMLNUM,BIGNUM]
394*
395 anrm = clange( 'M', m, n, a, lda, rwork )
396 iascl = 0
397 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
398*
399* Scale matrix norm up to SMLNUM
400*
401 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
402 iascl = 1
403 ELSE IF( anrm.GT.bignum ) THEN
404*
405* Scale matrix norm down to BIGNUM
406*
407 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
408 iascl = 2
409 ELSE IF( anrm.EQ.zero ) THEN
410*
411* Matrix all zero. Return zero solution.
412*
413 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
414 CALL slaset( 'F', minmn, 1, zero, zero, s, minmn )
415 rank = 0
416 GO TO 70
417 END IF
418*
419* Scale B if max element outside range [SMLNUM,BIGNUM]
420*
421 bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
422 ibscl = 0
423 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
424*
425* Scale matrix norm up to SMLNUM
426*
427 CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
428 ibscl = 1
429 ELSE IF( bnrm.GT.bignum ) THEN
430*
431* Scale matrix norm down to BIGNUM
432*
433 CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
434 ibscl = 2
435 END IF
436*
437* Overdetermined case
438*
439 IF( m.GE.n ) THEN
440*
441* Path 1 - overdetermined or exactly determined
442*
443 mm = m
444 IF( m.GE.mnthr ) THEN
445*
446* Path 1a - overdetermined, with many more rows than columns
447*
448 mm = n
449 itau = 1
450 iwork = itau + n
451*
452* Compute A=Q*R
453* (CWorkspace: need 2*N, prefer N+N*NB)
454* (RWorkspace: none)
455*
456 CALL cgeqrf( m, n, a, lda, work( itau ), work( iwork ),
457 $ lwork-iwork+1, info )
458*
459* Multiply B by transpose(Q)
460* (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
461* (RWorkspace: none)
462*
463 CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, work( itau ), b,
464 $ ldb, work( iwork ), lwork-iwork+1, info )
465*
466* Zero out below R
467*
468 IF( n.GT.1 )
469 $ CALL claset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
470 $ lda )
471 END IF
472*
473 ie = 1
474 itauq = 1
475 itaup = itauq + n
476 iwork = itaup + n
477*
478* Bidiagonalize R in A
479* (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
480* (RWorkspace: need N)
481*
482 CALL cgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
483 $ work( itaup ), work( iwork ), lwork-iwork+1,
484 $ info )
485*
486* Multiply B by transpose of left bidiagonalizing vectors of R
487* (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
488* (RWorkspace: none)
489*
490 CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, work( itauq ),
491 $ b, ldb, work( iwork ), lwork-iwork+1, info )
492*
493* Generate right bidiagonalizing vectors of R in A
494* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
495* (RWorkspace: none)
496*
497 CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
498 $ work( iwork ), lwork-iwork+1, info )
499 irwork = ie + n
500*
501* Perform bidiagonal QR iteration
502* multiply B by transpose of left singular vectors
503* compute right singular vectors in A
504* (CWorkspace: none)
505* (RWorkspace: need BDSPAC)
506*
507 CALL cbdsqr( 'U', n, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
508 $ 1, b, ldb, rwork( irwork ), info )
509 IF( info.NE.0 )
510 $ GO TO 70
511*
512* Multiply B by reciprocals of singular values
513*
514 thr = max( rcond*s( 1 ), sfmin )
515 IF( rcond.LT.zero )
516 $ thr = max( eps*s( 1 ), sfmin )
517 rank = 0
518 DO 10 i = 1, n
519 IF( s( i ).GT.thr ) THEN
520 CALL csrscl( nrhs, s( i ), b( i, 1 ), ldb )
521 rank = rank + 1
522 ELSE
523 CALL claset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
524 END IF
525 10 CONTINUE
526*
527* Multiply B by right singular vectors
528* (CWorkspace: need N, prefer N*NRHS)
529* (RWorkspace: none)
530*
531 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
532 CALL cgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
533 $ czero, work, ldb )
534 CALL clacpy( 'G', n, nrhs, work, ldb, b, ldb )
535 ELSE IF( nrhs.GT.1 ) THEN
536 chunk = lwork / n
537 DO 20 i = 1, nrhs, chunk
538 bl = min( nrhs-i+1, chunk )
539 CALL cgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
540 $ ldb, czero, work, n )
541 CALL clacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
542 20 CONTINUE
543 ELSE
544 CALL cgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
545 CALL ccopy( n, work, 1, b, 1 )
546 END IF
547*
548 ELSE IF( n.GE.mnthr .AND. lwork.GE.3*m+m*m+max( m, nrhs, n-2*m ) )
549 $ THEN
550*
551* Underdetermined case, M much less than N
552*
553* Path 2a - underdetermined, with many more columns than rows
554* and sufficient workspace for an efficient algorithm
555*
556 ldwork = m
557 IF( lwork.GE.3*m+m*lda+max( m, nrhs, n-2*m ) )
558 $ ldwork = lda
559 itau = 1
560 iwork = m + 1
561*
562* Compute A=L*Q
563* (CWorkspace: need 2*M, prefer M+M*NB)
564* (RWorkspace: none)
565*
566 CALL cgelqf( m, n, a, lda, work( itau ), work( iwork ),
567 $ lwork-iwork+1, info )
568 il = iwork
569*
570* Copy L to WORK(IL), zeroing out above it
571*
572 CALL clacpy( 'L', m, m, a, lda, work( il ), ldwork )
573 CALL claset( 'U', m-1, m-1, czero, czero, work( il+ldwork ),
574 $ ldwork )
575 ie = 1
576 itauq = il + ldwork*m
577 itaup = itauq + m
578 iwork = itaup + m
579*
580* Bidiagonalize L in WORK(IL)
581* (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
582* (RWorkspace: need M)
583*
584 CALL cgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
585 $ work( itauq ), work( itaup ), work( iwork ),
586 $ lwork-iwork+1, info )
587*
588* Multiply B by transpose of left bidiagonalizing vectors of L
589* (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
590* (RWorkspace: none)
591*
592 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, work( il ), ldwork,
593 $ work( itauq ), b, ldb, work( iwork ),
594 $ lwork-iwork+1, info )
595*
596* Generate right bidiagonalizing vectors of R in WORK(IL)
597* (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
598* (RWorkspace: none)
599*
600 CALL cungbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
601 $ work( iwork ), lwork-iwork+1, info )
602 irwork = ie + m
603*
604* Perform bidiagonal QR iteration, computing right singular
605* vectors of L in WORK(IL) and multiplying B by transpose of
606* left singular vectors
607* (CWorkspace: need M*M)
608* (RWorkspace: need BDSPAC)
609*
610 CALL cbdsqr( 'U', m, m, 0, nrhs, s, rwork( ie ), work( il ),
611 $ ldwork, a, lda, b, ldb, rwork( irwork ), info )
612 IF( info.NE.0 )
613 $ GO TO 70
614*
615* Multiply B by reciprocals of singular values
616*
617 thr = max( rcond*s( 1 ), sfmin )
618 IF( rcond.LT.zero )
619 $ thr = max( eps*s( 1 ), sfmin )
620 rank = 0
621 DO 30 i = 1, m
622 IF( s( i ).GT.thr ) THEN
623 CALL csrscl( nrhs, s( i ), b( i, 1 ), ldb )
624 rank = rank + 1
625 ELSE
626 CALL claset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
627 END IF
628 30 CONTINUE
629 iwork = il + m*ldwork
630*
631* Multiply B by right singular vectors of L in WORK(IL)
632* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
633* (RWorkspace: none)
634*
635 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
636 CALL cgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
637 $ b, ldb, czero, work( iwork ), ldb )
638 CALL clacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
639 ELSE IF( nrhs.GT.1 ) THEN
640 chunk = ( lwork-iwork+1 ) / m
641 DO 40 i = 1, nrhs, chunk
642 bl = min( nrhs-i+1, chunk )
643 CALL cgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
644 $ b( 1, i ), ldb, czero, work( iwork ), m )
645 CALL clacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
646 $ ldb )
647 40 CONTINUE
648 ELSE
649 CALL cgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
650 $ 1, czero, work( iwork ), 1 )
651 CALL ccopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
652 END IF
653*
654* Zero out below first M rows of B
655*
656 CALL claset( 'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
657 iwork = itau + m
658*
659* Multiply transpose(Q) by B
660* (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
661* (RWorkspace: none)
662*
663 CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, work( itau ), b,
664 $ ldb, work( iwork ), lwork-iwork+1, info )
665*
666 ELSE
667*
668* Path 2 - remaining underdetermined cases
669*
670 ie = 1
671 itauq = 1
672 itaup = itauq + m
673 iwork = itaup + m
674*
675* Bidiagonalize A
676* (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
677* (RWorkspace: need N)
678*
679 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
680 $ work( itaup ), work( iwork ), lwork-iwork+1,
681 $ info )
682*
683* Multiply B by transpose of left bidiagonalizing vectors
684* (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
685* (RWorkspace: none)
686*
687 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda, work( itauq ),
688 $ b, ldb, work( iwork ), lwork-iwork+1, info )
689*
690* Generate right bidiagonalizing vectors in A
691* (CWorkspace: need 3*M, prefer 2*M+M*NB)
692* (RWorkspace: none)
693*
694 CALL cungbr( 'P', m, n, m, a, lda, work( itaup ),
695 $ work( iwork ), lwork-iwork+1, info )
696 irwork = ie + m
697*
698* Perform bidiagonal QR iteration,
699* computing right singular vectors of A in A and
700* multiplying B by transpose of left singular vectors
701* (CWorkspace: none)
702* (RWorkspace: need BDSPAC)
703*
704 CALL cbdsqr( 'L', m, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
705 $ 1, b, ldb, rwork( irwork ), info )
706 IF( info.NE.0 )
707 $ GO TO 70
708*
709* Multiply B by reciprocals of singular values
710*
711 thr = max( rcond*s( 1 ), sfmin )
712 IF( rcond.LT.zero )
713 $ thr = max( eps*s( 1 ), sfmin )
714 rank = 0
715 DO 50 i = 1, m
716 IF( s( i ).GT.thr ) THEN
717 CALL csrscl( nrhs, s( i ), b( i, 1 ), ldb )
718 rank = rank + 1
719 ELSE
720 CALL claset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
721 END IF
722 50 CONTINUE
723*
724* Multiply B by right singular vectors of A
725* (CWorkspace: need N, prefer N*NRHS)
726* (RWorkspace: none)
727*
728 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
729 CALL cgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
730 $ czero, work, ldb )
731 CALL clacpy( 'G', n, nrhs, work, ldb, b, ldb )
732 ELSE IF( nrhs.GT.1 ) THEN
733 chunk = lwork / n
734 DO 60 i = 1, nrhs, chunk
735 bl = min( nrhs-i+1, chunk )
736 CALL cgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
737 $ ldb, czero, work, n )
738 CALL clacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
739 60 CONTINUE
740 ELSE
741 CALL cgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
742 CALL ccopy( n, work, 1, b, 1 )
743 END IF
744 END IF
745*
746* Undo scaling
747*
748 IF( iascl.EQ.1 ) THEN
749 CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
750 CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
751 $ info )
752 ELSE IF( iascl.EQ.2 ) THEN
753 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
754 CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
755 $ info )
756 END IF
757 IF( ibscl.EQ.1 ) THEN
758 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
759 ELSE IF( ibscl.EQ.2 ) THEN
760 CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
761 END IF
762 70 CONTINUE
763 work( 1 ) = maxwrk
764 RETURN
765*
766* End of CGELSS
767*
cungbr
subroutine cungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
CUNGBR
Definition cungbr.f:157
csrscl
subroutine csrscl(n, sa, sx, incx)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition csrscl.f:84
cbdsqr
subroutine cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
CBDSQR
Definition cbdsqr.f:222
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187

◆ cgelsx()

subroutine cgelsx ( integer m,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
real rcond,
integer rank,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CGELSX solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine CGELSY.
!>
!> CGELSX computes the minimum-norm solution to a complex linear least
!> squares problem:
!>     minimize || A * X - B ||
!> using a complete orthogonal factorization of A.  A is an M-by-N
!> matrix which may be rank-deficient.
!>
!> 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.
!>
!> The routine first computes a QR factorization with column pivoting:
!>     A * P = Q * [ R11 R12 ]
!>                 [  0  R22 ]
!> with R11 defined as the largest leading submatrix whose estimated
!> condition number is less than 1/RCOND.  The order of R11, RANK,
!> is the effective rank of A.
!>
!> Then, R22 is considered to be negligible, and R12 is annihilated
!> by unitary transformations from the right, arriving at the
!> complete orthogonal factorization:
!>    A * P = Q * [ T11 0 ] * Z
!>                [  0  0 ]
!> The minimum-norm solution is then
!>    X = P * Z**H [ inv(T11)*Q1**H*B ]
!>                 [        0         ]
!> where Q1 consists of the first RANK columns of Q.
!>
Parameters
[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 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 has been overwritten by details of its
!>          complete orthogonal factorization.
!>
[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 M-by-NRHS right hand side matrix B.
!>          On exit, the N-by-NRHS solution matrix X.
!>          If m >= n and RANK = n, the residual sum-of-squares for
!>          the solution in the i-th column is given by the sum of
!>          squares of elements N+1:M in that column.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,M,N).
!>
[in,out]JPVT
!>          JPVT is INTEGER array, dimension (N)
!>          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
!>          initial column, otherwise it is a free column.  Before
!>          the QR factorization of A, all initial columns are
!>          permuted to the leading positions; only the remaining
!>          free columns are moved as a result of column pivoting
!>          during the factorization.
!>          On exit, if JPVT(i) = k, then the i-th column of A*P
!>          was the k-th column of A.
!>
[in]RCOND
!>          RCOND is REAL
!>          RCOND is used to determine the effective rank of A, which
!>          is defined as the order of the largest leading triangular
!>          submatrix R11 in the QR factorization with pivoting of A,
!>          whose estimated condition number < 1/RCOND.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the order of the submatrix
!>          R11.  This is the same as the order of the submatrix T11
!>          in the complete orthogonal factorization of A.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension
!>                      (min(M,N) + max( N, 2*min(M,N)+NRHS )),
!>
[out]RWORK
!>          RWORK is REAL array, dimension (2*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 182 of file cgelsx.f.

184*
185* -- LAPACK driver routine --
186* -- LAPACK is a software package provided by Univ. of Tennessee, --
187* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*
189* .. Scalar Arguments ..
190 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
191 REAL RCOND
192* ..
193* .. Array Arguments ..
194 INTEGER JPVT( * )
195 REAL RWORK( * )
196 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 INTEGER IMAX, IMIN
203 parameter( imax = 1, imin = 2 )
204 REAL ZERO, ONE, DONE, NTDONE
205 parameter( zero = 0.0e+0, one = 1.0e+0, done = zero,
206 $ ntdone = one )
207 COMPLEX CZERO, CONE
208 parameter( czero = ( 0.0e+0, 0.0e+0 ),
209 $ cone = ( 1.0e+0, 0.0e+0 ) )
210* ..
211* .. Local Scalars ..
212 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
213 REAL ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
214 $ SMLNUM
215 COMPLEX C1, C2, S1, S2, T1, T2
216* ..
217* .. External Subroutines ..
218 EXTERNAL cgeqpf, claic1, clascl, claset, clatzm, ctrsm,
219 $ ctzrqf, cunm2r, slabad, xerbla
220* ..
221* .. External Functions ..
222 REAL CLANGE, SLAMCH
223 EXTERNAL clange, slamch
224* ..
225* .. Intrinsic Functions ..
226 INTRINSIC abs, conjg, max, min
227* ..
228* .. Executable Statements ..
229*
230 mn = min( m, n )
231 ismin = mn + 1
232 ismax = 2*mn + 1
233*
234* Test the input arguments.
235*
236 info = 0
237 IF( m.LT.0 ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( nrhs.LT.0 ) THEN
242 info = -3
243 ELSE IF( lda.LT.max( 1, m ) ) THEN
244 info = -5
245 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
246 info = -7
247 END IF
248*
249 IF( info.NE.0 ) THEN
250 CALL xerbla( 'CGELSX', -info )
251 RETURN
252 END IF
253*
254* Quick return if possible
255*
256 IF( min( m, n, nrhs ).EQ.0 ) THEN
257 rank = 0
258 RETURN
259 END IF
260*
261* Get machine parameters
262*
263 smlnum = slamch( 'S' ) / slamch( 'P' )
264 bignum = one / smlnum
265 CALL slabad( smlnum, bignum )
266*
267* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
268*
269 anrm = clange( 'M', m, n, a, lda, rwork )
270 iascl = 0
271 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
272*
273* Scale matrix norm up to SMLNUM
274*
275 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
276 iascl = 1
277 ELSE IF( anrm.GT.bignum ) THEN
278*
279* Scale matrix norm down to BIGNUM
280*
281 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
282 iascl = 2
283 ELSE IF( anrm.EQ.zero ) THEN
284*
285* Matrix all zero. Return zero solution.
286*
287 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
288 rank = 0
289 GO TO 100
290 END IF
291*
292 bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
293 ibscl = 0
294 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
295*
296* Scale matrix norm up to SMLNUM
297*
298 CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
299 ibscl = 1
300 ELSE IF( bnrm.GT.bignum ) THEN
301*
302* Scale matrix norm down to BIGNUM
303*
304 CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
305 ibscl = 2
306 END IF
307*
308* Compute QR factorization with column pivoting of A:
309* A * P = Q * R
310*
311 CALL cgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), rwork,
312 $ info )
313*
314* complex workspace MN+N. Real workspace 2*N. Details of Householder
315* rotations stored in WORK(1:MN).
316*
317* Determine RANK using incremental condition estimation
318*
319 work( ismin ) = cone
320 work( ismax ) = cone
321 smax = abs( a( 1, 1 ) )
322 smin = smax
323 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
324 rank = 0
325 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
326 GO TO 100
327 ELSE
328 rank = 1
329 END IF
330*
331 10 CONTINUE
332 IF( rank.LT.mn ) THEN
333 i = rank + 1
334 CALL claic1( imin, rank, work( ismin ), smin, a( 1, i ),
335 $ a( i, i ), sminpr, s1, c1 )
336 CALL claic1( imax, rank, work( ismax ), smax, a( 1, i ),
337 $ a( i, i ), smaxpr, s2, c2 )
338*
339 IF( smaxpr*rcond.LE.sminpr ) THEN
340 DO 20 i = 1, rank
341 work( ismin+i-1 ) = s1*work( ismin+i-1 )
342 work( ismax+i-1 ) = s2*work( ismax+i-1 )
343 20 CONTINUE
344 work( ismin+rank ) = c1
345 work( ismax+rank ) = c2
346 smin = sminpr
347 smax = smaxpr
348 rank = rank + 1
349 GO TO 10
350 END IF
351 END IF
352*
353* Logically partition R = [ R11 R12 ]
354* [ 0 R22 ]
355* where R11 = R(1:RANK,1:RANK)
356*
357* [R11,R12] = [ T11, 0 ] * Y
358*
359 IF( rank.LT.n )
360 $ CALL ctzrqf( rank, n, a, lda, work( mn+1 ), info )
361*
362* Details of Householder rotations stored in WORK(MN+1:2*MN)
363*
364* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
365*
366 CALL cunm2r( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
367 $ work( 1 ), b, ldb, work( 2*mn+1 ), info )
368*
369* workspace NRHS
370*
371* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
372*
373 CALL ctrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
374 $ nrhs, cone, a, lda, b, ldb )
375*
376 DO 40 i = rank + 1, n
377 DO 30 j = 1, nrhs
378 b( i, j ) = czero
379 30 CONTINUE
380 40 CONTINUE
381*
382* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
383*
384 IF( rank.LT.n ) THEN
385 DO 50 i = 1, rank
386 CALL clatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
387 $ conjg( work( mn+i ) ), b( i, 1 ),
388 $ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
389 50 CONTINUE
390 END IF
391*
392* workspace NRHS
393*
394* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
395*
396 DO 90 j = 1, nrhs
397 DO 60 i = 1, n
398 work( 2*mn+i ) = ntdone
399 60 CONTINUE
400 DO 80 i = 1, n
401 IF( work( 2*mn+i ).EQ.ntdone ) THEN
402 IF( jpvt( i ).NE.i ) THEN
403 k = i
404 t1 = b( k, j )
405 t2 = b( jpvt( k ), j )
406 70 CONTINUE
407 b( jpvt( k ), j ) = t1
408 work( 2*mn+k ) = done
409 t1 = t2
410 k = jpvt( k )
411 t2 = b( jpvt( k ), j )
412 IF( jpvt( k ).NE.i )
413 $ GO TO 70
414 b( i, j ) = t1
415 work( 2*mn+k ) = done
416 END IF
417 END IF
418 80 CONTINUE
419 90 CONTINUE
420*
421* Undo scaling
422*
423 IF( iascl.EQ.1 ) THEN
424 CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
425 CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
426 $ info )
427 ELSE IF( iascl.EQ.2 ) THEN
428 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
429 CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
430 $ info )
431 END IF
432 IF( ibscl.EQ.1 ) THEN
433 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
434 ELSE IF( ibscl.EQ.2 ) THEN
435 CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
436 END IF
437*
438 100 CONTINUE
439*
440 RETURN
441*
442* End of CGELSX
443*
cgeqpf
subroutine cgeqpf(m, n, a, lda, jpvt, tau, work, rwork, info)
CGEQPF
Definition cgeqpf.f:148
claic1
subroutine claic1(job, j, x, sest, w, gamma, sestpr, s, c)
CLAIC1 applies one step of incremental condition estimation.
Definition claic1.f:135
ctzrqf
subroutine ctzrqf(m, n, a, lda, tau, info)
CTZRQF
Definition ctzrqf.f:138
clatzm
subroutine clatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
CLATZM
Definition clatzm.f:152
cunm2r
subroutine cunm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition cunm2r.f:159
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180

◆ cgelsy()

subroutine cgelsy ( integer m,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
real rcond,
integer rank,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

CGELSY solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> CGELSY computes the minimum-norm solution to a complex linear least
!> squares problem:
!>     minimize || A * X - B ||
!> using a complete orthogonal factorization of A.  A is an M-by-N
!> matrix which may be rank-deficient.
!>
!> 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.
!>
!> The routine first computes a QR factorization with column pivoting:
!>     A * P = Q * [ R11 R12 ]
!>                 [  0  R22 ]
!> with R11 defined as the largest leading submatrix whose estimated
!> condition number is less than 1/RCOND.  The order of R11, RANK,
!> is the effective rank of A.
!>
!> Then, R22 is considered to be negligible, and R12 is annihilated
!> by unitary transformations from the right, arriving at the
!> complete orthogonal factorization:
!>    A * P = Q * [ T11 0 ] * Z
!>                [  0  0 ]
!> The minimum-norm solution is then
!>    X = P * Z**H [ inv(T11)*Q1**H*B ]
!>                 [        0         ]
!> where Q1 consists of the first RANK columns of Q.
!>
!> This routine is basically identical to the original xGELSX except
!> three differences:
!>   o The permutation of matrix B (the right hand side) is faster and
!>     more simple.
!>   o The call to the subroutine xGEQPF has been substituted by the
!>     the call to the subroutine xGEQP3. This subroutine is a Blas-3
!>     version of the QR factorization with column pivoting.
!>   o Matrix B (the right hand side) is updated with Blas-3.
!>
Parameters
[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 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 has been overwritten by details of its
!>          complete orthogonal factorization.
!>
[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 M-by-NRHS right hand side matrix B.
!>          On exit, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,M,N).
!>
[in,out]JPVT
!>          JPVT is INTEGER array, dimension (N)
!>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
!>          to the front of AP, otherwise column i is a free column.
!>          On exit, if JPVT(i) = k, then the i-th column of A*P
!>          was the k-th column of A.
!>
[in]RCOND
!>          RCOND is REAL
!>          RCOND is used to determine the effective rank of A, which
!>          is defined as the order of the largest leading triangular
!>          submatrix R11 in the QR factorization with pivoting of A,
!>          whose estimated condition number < 1/RCOND.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the order of the submatrix
!>          R11.  This is the same as the order of the submatrix T11
!>          in the complete orthogonal factorization of A.
!>
[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.
!>          The unblocked strategy requires that:
!>            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
!>          where MN = min(M,N).
!>          The block algorithm requires that:
!>            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
!>          where NB is an upper bound on the blocksize returned
!>          by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR,
!>          and CUNMRZ.
!>
!>          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]RWORK
!>          RWORK is REAL array, dimension (2*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

Definition at line 208 of file cgelsy.f.

210*
211* -- LAPACK driver routine --
212* -- LAPACK is a software package provided by Univ. of Tennessee, --
213* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214*
215* .. Scalar Arguments ..
216 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
217 REAL RCOND
218* ..
219* .. Array Arguments ..
220 INTEGER JPVT( * )
221 REAL RWORK( * )
222 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228 INTEGER IMAX, IMIN
229 parameter( imax = 1, imin = 2 )
230 REAL ZERO, ONE
231 parameter( zero = 0.0e+0, one = 1.0e+0 )
232 COMPLEX CZERO, CONE
233 parameter( czero = ( 0.0e+0, 0.0e+0 ),
234 $ cone = ( 1.0e+0, 0.0e+0 ) )
235* ..
236* .. Local Scalars ..
237 LOGICAL LQUERY
238 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
239 $ NB, NB1, NB2, NB3, NB4
240 REAL ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
241 $ SMLNUM, WSIZE
242 COMPLEX C1, C2, S1, S2
243* ..
244* .. External Subroutines ..
245 EXTERNAL ccopy, cgeqp3, claic1, clascl, claset, ctrsm,
246 $ ctzrzf, cunmqr, cunmrz, slabad, xerbla
247* ..
248* .. External Functions ..
249 INTEGER ILAENV
250 REAL CLANGE, SLAMCH
251 EXTERNAL clange, ilaenv, slamch
252* ..
253* .. Intrinsic Functions ..
254 INTRINSIC abs, max, min, real, cmplx
255* ..
256* .. Executable Statements ..
257*
258 mn = min( m, n )
259 ismin = mn + 1
260 ismax = 2*mn + 1
261*
262* Test the input arguments.
263*
264 info = 0
265 nb1 = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
266 nb2 = ilaenv( 1, 'CGERQF', ' ', m, n, -1, -1 )
267 nb3 = ilaenv( 1, 'CUNMQR', ' ', m, n, nrhs, -1 )
268 nb4 = ilaenv( 1, 'CUNMRQ', ' ', m, n, nrhs, -1 )
269 nb = max( nb1, nb2, nb3, nb4 )
270 lwkopt = max( 1, mn+2*n+nb*(n+1), 2*mn+nb*nrhs )
271 work( 1 ) = cmplx( lwkopt )
272 lquery = ( lwork.EQ.-1 )
273 IF( m.LT.0 ) THEN
274 info = -1
275 ELSE IF( n.LT.0 ) THEN
276 info = -2
277 ELSE IF( nrhs.LT.0 ) THEN
278 info = -3
279 ELSE IF( lda.LT.max( 1, m ) ) THEN
280 info = -5
281 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
282 info = -7
283 ELSE IF( lwork.LT.( mn+max( 2*mn, n+1, mn+nrhs ) ) .AND.
284 $ .NOT.lquery ) THEN
285 info = -12
286 END IF
287*
288 IF( info.NE.0 ) THEN
289 CALL xerbla( 'CGELSY', -info )
290 RETURN
291 ELSE IF( lquery ) THEN
292 RETURN
293 END IF
294*
295* Quick return if possible
296*
297 IF( min( m, n, nrhs ).EQ.0 ) THEN
298 rank = 0
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 entries outside range [SMLNUM,BIGNUM]
309*
310 anrm = clange( 'M', m, n, a, lda, rwork )
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', max( m, n ), nrhs, czero, czero, b, ldb )
329 rank = 0
330 GO TO 70
331 END IF
332*
333 bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
334 ibscl = 0
335 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
336*
337* Scale matrix norm up to SMLNUM
338*
339 CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
340 ibscl = 1
341 ELSE IF( bnrm.GT.bignum ) THEN
342*
343* Scale matrix norm down to BIGNUM
344*
345 CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
346 ibscl = 2
347 END IF
348*
349* Compute QR factorization with column pivoting of A:
350* A * P = Q * R
351*
352 CALL cgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
353 $ lwork-mn, rwork, info )
354 wsize = mn + real( work( mn+1 ) )
355*
356* complex workspace: MN+NB*(N+1). real workspace 2*N.
357* Details of Householder rotations stored in WORK(1:MN).
358*
359* Determine RANK using incremental condition estimation
360*
361 work( ismin ) = cone
362 work( ismax ) = cone
363 smax = abs( a( 1, 1 ) )
364 smin = smax
365 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
366 rank = 0
367 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
368 GO TO 70
369 ELSE
370 rank = 1
371 END IF
372*
373 10 CONTINUE
374 IF( rank.LT.mn ) THEN
375 i = rank + 1
376 CALL claic1( imin, rank, work( ismin ), smin, a( 1, i ),
377 $ a( i, i ), sminpr, s1, c1 )
378 CALL claic1( imax, rank, work( ismax ), smax, a( 1, i ),
379 $ a( i, i ), smaxpr, s2, c2 )
380*
381 IF( smaxpr*rcond.LE.sminpr ) THEN
382 DO 20 i = 1, rank
383 work( ismin+i-1 ) = s1*work( ismin+i-1 )
384 work( ismax+i-1 ) = s2*work( ismax+i-1 )
385 20 CONTINUE
386 work( ismin+rank ) = c1
387 work( ismax+rank ) = c2
388 smin = sminpr
389 smax = smaxpr
390 rank = rank + 1
391 GO TO 10
392 END IF
393 END IF
394*
395* complex workspace: 3*MN.
396*
397* Logically partition R = [ R11 R12 ]
398* [ 0 R22 ]
399* where R11 = R(1:RANK,1:RANK)
400*
401* [R11,R12] = [ T11, 0 ] * Y
402*
403 IF( rank.LT.n )
404 $ CALL ctzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
405 $ lwork-2*mn, info )
406*
407* complex workspace: 2*MN.
408* Details of Householder rotations stored in WORK(MN+1:2*MN)
409*
410* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
411*
412 CALL cunmqr( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
413 $ work( 1 ), b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
414 wsize = max( wsize, 2*mn+real( work( 2*mn+1 ) ) )
415*
416* complex workspace: 2*MN+NB*NRHS.
417*
418* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
419*
420 CALL ctrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
421 $ nrhs, cone, a, lda, b, ldb )
422*
423 DO 40 j = 1, nrhs
424 DO 30 i = rank + 1, n
425 b( i, j ) = czero
426 30 CONTINUE
427 40 CONTINUE
428*
429* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
430*
431 IF( rank.LT.n ) THEN
432 CALL cunmrz( 'Left', 'Conjugate transpose', n, nrhs, rank,
433 $ n-rank, a, lda, work( mn+1 ), b, ldb,
434 $ work( 2*mn+1 ), lwork-2*mn, info )
435 END IF
436*
437* complex workspace: 2*MN+NRHS.
438*
439* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
440*
441 DO 60 j = 1, nrhs
442 DO 50 i = 1, n
443 work( jpvt( i ) ) = b( i, j )
444 50 CONTINUE
445 CALL ccopy( n, work( 1 ), 1, b( 1, j ), 1 )
446 60 CONTINUE
447*
448* complex workspace: N.
449*
450* Undo scaling
451*
452 IF( iascl.EQ.1 ) THEN
453 CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
454 CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
455 $ info )
456 ELSE IF( iascl.EQ.2 ) THEN
457 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
458 CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
459 $ info )
460 END IF
461 IF( ibscl.EQ.1 ) THEN
462 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
463 ELSE IF( ibscl.EQ.2 ) THEN
464 CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
465 END IF
466*
467 70 CONTINUE
468 work( 1 ) = cmplx( lwkopt )
469*
470 RETURN
471*
472* End of CGELSY
473*
cmplx
float cmplx[2]
Definition pblas.h:136
cgeqp3
subroutine cgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
CGEQP3
Definition cgeqp3.f:159
cunmrz
subroutine cunmrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
CUNMRZ
Definition cunmrz.f:187
ctzrzf
subroutine ctzrzf(m, n, a, lda, tau, work, lwork, info)
CTZRZF
Definition ctzrzf.f:151

◆ cgesv()

subroutine cgesv ( integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
complex, dimension( ldb, * ) b,
integer ldb,
integer info )

CGESV computes the solution to system of linear equations A * X = B for GE matrices (simple driver)

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

Purpose:
!>
!> CGESV computes the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!> The LU decomposition with partial pivoting and row interchanges is
!> used to factor A as
!>    A = P * L * U,
!> where P is a permutation matrix, L is unit lower triangular, and U is
!> upper triangular.  The factored form of A is then used to solve the
!> system of equations A * X = B.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order 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 matrix B.  NRHS >= 0.
!>
[in,out]A
!>          A is COMPLEX array, dimension (LDA,N)
!>          On entry, the N-by-N coefficient matrix A.
!>          On exit, the factors L and U from the factorization
!>          A = P*L*U; the unit diagonal elements of L are not stored.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices that define the permutation matrix P;
!>          row i of the matrix was interchanged with row IPIV(i).
!>
[in,out]B
!>          B is COMPLEX array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS matrix of right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
!>                has been completed, but the factor U is exactly
!>                singular, so the solution could not be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 121 of file cgesv.f.

122*
123* -- LAPACK driver routine --
124* -- LAPACK is a software package provided by Univ. of Tennessee, --
125* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126*
127* .. Scalar Arguments ..
128 INTEGER INFO, LDA, LDB, N, NRHS
129* ..
130* .. Array Arguments ..
131 INTEGER IPIV( * )
132 COMPLEX A( LDA, * ), B( LDB, * )
133* ..
134*
135* =====================================================================
136*
137* .. External Subroutines ..
138 EXTERNAL cgetrf, cgetrs, xerbla
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC max
142* ..
143* .. Executable Statements ..
144*
145* Test the input parameters.
146*
147 info = 0
148 IF( n.LT.0 ) THEN
149 info = -1
150 ELSE IF( nrhs.LT.0 ) THEN
151 info = -2
152 ELSE IF( lda.LT.max( 1, n ) ) THEN
153 info = -4
154 ELSE IF( ldb.LT.max( 1, n ) ) THEN
155 info = -7
156 END IF
157 IF( info.NE.0 ) THEN
158 CALL xerbla( 'CGESV ', -info )
159 RETURN
160 END IF
161*
162* Compute the LU factorization of A.
163*
164 CALL cgetrf( n, n, a, lda, ipiv, info )
165 IF( info.EQ.0 ) THEN
166*
167* Solve the system A*X = B, overwriting B with X.
168*
169 CALL cgetrs( 'No transpose', n, nrhs, a, lda, ipiv, b, ldb,
170 $ info )
171 END IF
172 RETURN
173*
174* End of CGESV
175*
cgetrs
subroutine cgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
CGETRS
Definition cgetrs.f:121
cgetrf
subroutine cgetrf(m, n, a, lda, ipiv, info)
CGETRF
Definition cgetrf.f:108

◆ cgesvx()

subroutine cgesvx ( character fact,
character trans,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldaf, * ) af,
integer ldaf,
integer, dimension( * ) ipiv,
character equed,
real, dimension( * ) r,
real, dimension( * ) c,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( ldx, * ) x,
integer ldx,
real rcond,
real, dimension( * ) ferr,
real, dimension( * ) berr,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CGESVX computes the solution to system of linear equations A * X = B for GE matrices

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

Purpose:
!>
!> CGESVX uses the LU factorization to compute the solution to a complex
!> system of linear equations
!>    A * X = B,
!> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>       TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
!>       TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
!>       TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
!>    or diag(C)*B (if TRANS = 'T' or 'C').
!>
!> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
!>    matrix A (after equilibration if FACT = 'E') as
!>       A = P * L * U,
!>    where P is a permutation matrix, L is a unit lower triangular
!>    matrix, and U is upper triangular.
!>
!> 3. If some U(i,i)=0, so that U is exactly singular, then the routine
!>    returns with INFO = i. Otherwise, the factored form of A is used
!>    to estimate the condition number of the matrix A.  If the
!>    reciprocal of the condition number is less than machine precision,
!>    INFO = N+1 is returned as a warning, but the routine still goes on
!>    to solve for X and compute error bounds as described below.
!>
!> 4. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 5. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
!> 6. If equilibration was used, the matrix X is premultiplied by
!>    diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
!>    that it solves the original system before equilibration.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AF and IPIV contain the factored form of A.
!>                  If EQUED is not 'N', the matrix A has been
!>                  equilibrated with scaling factors given by R and C.
!>                  A, AF, and IPIV are not modified.
!>          = 'N':  The matrix A will be copied to AF and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AF and factored.
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations:
!>          = 'N':  A * X = B     (No transpose)
!>          = 'T':  A**T * X = B  (Transpose)
!>          = 'C':  A**H * X = B  (Conjugate transpose)
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order 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 N-by-N matrix A.  If FACT = 'F' and EQUED is
!>          not 'N', then A must have been equilibrated by the scaling
!>          factors in R and/or C.  A is not modified if FACT = 'F' or
!>          'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>          On exit, if EQUED .ne. 'N', A is scaled as follows:
!>          EQUED = 'R':  A := diag(R) * A
!>          EQUED = 'C':  A := A * diag(C)
!>          EQUED = 'B':  A := diag(R) * A * diag(C).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]AF
!>          AF is COMPLEX array, dimension (LDAF,N)
!>          If FACT = 'F', then AF is an input argument and on entry
!>          contains the factors L and U from the factorization
!>          A = P*L*U as computed by CGETRF.  If EQUED .ne. 'N', then
!>          AF is the factored form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AF is an output argument and on exit
!>          returns the factors L and U from the factorization A = P*L*U
!>          of the original matrix A.
!>
!>          If FACT = 'E', then AF is an output argument and on exit
!>          returns the factors L and U from the factorization A = P*L*U
!>          of the equilibrated matrix A (see the description of A for
!>          the form of the equilibrated matrix).
!>
[in]LDAF
!>          LDAF is INTEGER
!>          The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          If FACT = 'F', then IPIV is an input argument and on entry
!>          contains the pivot indices from the factorization A = P*L*U
!>          as computed by CGETRF; row i of the matrix was interchanged
!>          with row IPIV(i).
!>
!>          If FACT = 'N', then IPIV is an output argument and on exit
!>          contains the pivot indices from the factorization A = P*L*U
!>          of the original matrix A.
!>
!>          If FACT = 'E', then IPIV is an output argument and on exit
!>          contains the pivot indices from the factorization A = P*L*U
!>          of the equilibrated matrix A.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'R':  Row equilibration, i.e., A has been premultiplied by
!>                  diag(R).
!>          = 'C':  Column equilibration, i.e., A has been postmultiplied
!>                  by diag(C).
!>          = 'B':  Both row and column equilibration, i.e., A has been
!>                  replaced by diag(R) * A * diag(C).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]R
!>          R is REAL array, dimension (N)
!>          The row scale factors for A.  If EQUED = 'R' or 'B', A is
!>          multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
!>          is not accessed.  R is an input argument if FACT = 'F';
!>          otherwise, R is an output argument.  If FACT = 'F' and
!>          EQUED = 'R' or 'B', each element of R must be positive.
!>
[in,out]C
!>          C is REAL array, dimension (N)
!>          The column scale factors for A.  If EQUED = 'C' or 'B', A is
!>          multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
!>          is not accessed.  C is an input argument if FACT = 'F';
!>          otherwise, C is an output argument.  If FACT = 'F' and
!>          EQUED = 'C' or 'B', each element of C must be positive.
!>
[in,out]B
!>          B is COMPLEX array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit,
!>          if EQUED = 'N', B is not modified;
!>          if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
!>          diag(R)*B;
!>          if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
!>          overwritten by diag(C)*B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is COMPLEX array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
!>          to the original system of equations.  Note that A and B are
!>          modified on exit if EQUED .ne. 'N', and the solution to the
!>          equilibrated system is inv(diag(C))*X if TRANS = 'N' and
!>          EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
!>          and EQUED = 'R' or 'B'.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>
[out]FERR
!>          FERR is REAL array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (2*N)
!>          On exit, RWORK(1) contains the reciprocal pivot growth
!>          factor norm(A)/norm(U). The  norm is
!>          used. If RWORK(1) is much less than 1, then the stability
!>          of the LU factorization of the (equilibrated) matrix A
!>          could be poor. This also means that the solution X, condition
!>          estimator RCOND, and forward error bound FERR could be
!>          unreliable. If factorization fails with 0<INFO<=N, then
!>          RWORK(1) contains the reciprocal pivot growth factor for the
!>          leading INFO columns of A.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  U(i,i) is exactly zero.  The factorization has
!>                       been completed, but the factor U is exactly
!>                       singular, so the solution and error bounds
!>                       could not be computed. RCOND = 0 is returned.
!>                = N+1: U is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 347 of file cgesvx.f.

350*
351* -- LAPACK driver routine --
352* -- LAPACK is a software package provided by Univ. of Tennessee, --
353* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
354*
355* .. Scalar Arguments ..
356 CHARACTER EQUED, FACT, TRANS
357 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
358 REAL RCOND
359* ..
360* .. Array Arguments ..
361 INTEGER IPIV( * )
362 REAL BERR( * ), C( * ), FERR( * ), R( * ),
363 $ RWORK( * )
364 COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
365 $ WORK( * ), X( LDX, * )
366* ..
367*
368* =====================================================================
369*
370* .. Parameters ..
371 REAL ZERO, ONE
372 parameter( zero = 0.0e+0, one = 1.0e+0 )
373* ..
374* .. Local Scalars ..
375 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
376 CHARACTER NORM
377 INTEGER I, INFEQU, J
378 REAL AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
379 $ ROWCND, RPVGRW, SMLNUM
380* ..
381* .. External Functions ..
382 LOGICAL LSAME
383 REAL CLANGE, CLANTR, SLAMCH
384 EXTERNAL lsame, clange, clantr, slamch
385* ..
386* .. External Subroutines ..
387 EXTERNAL cgecon, cgeequ, cgerfs, cgetrf, cgetrs, clacpy,
388 $ claqge, xerbla
389* ..
390* .. Intrinsic Functions ..
391 INTRINSIC max, min
392* ..
393* .. Executable Statements ..
394*
395 info = 0
396 nofact = lsame( fact, 'N' )
397 equil = lsame( fact, 'E' )
398 notran = lsame( trans, 'N' )
399 IF( nofact .OR. equil ) THEN
400 equed = 'N'
401 rowequ = .false.
402 colequ = .false.
403 ELSE
404 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
405 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
406 smlnum = slamch( 'Safe minimum' )
407 bignum = one / smlnum
408 END IF
409*
410* Test the input parameters.
411*
412 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
413 $ THEN
414 info = -1
415 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
416 $ lsame( trans, 'C' ) ) THEN
417 info = -2
418 ELSE IF( n.LT.0 ) THEN
419 info = -3
420 ELSE IF( nrhs.LT.0 ) THEN
421 info = -4
422 ELSE IF( lda.LT.max( 1, n ) ) THEN
423 info = -6
424 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
425 info = -8
426 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
427 $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
428 info = -10
429 ELSE
430 IF( rowequ ) THEN
431 rcmin = bignum
432 rcmax = zero
433 DO 10 j = 1, n
434 rcmin = min( rcmin, r( j ) )
435 rcmax = max( rcmax, r( j ) )
436 10 CONTINUE
437 IF( rcmin.LE.zero ) THEN
438 info = -11
439 ELSE IF( n.GT.0 ) THEN
440 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
441 ELSE
442 rowcnd = one
443 END IF
444 END IF
445 IF( colequ .AND. info.EQ.0 ) THEN
446 rcmin = bignum
447 rcmax = zero
448 DO 20 j = 1, n
449 rcmin = min( rcmin, c( j ) )
450 rcmax = max( rcmax, c( j ) )
451 20 CONTINUE
452 IF( rcmin.LE.zero ) THEN
453 info = -12
454 ELSE IF( n.GT.0 ) THEN
455 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
456 ELSE
457 colcnd = one
458 END IF
459 END IF
460 IF( info.EQ.0 ) THEN
461 IF( ldb.LT.max( 1, n ) ) THEN
462 info = -14
463 ELSE IF( ldx.LT.max( 1, n ) ) THEN
464 info = -16
465 END IF
466 END IF
467 END IF
468*
469 IF( info.NE.0 ) THEN
470 CALL xerbla( 'CGESVX', -info )
471 RETURN
472 END IF
473*
474 IF( equil ) THEN
475*
476* Compute row and column scalings to equilibrate the matrix A.
477*
478 CALL cgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
479 IF( infequ.EQ.0 ) THEN
480*
481* Equilibrate the matrix.
482*
483 CALL claqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
484 $ equed )
485 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
486 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
487 END IF
488 END IF
489*
490* Scale the right hand side.
491*
492 IF( notran ) THEN
493 IF( rowequ ) THEN
494 DO 40 j = 1, nrhs
495 DO 30 i = 1, n
496 b( i, j ) = r( i )*b( i, j )
497 30 CONTINUE
498 40 CONTINUE
499 END IF
500 ELSE IF( colequ ) THEN
501 DO 60 j = 1, nrhs
502 DO 50 i = 1, n
503 b( i, j ) = c( i )*b( i, j )
504 50 CONTINUE
505 60 CONTINUE
506 END IF
507*
508 IF( nofact .OR. equil ) THEN
509*
510* Compute the LU factorization of A.
511*
512 CALL clacpy( 'Full', n, n, a, lda, af, ldaf )
513 CALL cgetrf( n, n, af, ldaf, ipiv, info )
514*
515* Return if INFO is non-zero.
516*
517 IF( info.GT.0 ) THEN
518*
519* Compute the reciprocal pivot growth factor of the
520* leading rank-deficient INFO columns of A.
521*
522 rpvgrw = clantr( 'M', 'U', 'N', info, info, af, ldaf,
523 $ rwork )
524 IF( rpvgrw.EQ.zero ) THEN
525 rpvgrw = one
526 ELSE
527 rpvgrw = clange( 'M', n, info, a, lda, rwork ) /
528 $ rpvgrw
529 END IF
530 rwork( 1 ) = rpvgrw
531 rcond = zero
532 RETURN
533 END IF
534 END IF
535*
536* Compute the norm of the matrix A and the
537* reciprocal pivot growth factor RPVGRW.
538*
539 IF( notran ) THEN
540 norm = '1'
541 ELSE
542 norm = 'I'
543 END IF
544 anorm = clange( norm, n, n, a, lda, rwork )
545 rpvgrw = clantr( 'M', 'U', 'N', n, n, af, ldaf, rwork )
546 IF( rpvgrw.EQ.zero ) THEN
547 rpvgrw = one
548 ELSE
549 rpvgrw = clange( 'M', n, n, a, lda, rwork ) / rpvgrw
550 END IF
551*
552* Compute the reciprocal of the condition number of A.
553*
554 CALL cgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, info )
555*
556* Compute the solution matrix X.
557*
558 CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
559 CALL cgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
560*
561* Use iterative refinement to improve the computed solution and
562* compute error bounds and backward error estimates for it.
563*
564 CALL cgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
565 $ ldx, ferr, berr, work, rwork, info )
566*
567* Transform the solution matrix X to a solution of the original
568* system.
569*
570 IF( notran ) THEN
571 IF( colequ ) THEN
572 DO 80 j = 1, nrhs
573 DO 70 i = 1, n
574 x( i, j ) = c( i )*x( i, j )
575 70 CONTINUE
576 80 CONTINUE
577 DO 90 j = 1, nrhs
578 ferr( j ) = ferr( j ) / colcnd
579 90 CONTINUE
580 END IF
581 ELSE IF( rowequ ) THEN
582 DO 110 j = 1, nrhs
583 DO 100 i = 1, n
584 x( i, j ) = r( i )*x( i, j )
585 100 CONTINUE
586 110 CONTINUE
587 DO 120 j = 1, nrhs
588 ferr( j ) = ferr( j ) / rowcnd
589 120 CONTINUE
590 END IF
591*
592* Set INFO = N+1 if the matrix is singular to working precision.
593*
594 IF( rcond.LT.slamch( 'Epsilon' ) )
595 $ info = n + 1
596*
597 rwork( 1 ) = rpvgrw
598 RETURN
599*
600* End of CGESVX
601*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
claqge
subroutine claqge(m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
CLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
Definition claqge.f:143
cgeequ
subroutine cgeequ(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
CGEEQU
Definition cgeequ.f:140
cgecon
subroutine cgecon(norm, n, a, lda, anorm, rcond, work, rwork, info)
CGECON
Definition cgecon.f:124
cgerfs
subroutine cgerfs(trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CGERFS
Definition cgerfs.f:186
clantr
real function clantr(norm, uplo, diag, m, n, a, lda, work)
CLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clantr.f:142

◆ cgesvxx()

subroutine cgesvxx ( character fact,
character trans,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldaf, * ) af,
integer ldaf,
integer, dimension( * ) ipiv,
character equed,
real, dimension( * ) r,
real, dimension( * ) c,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( ldx , * ) x,
integer ldx,
real rcond,
real rpvgrw,
real, dimension( * ) berr,
integer n_err_bnds,
real, dimension( nrhs, * ) err_bnds_norm,
real, dimension( nrhs, * ) err_bnds_comp,
integer nparams,
real, dimension( * ) params,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CGESVXX computes the solution to system of linear equations A * X = B for GE matrices

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

Purpose:
!>
!>    CGESVXX uses the LU factorization to compute the solution to a
!>    complex system of linear equations  A * X = B,  where A is an
!>    N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!>    If requested, both normwise and maximum componentwise error bounds
!>    are returned. CGESVXX will return a solution with a tiny
!>    guaranteed error (O(eps) where eps is the working machine
!>    precision) unless the matrix is very ill-conditioned, in which
!>    case a warning is returned. Relevant condition numbers also are
!>    calculated and returned.
!>
!>    CGESVXX accepts user-provided factorizations and equilibration
!>    factors; see the definitions of the FACT and EQUED options.
!>    Solving with refinement and using a factorization from a previous
!>    CGESVXX call will also produce a solution with either O(eps)
!>    errors or warnings, but we cannot make that claim for general
!>    user-provided factorizations and equilibration factors if they
!>    differ from what CGESVXX would itself produce.
!>
Description:
!>
!>    The following steps are performed:
!>
!>    1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>
!>      TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
!>      TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
!>      TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
!>
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
!>    or diag(C)*B (if TRANS = 'T' or 'C').
!>
!>    2. If FACT = 'N' or 'E', the LU decomposition is used to factor
!>    the matrix A (after equilibration if FACT = 'E') as
!>
!>      A = P * L * U,
!>
!>    where P is a permutation matrix, L is a unit lower triangular
!>    matrix, and U is upper triangular.
!>
!>    3. If some U(i,i)=0, so that U is exactly singular, then the
!>    routine returns with INFO = i. Otherwise, the factored form of A
!>    is used to estimate the condition number of the matrix A (see
!>    argument RCOND). If the reciprocal of the condition number is less
!>    than machine precision, the routine still goes on to solve for X
!>    and compute error bounds as described below.
!>
!>    4. The system of equations is solved for X using the factored form
!>    of A.
!>
!>    5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
!>    the routine will use iterative refinement to try to get a small
!>    error and error bounds.  Refinement calculates the residual to at
!>    least twice the working precision.
!>
!>    6. If equilibration was used, the matrix X is premultiplied by
!>    diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
!>    that it solves the original system before equilibration.
!> 
!>     Some optional parameters are bundled in the PARAMS array.  These
!>     settings determine how refinement is performed, but often the
!>     defaults are acceptable.  If the defaults are acceptable, users
!>     can pass NPARAMS = 0 which prevents the source code from accessing
!>     the PARAMS argument.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>     Specifies whether or not the factored form of the matrix A is
!>     supplied on entry, and if not, whether the matrix A should be
!>     equilibrated before it is factored.
!>       = 'F':  On entry, AF and IPIV contain the factored form of A.
!>               If EQUED is not 'N', the matrix A has been
!>               equilibrated with scaling factors given by R and C.
!>               A, AF, and IPIV are not modified.
!>       = 'N':  The matrix A will be copied to AF and factored.
!>       = 'E':  The matrix A will be equilibrated if necessary, then
!>               copied to AF and factored.
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>     Specifies the form of the system of equations:
!>       = 'N':  A * X = B     (No transpose)
!>       = 'T':  A**T * X = B  (Transpose)
!>       = 'C':  A**H * X = B  (Conjugate Transpose)
!>
[in]N
!>          N is INTEGER
!>     The number of linear equations, i.e., the order 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 N-by-N matrix A.  If FACT = 'F' and EQUED is
!>     not 'N', then A must have been equilibrated by the scaling
!>     factors in R and/or C.  A is not modified if FACT = 'F' or
!>     'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>     On exit, if EQUED .ne. 'N', A is scaled as follows:
!>     EQUED = 'R':  A := diag(R) * A
!>     EQUED = 'C':  A := A * diag(C)
!>     EQUED = 'B':  A := diag(R) * A * diag(C).
!>
[in]LDA
!>          LDA is INTEGER
!>     The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]AF
!>          AF is COMPLEX array, dimension (LDAF,N)
!>     If FACT = 'F', then AF is an input argument and on entry
!>     contains the factors L and U from the factorization
!>     A = P*L*U as computed by CGETRF.  If EQUED .ne. 'N', then
!>     AF is the factored form of the equilibrated matrix A.
!>
!>     If FACT = 'N', then AF is an output argument and on exit
!>     returns the factors L and U from the factorization A = P*L*U
!>     of the original matrix A.
!>
!>     If FACT = 'E', then AF is an output argument and on exit
!>     returns the factors L and U from the factorization A = P*L*U
!>     of the equilibrated matrix A (see the description of A for
!>     the form of the equilibrated matrix).
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>     If FACT = 'F', then IPIV is an input argument and on entry
!>     contains the pivot indices from the factorization A = P*L*U
!>     as computed by CGETRF; row i of the matrix was interchanged
!>     with row IPIV(i).
!>
!>     If FACT = 'N', then IPIV is an output argument and on exit
!>     contains the pivot indices from the factorization A = P*L*U
!>     of the original matrix A.
!>
!>     If FACT = 'E', then IPIV is an output argument and on exit
!>     contains the pivot indices from the factorization A = P*L*U
!>     of the equilibrated matrix A.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>     Specifies the form of equilibration that was done.
!>       = 'N':  No equilibration (always true if FACT = 'N').
!>       = 'R':  Row equilibration, i.e., A has been premultiplied by
!>               diag(R).
!>       = 'C':  Column equilibration, i.e., A has been postmultiplied
!>               by diag(C).
!>       = 'B':  Both row and column equilibration, i.e., A has been
!>               replaced by diag(R) * A * diag(C).
!>     EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>     output argument.
!>
[in,out]R
!>          R is REAL array, dimension (N)
!>     The row scale factors for A.  If EQUED = 'R' or 'B', A is
!>     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
!>     is not accessed.  R is an input argument if FACT = 'F';
!>     otherwise, R is an output argument.  If FACT = 'F' and
!>     EQUED = 'R' or 'B', each element of R must be positive.
!>     If R is output, each element of R is a power of the radix.
!>     If R is input, each element of R should be a power of the radix
!>     to ensure a reliable solution and error estimates. Scaling by
!>     powers of the radix does not cause rounding errors unless the
!>     result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in,out]C
!>          C is REAL array, dimension (N)
!>     The column scale factors for A.  If EQUED = 'C' or 'B', A is
!>     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
!>     is not accessed.  C is an input argument if FACT = 'F';
!>     otherwise, C is an output argument.  If FACT = 'F' and
!>     EQUED = 'C' or 'B', each element of C must be positive.
!>     If C is output, each element of C is a power of the radix.
!>     If C is input, each element of C should be a power of the radix
!>     to ensure a reliable solution and error estimates. Scaling by
!>     powers of the radix does not cause rounding errors unless the
!>     result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in,out]B
!>          B is COMPLEX array, dimension (LDB,NRHS)
!>     On entry, the N-by-NRHS right hand side matrix B.
!>     On exit,
!>     if EQUED = 'N', B is not modified;
!>     if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
!>        diag(R)*B;
!>     if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
!>        overwritten by diag(C)*B.
!>
[in]LDB
!>          LDB is INTEGER
!>     The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is COMPLEX array, dimension (LDX,NRHS)
!>     If INFO = 0, the N-by-NRHS solution matrix X to the original
!>     system of equations.  Note that A and B are modified on exit
!>     if EQUED .ne. 'N', and the solution to the equilibrated system is
!>     inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
!>     inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
!>
[in]LDX
!>          LDX is INTEGER
!>     The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[out]RPVGRW
!>          RPVGRW is REAL
!>     Reciprocal pivot growth.  On exit, this contains the reciprocal
!>     pivot growth factor norm(A)/norm(U). The 
!>     norm is used.  If this is much less than 1, then the stability of
!>     the LU factorization of the (equilibrated) matrix A could be poor.
!>     This also means that the solution X, estimated condition numbers,
!>     and error bounds could be unreliable. If factorization fails with
!>     0<INFO<=N, then this contains the reciprocal pivot growth factor
!>     for the leading INFO columns of A.  In CGESVX, this quantity is
!>     returned in WORK(1).
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>     Componentwise relative backward error.  This is the
!>     componentwise relative backward error of each solution vector X(j)
!>     (i.e., the smallest relative change in any element of A or B that
!>     makes X(j) an exact solution).
!>
[in]N_ERR_BNDS
!>          N_ERR_BNDS is INTEGER
!>     Number of error bounds to return for each right hand side
!>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
!>     ERR_BNDS_COMP below.
!>
[out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is REAL array, dimension NPARAMS
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the double-precision refinement algorithm,
!>                    possibly with doubled-single computations if the
!>                    compilation environment does not support DOUBLE
!>                    PRECISION.
!>              (other values are reserved for future use)
!>
!>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
!>            computations allowed for refinement.
!>         Default: 10
!>         Aggressive: Set to 100 to permit convergence using approximate
!>                     factorizations or factorizations other than LU. If
!>                     the factorization uses a technique other than
!>                     Gaussian elimination, the guarantees in
!>                     err_bnds_norm and err_bnds_comp may no longer be
!>                     trustworthy.
!>
!>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
!>            will attempt to find a solution with small componentwise
!>            relative error in the double-precision algorithm.  Positive
!>            is true, 0.0 is false.
!>         Default: 1.0 (attempt componentwise convergence)
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (2*N)
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit. The solution to every right-hand side is
!>         guaranteed.
!>       < 0:  If INFO = -i, the i-th argument had an illegal value
!>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
!>         has been completed, but the factor U is exactly singular, so
!>         the solution and error bounds could not be computed. RCOND = 0
!>         is returned.
!>       = N+J: The solution corresponding to the Jth right-hand side is
!>         not guaranteed. The solutions corresponding to other right-
!>         hand sides K with K > J may not be guaranteed as well, but
!>         only the first such right-hand side is reported. If a small
!>         componentwise error is not requested (PARAMS(3) = 0.0) then
!>         the Jth right-hand side is the first with a normwise error
!>         bound that is not guaranteed (the smallest J such
!>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
!>         the Jth right-hand side is the first with either a normwise or
!>         componentwise error bound that is not guaranteed (the smallest
!>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
!>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
!>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
!>         about all of the right-hand sides check ERR_BNDS_NORM or
!>         ERR_BNDS_COMP.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 538 of file cgesvxx.f.

543*
544* -- LAPACK driver routine --
545* -- LAPACK is a software package provided by Univ. of Tennessee, --
546* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
547*
548* .. Scalar Arguments ..
549 CHARACTER EQUED, FACT, TRANS
550 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
551 $ N_ERR_BNDS
552 REAL RCOND, RPVGRW
553* ..
554* .. Array Arguments ..
555 INTEGER IPIV( * )
556 COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
557 $ X( LDX , * ),WORK( * )
558 REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
559 $ ERR_BNDS_NORM( NRHS, * ),
560 $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
561* ..
562*
563* ==================================================================
564*
565* .. Parameters ..
566 REAL ZERO, ONE
567 parameter( zero = 0.0e+0, one = 1.0e+0 )
568 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
569 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
570 INTEGER CMP_ERR_I, PIV_GROWTH_I
571 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
572 $ berr_i = 3 )
573 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
574 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
575 $ piv_growth_i = 9 )
576* ..
577* .. Local Scalars ..
578 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
579 INTEGER INFEQU, J
580 REAL AMAX, BIGNUM, COLCND, RCMAX, RCMIN,
581 $ ROWCND, SMLNUM
582* ..
583* .. External Functions ..
584 EXTERNAL lsame, slamch, cla_gerpvgrw
585 LOGICAL LSAME
586 REAL SLAMCH, CLA_GERPVGRW
587* ..
588* .. External Subroutines ..
589 EXTERNAL cgeequb, cgetrf, cgetrs, clacpy, claqge,
590 $ xerbla, clascl2, cgerfsx
591* ..
592* .. Intrinsic Functions ..
593 INTRINSIC max, min
594* ..
595* .. Executable Statements ..
596*
597 info = 0
598 nofact = lsame( fact, 'N' )
599 equil = lsame( fact, 'E' )
600 notran = lsame( trans, 'N' )
601 smlnum = slamch( 'Safe minimum' )
602 bignum = one / smlnum
603 IF( nofact .OR. equil ) THEN
604 equed = 'N'
605 rowequ = .false.
606 colequ = .false.
607 ELSE
608 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
609 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
610 END IF
611*
612* Default is failure. If an input parameter is wrong or
613* factorization fails, make everything look horrible. Only the
614* pivot growth is set here, the rest is initialized in CGERFSX.
615*
616 rpvgrw = zero
617*
618* Test the input parameters. PARAMS is not tested until CGERFSX.
619*
620 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
621 $ lsame( fact, 'F' ) ) THEN
622 info = -1
623 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
624 $ lsame( trans, 'C' ) ) THEN
625 info = -2
626 ELSE IF( n.LT.0 ) THEN
627 info = -3
628 ELSE IF( nrhs.LT.0 ) THEN
629 info = -4
630 ELSE IF( lda.LT.max( 1, n ) ) THEN
631 info = -6
632 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
633 info = -8
634 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
635 $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
636 info = -10
637 ELSE
638 IF( rowequ ) THEN
639 rcmin = bignum
640 rcmax = zero
641 DO 10 j = 1, n
642 rcmin = min( rcmin, r( j ) )
643 rcmax = max( rcmax, r( j ) )
644 10 CONTINUE
645 IF( rcmin.LE.zero ) THEN
646 info = -11
647 ELSE IF( n.GT.0 ) THEN
648 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
649 ELSE
650 rowcnd = one
651 END IF
652 END IF
653 IF( colequ .AND. info.EQ.0 ) THEN
654 rcmin = bignum
655 rcmax = zero
656 DO 20 j = 1, n
657 rcmin = min( rcmin, c( j ) )
658 rcmax = max( rcmax, c( j ) )
659 20 CONTINUE
660 IF( rcmin.LE.zero ) THEN
661 info = -12
662 ELSE IF( n.GT.0 ) THEN
663 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
664 ELSE
665 colcnd = one
666 END IF
667 END IF
668 IF( info.EQ.0 ) THEN
669 IF( ldb.LT.max( 1, n ) ) THEN
670 info = -14
671 ELSE IF( ldx.LT.max( 1, n ) ) THEN
672 info = -16
673 END IF
674 END IF
675 END IF
676*
677 IF( info.NE.0 ) THEN
678 CALL xerbla( 'CGESVXX', -info )
679 RETURN
680 END IF
681*
682 IF( equil ) THEN
683*
684* Compute row and column scalings to equilibrate the matrix A.
685*
686 CALL cgeequb( n, n, a, lda, r, c, rowcnd, colcnd, amax,
687 $ infequ )
688 IF( infequ.EQ.0 ) THEN
689*
690* Equilibrate the matrix.
691*
692 CALL claqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
693 $ equed )
694 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
695 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
696 END IF
697*
698* If the scaling factors are not applied, set them to 1.0.
699*
700 IF ( .NOT.rowequ ) THEN
701 DO j = 1, n
702 r( j ) = 1.0
703 END DO
704 END IF
705 IF ( .NOT.colequ ) THEN
706 DO j = 1, n
707 c( j ) = 1.0
708 END DO
709 END IF
710 END IF
711*
712* Scale the right-hand side.
713*
714 IF( notran ) THEN
715 IF( rowequ ) CALL clascl2( n, nrhs, r, b, ldb )
716 ELSE
717 IF( colequ ) CALL clascl2( n, nrhs, c, b, ldb )
718 END IF
719*
720 IF( nofact .OR. equil ) THEN
721*
722* Compute the LU factorization of A.
723*
724 CALL clacpy( 'Full', n, n, a, lda, af, ldaf )
725 CALL cgetrf( n, n, af, ldaf, ipiv, info )
726*
727* Return if INFO is non-zero.
728*
729 IF( info.GT.0 ) THEN
730*
731* Pivot in column INFO is exactly 0
732* Compute the reciprocal pivot growth factor of the
733* leading rank-deficient INFO columns of A.
734*
735 rpvgrw = cla_gerpvgrw( n, info, a, lda, af, ldaf )
736 RETURN
737 END IF
738 END IF
739*
740* Compute the reciprocal pivot growth factor RPVGRW.
741*
742 rpvgrw = cla_gerpvgrw( n, n, a, lda, af, ldaf )
743*
744* Compute the solution matrix X.
745*
746 CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
747 CALL cgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
748*
749* Use iterative refinement to improve the computed solution and
750* compute error bounds and backward error estimates for it.
751*
752 CALL cgerfsx( trans, equed, n, nrhs, a, lda, af, ldaf,
753 $ ipiv, r, c, b, ldb, x, ldx, rcond, berr,
754 $ n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params,
755 $ work, rwork, info )
756*
757* Scale solutions.
758*
759 IF ( colequ .AND. notran ) THEN
760 CALL clascl2 ( n, nrhs, c, x, ldx )
761 ELSE IF ( rowequ .AND. .NOT.notran ) THEN
762 CALL clascl2 ( n, nrhs, r, x, ldx )
763 END IF
764*
765 RETURN
766*
767* End of CGESVXX
768*
cla_gerpvgrw
real function cla_gerpvgrw(n, ncols, a, lda, af, ldaf)
CLA_GERPVGRW multiplies a square real matrix by a complex matrix.
Definition cla_gerpvgrw.f:98
cgeequb
subroutine cgeequb(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
CGEEQUB
Definition cgeequb.f:147
cgerfsx
subroutine cgerfsx(trans, equed, n, nrhs, a, lda, af, ldaf, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, rwork, info)
CGERFSX
Definition cgerfsx.f:414
clascl2
subroutine clascl2(m, n, d, x, ldx)
CLASCL2 performs diagonal scaling on a vector.
Definition clascl2.f:91

◆ cgetsls()

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 )

CGETSLS

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**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,
!>          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.

Definition at line 160 of file cgetsls.f.

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