LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for complex:

Functions

subroutine cgels (TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
  CGELS solves overdetermined or underdetermined systems for GE matrices More...
 
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 More...
 
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 More...
 
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 More...
 
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) More...
 
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 More...
 
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 More...
 
subroutine cgelsx (M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, RWORK, INFO)
  CGELSX solves overdetermined or underdetermined systems for GE matrices More...
 

Detailed Description

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

Function Documentation

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 undetermined 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.
Date
November 2011

Definition at line 184 of file cgels.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 tranformations, reducing the original problem
     into a "bidiagonal least squares problem" (BLS)
 (2) Solve the BLS using a divide and conquer approach.
 (3) Apply back all the Householder tranformations 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.
Date
November 2011
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 227 of file cgelsd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 180 of file cgelss.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 186 of file cgelsx.f.

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

Here is the call graph for this function:

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.
Date
November 2011
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 212 of file cgelsy.f.

212 *
213 * -- LAPACK driver routine (version 3.4.0) --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 * November 2011
217 *
218 * .. Scalar Arguments ..
219  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
220  REAL rcond
221 * ..
222 * .. Array Arguments ..
223  INTEGER jpvt( * )
224  REAL rwork( * )
225  COMPLEX a( lda, * ), b( ldb, * ), work( * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  INTEGER imax, imin
232  parameter( imax = 1, imin = 2 )
233  REAL zero, one
234  parameter( zero = 0.0e+0, one = 1.0e+0 )
235  COMPLEX czero, cone
236  parameter( czero = ( 0.0e+0, 0.0e+0 ),
237  $ cone = ( 1.0e+0, 0.0e+0 ) )
238 * ..
239 * .. Local Scalars ..
240  LOGICAL lquery
241  INTEGER i, iascl, ibscl, ismax, ismin, j, lwkopt, mn,
242  $ nb, nb1, nb2, nb3, nb4
243  REAL anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
244  $ smlnum, wsize
245  COMPLEX c1, c2, s1, s2
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL ccopy, cgeqp3, claic1, clascl, claset, ctrsm,
250 * ..
251 * .. External Functions ..
252  INTEGER ilaenv
253  REAL clange, slamch
254  EXTERNAL clange, ilaenv, slamch
255 * ..
256 * .. Intrinsic Functions ..
257  INTRINSIC abs, max, min, REAL, cmplx
258 * ..
259 * .. Executable Statements ..
260 *
261  mn = min( m, n )
262  ismin = mn + 1
263  ismax = 2*mn + 1
264 *
265 * Test the input arguments.
266 *
267  info = 0
268  nb1 = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
269  nb2 = ilaenv( 1, 'CGERQF', ' ', m, n, -1, -1 )
270  nb3 = ilaenv( 1, 'CUNMQR', ' ', m, n, nrhs, -1 )
271  nb4 = ilaenv( 1, 'CUNMRQ', ' ', m, n, nrhs, -1 )
272  nb = max( nb1, nb2, nb3, nb4 )
273  lwkopt = max( 1, mn+2*n+nb*(n+1), 2*mn+nb*nrhs )
274  work( 1 ) = cmplx( lwkopt )
275  lquery = ( lwork.EQ.-1 )
276  IF( m.LT.0 ) THEN
277  info = -1
278  ELSE IF( n.LT.0 ) THEN
279  info = -2
280  ELSE IF( nrhs.LT.0 ) THEN
281  info = -3
282  ELSE IF( lda.LT.max( 1, m ) ) THEN
283  info = -5
284  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
285  info = -7
286  ELSE IF( lwork.LT.( mn+max( 2*mn, n+1, mn+nrhs ) ) .AND.
287  $ .NOT.lquery ) THEN
288  info = -12
289  END IF
290 *
291  IF( info.NE.0 ) THEN
292  CALL xerbla( 'CGELSY', -info )
293  RETURN
294  ELSE IF( lquery ) THEN
295  RETURN
296  END IF
297 *
298 * Quick return if possible
299 *
300  IF( min( m, n, nrhs ).EQ.0 ) THEN
301  rank = 0
302  RETURN
303  END IF
304 *
305 * Get machine parameters
306 *
307  smlnum = slamch( 'S' ) / slamch( 'P' )
308  bignum = one / smlnum
309  CALL slabad( smlnum, bignum )
310 *
311 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
312 *
313  anrm = clange( 'M', m, n, a, lda, rwork )
314  iascl = 0
315  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
316 *
317 * Scale matrix norm up to SMLNUM
318 *
319  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
320  iascl = 1
321  ELSE IF( anrm.GT.bignum ) THEN
322 *
323 * Scale matrix norm down to BIGNUM
324 *
325  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
326  iascl = 2
327  ELSE IF( anrm.EQ.zero ) THEN
328 *
329 * Matrix all zero. Return zero solution.
330 *
331  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
332  rank = 0
333  GO TO 70
334  END IF
335 *
336  bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
337  ibscl = 0
338  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
339 *
340 * Scale matrix norm up to SMLNUM
341 *
342  CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
343  ibscl = 1
344  ELSE IF( bnrm.GT.bignum ) THEN
345 *
346 * Scale matrix norm down to BIGNUM
347 *
348  CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
349  ibscl = 2
350  END IF
351 *
352 * Compute QR factorization with column pivoting of A:
353 * A * P = Q * R
354 *
355  CALL cgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356  $ lwork-mn, rwork, info )
357  wsize = mn + REAL( WORK( MN+1 ) )
358 *
359 * complex workspace: MN+NB*(N+1). real workspace 2*N.
360 * Details of Householder rotations stored in WORK(1:MN).
361 *
362 * Determine RANK using incremental condition estimation
363 *
364  work( ismin ) = cone
365  work( ismax ) = cone
366  smax = abs( a( 1, 1 ) )
367  smin = smax
368  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
369  rank = 0
370  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
371  GO TO 70
372  ELSE
373  rank = 1
374  END IF
375 *
376  10 CONTINUE
377  IF( rank.LT.mn ) THEN
378  i = rank + 1
379  CALL claic1( imin, rank, work( ismin ), smin, a( 1, i ),
380  $ a( i, i ), sminpr, s1, c1 )
381  CALL claic1( imax, rank, work( ismax ), smax, a( 1, i ),
382  $ a( i, i ), smaxpr, s2, c2 )
383 *
384  IF( smaxpr*rcond.LE.sminpr ) THEN
385  DO 20 i = 1, rank
386  work( ismin+i-1 ) = s1*work( ismin+i-1 )
387  work( ismax+i-1 ) = s2*work( ismax+i-1 )
388  20 CONTINUE
389  work( ismin+rank ) = c1
390  work( ismax+rank ) = c2
391  smin = sminpr
392  smax = smaxpr
393  rank = rank + 1
394  GO TO 10
395  END IF
396  END IF
397 *
398 * complex workspace: 3*MN.
399 *
400 * Logically partition R = [ R11 R12 ]
401 * [ 0 R22 ]
402 * where R11 = R(1:RANK,1:RANK)
403 *
404 * [R11,R12] = [ T11, 0 ] * Y
405 *
406  IF( rank.LT.n )
407  $ CALL ctzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
408  $ lwork-2*mn, info )
409 *
410 * complex workspace: 2*MN.
411 * Details of Householder rotations stored in WORK(MN+1:2*MN)
412 *
413 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
414 *
415  CALL cunmqr( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
416  $ work( 1 ), b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
417  wsize = max( wsize, 2*mn+REAL( WORK( 2*MN+1 ) ) )
418 *
419 * complex workspace: 2*MN+NB*NRHS.
420 *
421 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
422 *
423  CALL ctrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
424  $ nrhs, cone, a, lda, b, ldb )
425 *
426  DO 40 j = 1, nrhs
427  DO 30 i = rank + 1, n
428  b( i, j ) = czero
429  30 CONTINUE
430  40 CONTINUE
431 *
432 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
433 *
434  IF( rank.LT.n ) THEN
435  CALL cunmrz( 'Left', 'Conjugate transpose', n, nrhs, rank,
436  $ n-rank, a, lda, work( mn+1 ), b, ldb,
437  $ work( 2*mn+1 ), lwork-2*mn, info )
438  END IF
439 *
440 * complex workspace: 2*MN+NRHS.
441 *
442 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
443 *
444  DO 60 j = 1, nrhs
445  DO 50 i = 1, n
446  work( jpvt( i ) ) = b( i, j )
447  50 CONTINUE
448  CALL ccopy( n, work( 1 ), 1, b( 1, j ), 1 )
449  60 CONTINUE
450 *
451 * complex workspace: N.
452 *
453 * Undo scaling
454 *
455  IF( iascl.EQ.1 ) THEN
456  CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
457  CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
458  $ info )
459  ELSE IF( iascl.EQ.2 ) THEN
460  CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
461  CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
462  $ info )
463  END IF
464  IF( ibscl.EQ.1 ) THEN
465  CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
466  ELSE IF( ibscl.EQ.2 ) THEN
467  CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
468  END IF
469 *
470  70 CONTINUE
471  work( 1 ) = cmplx( lwkopt )
472 *
473  RETURN
474 *
475 * End of CGELSY
476 *
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:161
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine claic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
CLAIC1 applies one step of incremental condition estimation.
Definition: claic1.f:137
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
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:141
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
subroutine cunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMRZ
Definition: cunmrz.f:189
subroutine ctzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CTZRZF
Definition: ctzrzf.f:153

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 124 of file cgesv.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "max absolute element" 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.
Date
April 2012

Definition at line 352 of file cgesvx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "max absolute element"
     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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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 .LT. 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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 .LE. 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 .LT. 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.
Date
April 2012

Definition at line 545 of file cgesvxx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: