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


This browser is not able to show SVG: try Firefox, Chrome, Safari, or Opera instead.

## Functions

subroutine zgelsx (M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, RWORK, INFO)
ZGELSX solves overdetermined or underdetermined systems for GE matrices More...

subroutine zcgesv (N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK, SWORK, RWORK, ITER, INFO)
ZCGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement) More...

subroutine zgels (TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
ZGELS solves overdetermined or underdetermined systems for GE matrices More...

subroutine zgelsd (M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, IWORK, INFO)
ZGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices More...

subroutine zgelss (M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSS solves overdetermined or underdetermined systems for GE matrices More...

subroutine zgelsy (M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSY solves overdetermined or underdetermined systems for GE matrices More...

subroutine zgesv (N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGESV computes the solution to system of linear equations A * X = B for GE matrices (simple driver) More...

subroutine zgesvx (FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
ZGESVX computes the solution to system of linear equations A * X = B for GE matrices More...

subroutine zgesvxx (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)
ZGESVXX computes the solution to system of linear equations A * X = B for GE matrices More...

## Detailed Description

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

## Function Documentation

 subroutine zcgesv ( integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, integer, dimension( * ) IPIV, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldx, * ) X, integer LDX, complex*16, dimension( n, * ) WORK, complex, dimension( * ) SWORK, double precision, dimension( * ) RWORK, integer ITER, integer INFO )

ZCGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement)

Purpose:
``` ZCGESV 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.

ZCGESV first attempts to factorize the matrix in COMPLEX and use this
factorization within an iterative refinement procedure to produce a
solution with COMPLEX*16 normwise backward error quality (see below).
If the approach fails the method switches to a COMPLEX*16
factorization and solve.

The iterative refinement is not going to be a winning strategy if
the ratio COMPLEX performance over COMPLEX*16 performance is too
small. A reasonable strategy should take the number of right-hand
sides and the size of the matrix into account. This might be done
with a call to ILAENV in the future. Up to now, we always try
iterative refinement.

The iterative refinement process is stopped if
ITER > ITERMAX
or for all the RHS we have:
RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
where
o ITER is the number of the current iteration in the iterative
refinement process
o RNRM is the infinity-norm of the residual
o XNRM is the infinity-norm of the solution
o ANRM is the infinity-operator-norm of the matrix A
o EPS is the machine epsilon returned by DLAMCH('Epsilon')
The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
respectively.```
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*16 array, dimension (LDA,N) On entry, the N-by-N coefficient matrix A. On exit, if iterative refinement has been successfully used (INFO.EQ.0 and ITER.GE.0, see description below), then A is unchanged, if double precision factorization has been used (INFO.EQ.0 and ITER.LT.0, see description below), then the array A contains 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). Corresponds either to the single precision factorization (if INFO.EQ.0 and ITER.GE.0) or the double precision factorization (if INFO.EQ.0 and ITER.LT.0).``` [in] B ``` B is COMPLEX*16 array, dimension (LDB,NRHS) The N-by-NRHS right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [out] X ``` X is COMPLEX*16 array, dimension (LDX,NRHS) If INFO = 0, the N-by-NRHS solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] WORK ``` WORK is COMPLEX*16 array, dimension (N*NRHS) This array is used to hold the residual vectors.``` [out] SWORK ``` SWORK is COMPLEX array, dimension (N*(N+NRHS)) This array is used to use the single precision matrix and the right-hand sides or solutions in single precision.``` [out] RWORK ` RWORK is DOUBLE PRECISION array, dimension (N)` [out] ITER ``` ITER is INTEGER < 0: iterative refinement has failed, COMPLEX*16 factorization has been performed -1 : the routine fell back to full precision for implementation- or machine-specific reasons -2 : narrowing the precision induced an overflow, the routine fell back to full precision -3 : failure of CGETRF -31: stop the iterative refinement after the 30th iterations > 0: iterative refinement has been sucessfully used. Returns the number of iterations``` [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) computed in COMPLEX*16 is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.```
Date
November 2011

Definition at line 203 of file zcgesv.f.

203 *
204 * -- LAPACK driver routine (version 3.4.0) --
205 * -- LAPACK is a software package provided by Univ. of Tennessee, --
206 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
207 * November 2011
208 *
209 * .. Scalar Arguments ..
210  INTEGER info, iter, lda, ldb, ldx, n, nrhs
211 * ..
212 * .. Array Arguments ..
213  INTEGER ipiv( * )
214  DOUBLE PRECISION rwork( * )
215  COMPLEX swork( * )
216  COMPLEX*16 a( lda, * ), b( ldb, * ), work( n, * ),
217  \$ x( ldx, * )
218 * ..
219 *
220 * =====================================================================
221 *
222 * .. Parameters ..
223  LOGICAL doitref
224  parameter( doitref = .true. )
225 *
226  INTEGER itermax
227  parameter( itermax = 30 )
228 *
229  DOUBLE PRECISION bwdmax
230  parameter( bwdmax = 1.0e+00 )
231 *
232  COMPLEX*16 negone, one
233  parameter( negone = ( -1.0d+00, 0.0d+00 ),
234  \$ one = ( 1.0d+00, 0.0d+00 ) )
235 *
236 * .. Local Scalars ..
237  INTEGER i, iiter, ptsa, ptsx
238  DOUBLE PRECISION anrm, cte, eps, rnrm, xnrm
239  COMPLEX*16 zdum
240 *
241 * .. External Subroutines ..
242  EXTERNAL cgetrs, cgetrf, clag2z, xerbla, zaxpy, zgemm,
243  \$ zlacpy, zlag2c
244 * ..
245 * .. External Functions ..
246  INTEGER izamax
247  DOUBLE PRECISION dlamch, zlange
248  EXTERNAL izamax, dlamch, zlange
249 * ..
250 * .. Intrinsic Functions ..
251  INTRINSIC abs, dble, max, sqrt
252 * ..
253 * .. Statement Functions ..
254  DOUBLE PRECISION cabs1
255 * ..
256 * .. Statement Function definitions ..
257  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
258 * ..
259 * .. Executable Statements ..
260 *
261  info = 0
262  iter = 0
263 *
264 * Test the input parameters.
265 *
266  IF( n.LT.0 ) THEN
267  info = -1
268  ELSE IF( nrhs.LT.0 ) THEN
269  info = -2
270  ELSE IF( lda.LT.max( 1, n ) ) THEN
271  info = -4
272  ELSE IF( ldb.LT.max( 1, n ) ) THEN
273  info = -7
274  ELSE IF( ldx.LT.max( 1, n ) ) THEN
275  info = -9
276  END IF
277  IF( info.NE.0 ) THEN
278  CALL xerbla( 'ZCGESV', -info )
279  RETURN
280  END IF
281 *
282 * Quick return if (N.EQ.0).
283 *
284  IF( n.EQ.0 )
285  \$ RETURN
286 *
287 * Skip single precision iterative refinement if a priori slower
288 * than double precision factorization.
289 *
290  IF( .NOT.doitref ) THEN
291  iter = -1
292  GO TO 40
293  END IF
294 *
295 * Compute some constants.
296 *
297  anrm = zlange( 'I', n, n, a, lda, rwork )
298  eps = dlamch( 'Epsilon' )
299  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
300 *
301 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
302 *
303  ptsa = 1
304  ptsx = ptsa + n*n
305 *
306 * Convert B from double precision to single precision and store the
307 * result in SX.
308 *
309  CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
310 *
311  IF( info.NE.0 ) THEN
312  iter = -2
313  GO TO 40
314  END IF
315 *
316 * Convert A from double precision to single precision and store the
317 * result in SA.
318 *
319  CALL zlag2c( n, n, a, lda, swork( ptsa ), n, info )
320 *
321  IF( info.NE.0 ) THEN
322  iter = -2
323  GO TO 40
324  END IF
325 *
326 * Compute the LU factorization of SA.
327 *
328  CALL cgetrf( n, n, swork( ptsa ), n, ipiv, info )
329 *
330  IF( info.NE.0 ) THEN
331  iter = -3
332  GO TO 40
333  END IF
334 *
335 * Solve the system SA*SX = SB.
336 *
337  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
338  \$ swork( ptsx ), n, info )
339 *
340 * Convert SX back to double precision
341 *
342  CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
343 *
344 * Compute R = B - AX (R is WORK).
345 *
346  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
347 *
348  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone, a,
349  \$ lda, x, ldx, one, work, n )
350 *
351 * Check whether the NRHS normwise backward errors satisfy the
352 * stopping criterion. If yes, set ITER=0 and return.
353 *
354  DO i = 1, nrhs
355  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
356  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
357  IF( rnrm.GT.xnrm*cte )
358  \$ GO TO 10
359  END DO
360 *
361 * If we are here, the NRHS normwise backward errors satisfy the
362 * stopping criterion. We are good to exit.
363 *
364  iter = 0
365  RETURN
366 *
367  10 CONTINUE
368 *
369  DO 30 iiter = 1, itermax
370 *
371 * Convert R (in WORK) from double precision to single precision
372 * and store the result in SX.
373 *
374  CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
375 *
376  IF( info.NE.0 ) THEN
377  iter = -2
378  GO TO 40
379  END IF
380 *
381 * Solve the system SA*SX = SR.
382 *
383  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
384  \$ swork( ptsx ), n, info )
385 *
386 * Convert SX back to double precision and update the current
387 * iterate.
388 *
389  CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
390 *
391  DO i = 1, nrhs
392  CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
393  END DO
394 *
395 * Compute R = B - AX (R is WORK).
396 *
397  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
398 *
399  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone,
400  \$ a, lda, x, ldx, one, work, n )
401 *
402 * Check whether the NRHS normwise backward errors satisfy the
403 * stopping criterion. If yes, set ITER=IITER>0 and return.
404 *
405  DO i = 1, nrhs
406  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
407  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
408  IF( rnrm.GT.xnrm*cte )
409  \$ GO TO 20
410  END DO
411 *
412 * If we are here, the NRHS normwise backward errors satisfy the
413 * stopping criterion, we are good to exit.
414 *
415  iter = iiter
416 *
417  RETURN
418 *
419  20 CONTINUE
420 *
421  30 CONTINUE
422 *
423 * If we are at this place of the code, this is because we have
424 * performed ITER=ITERMAX iterations and never satisified the stopping
425 * criterion, set up the ITER flag accordingly and follow up on double
426 * precision routine.
427 *
428  iter = -itermax - 1
429 *
430  40 CONTINUE
431 *
432 * Single-precision iterative refinement failed to converge to a
433 * satisfactory solution, so we resort to double precision.
434 *
435  CALL zgetrf( n, n, a, lda, ipiv, info )
436 *
437  IF( info.NE.0 )
438  \$ RETURN
439 *
440  CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
441  CALL zgetrs( 'No transpose', n, nrhs, a, lda, ipiv, x, ldx,
442  \$ info )
443 *
444  RETURN
445 *
446 * End of ZCGESV.
447 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
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 zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine clag2z(M, N, SA, LDSA, A, LDA, INFO)
CLAG2Z converts a complex single precision matrix to a complex double precision matrix.
Definition: clag2z.f:105
subroutine zlag2c(M, N, A, LDA, SA, LDSA, INFO)
ZLAG2C converts a complex double precision matrix to a complex single precision matrix.
Definition: zlag2c.f:109
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgels ( character TRANS, integer M, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( * ) WORK, integer LWORK, integer INFO )

ZGELS solves overdetermined or underdetermined systems for GE matrices

Purpose:
``` ZGELS 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*16 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 ZGEQRF; if M < N, A is overwritten by details of its LQ factorization as returned by ZGELQF.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).``` [in,out] B ``` B is COMPLEX*16 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*16 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.```
Date
November 2011

Definition at line 184 of file zgels.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*16 a( lda, * ), b( ldb, * ), work( * )
196 * ..
197 *
198 * =====================================================================
199 *
200 * .. Parameters ..
201  DOUBLE PRECISION zero, one
202  parameter( zero = 0.0d+0, one = 1.0d+0 )
203  COMPLEX*16 czero
204  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
205 * ..
206 * .. Local Scalars ..
207  LOGICAL lquery, tpsd
208  INTEGER brow, i, iascl, ibscl, j, mn, nb, scllen, wsize
209  DOUBLE PRECISION anrm, bignum, bnrm, smlnum
210 * ..
211 * .. Local Arrays ..
212  DOUBLE PRECISION rwork( 1 )
213 * ..
214 * .. External Functions ..
215  LOGICAL lsame
216  INTEGER ilaenv
217  DOUBLE PRECISION dlamch, zlange
218  EXTERNAL lsame, ilaenv, dlamch, zlange
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL dlabad, xerbla, zgelqf, zgeqrf, zlascl, zlaset,
222  \$ ztrtrs, zunmlq, zunmqr
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC dble, max, min
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. .NOT.lquery )
247  \$ 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, 'ZGEQRF', ' ', m, n, -1, -1 )
261  IF( tpsd ) THEN
262  nb = max( nb, ilaenv( 1, 'ZUNMQR', 'LN', m, nrhs, n,
263  \$ -1 ) )
264  ELSE
265  nb = max( nb, ilaenv( 1, 'ZUNMQR', 'LC', m, nrhs, n,
266  \$ -1 ) )
267  END IF
268  ELSE
269  nb = ilaenv( 1, 'ZGELQF', ' ', m, n, -1, -1 )
270  IF( tpsd ) THEN
271  nb = max( nb, ilaenv( 1, 'ZUNMLQ', 'LC', n, nrhs, m,
272  \$ -1 ) )
273  ELSE
274  nb = max( nb, ilaenv( 1, 'ZUNMLQ', '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 ) = dble( wsize )
281 *
282  END IF
283 *
284  IF( info.NE.0 ) THEN
285  CALL xerbla( 'ZGELS ', -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 zlaset( 'Full', max( m, n ), nrhs, czero, czero, b, ldb )
295  RETURN
296  END IF
297 *
298 * Get machine parameters
299 *
300  smlnum = dlamch( 'S' ) / dlamch( 'P' )
301  bignum = one / smlnum
302  CALL dlabad( smlnum, bignum )
303 *
304 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
305 *
306  anrm = zlange( '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 zlascl( '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 zlascl( '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 zlaset( '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 = zlange( '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 zlascl( '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 zlascl( '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 zgeqrf( 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 zunmqr( '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 ztrtrs( '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 ztrtrs( '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 zunmqr( '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 zgelqf( 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 ztrtrs( '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 zunmlq( '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 zunmlq( '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 ztrtrs( '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 zlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
485  \$ info )
486  ELSE IF( iascl.EQ.2 ) THEN
487  CALL zlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
488  \$ info )
489  END IF
490  IF( ibscl.EQ.1 ) THEN
491  CALL zlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
492  \$ info )
493  ELSE IF( ibscl.EQ.2 ) THEN
494  CALL zlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
495  \$ info )
496  END IF
497 *
498  50 CONTINUE
499  work( 1 ) = dble( wsize )
500 *
501  RETURN
502 *
503 * End of ZGELS
504 *
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
Definition: zunmlq.f:169
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:137
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine ztrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
ZTRTRS
Definition: ztrtrs.f:142
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgelsd ( integer M, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, double precision, dimension( * ) S, double precision RCOND, integer RANK, complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension( * ) RWORK, integer, dimension( * ) IWORK, integer INFO )

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgelss ( integer M, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, double precision, dimension( * ) S, double precision RCOND, integer RANK, complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension( * ) RWORK, integer INFO )

ZGELSS solves overdetermined or underdetermined systems for GE matrices

Purpose:
``` ZGELSS 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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.```
Date
November 2011

Definition at line 180 of file zgelss.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  DOUBLE PRECISION rcond
189 * ..
190 * .. Array Arguments ..
191  DOUBLE PRECISION rwork( * ), s( * )
192  COMPLEX*16 a( lda, * ), b( ldb, * ), work( * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  DOUBLE PRECISION zero, one
199  parameter( zero = 0.0d+0, one = 1.0d+0 )
200  COMPLEX*16 czero, cone
201  parameter( czero = ( 0.0d+0, 0.0d+0 ),
202  \$ cone = ( 1.0d+0, 0.0d+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_zgeqrf, lwork_zunmqr, lwork_zgebrd,
210  \$ lwork_zunmbr, lwork_zungbr, lwork_zunmlq,
211  \$ lwork_zgelqf
212  DOUBLE PRECISION anrm, bignum, bnrm, eps, sfmin, smlnum, thr
213 * ..
214 * .. Local Arrays ..
215  COMPLEX*16 dum( 1 )
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL dlabad, dlascl, dlaset, xerbla, zbdsqr, zcopy,
221  \$ zunmqr
222 * ..
223 * .. External Functions ..
224  INTEGER ilaenv
225  DOUBLE PRECISION dlamch, zlange
226  EXTERNAL ilaenv, dlamch, zlange
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, 'ZGELSS', ' ', 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 ZGEQRF
271  CALL zgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
272  lwork_zgeqrf=dum(1)
273 * Compute space needed for ZUNMQR
274  CALL zunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
275  \$ ldb, dum(1), -1, info )
276  lwork_zunmqr=dum(1)
277  mm = n
278  maxwrk = max( maxwrk, n + n*ilaenv( 1, 'ZGEQRF', ' ', m,
279  \$ n, -1, -1 ) )
280  maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'ZUNMQR', '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 ZGEBRD
288  CALL zgebrd( mm, n, a, lda, s, dum(1), dum(1),
289  \$ dum(1), dum(1), -1, info )
290  lwork_zgebrd=dum(1)
291 * Compute space needed for ZUNMBR
292  CALL zunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
293  \$ b, ldb, dum(1), -1, info )
294  lwork_zunmbr=dum(1)
295 * Compute space needed for ZUNGBR
296  CALL zungbr( 'P', n, n, n, a, lda, dum(1),
297  \$ dum(1), -1, info )
298  lwork_zungbr=dum(1)
299 * Compute total workspace needed
300  maxwrk = max( maxwrk, 2*n + lwork_zgebrd )
301  maxwrk = max( maxwrk, 2*n + lwork_zunmbr )
302  maxwrk = max( maxwrk, 2*n + lwork_zungbr )
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 ZGELQF
314  CALL zgelqf( m, n, a, lda, dum(1), dum(1),
315  \$ -1, info )
316  lwork_zgelqf=dum(1)
317 * Compute space needed for ZGEBRD
318  CALL zgebrd( m, m, a, lda, s, dum(1), dum(1),
319  \$ dum(1), dum(1), -1, info )
320  lwork_zgebrd=dum(1)
321 * Compute space needed for ZUNMBR
322  CALL zunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
323  \$ dum(1), b, ldb, dum(1), -1, info )
324  lwork_zunmbr=dum(1)
325 * Compute space needed for ZUNGBR
326  CALL zungbr( 'P', m, m, m, a, lda, dum(1),
327  \$ dum(1), -1, info )
328  lwork_zungbr=dum(1)
329 * Compute space needed for ZUNMLQ
330  CALL zunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
331  \$ b, ldb, dum(1), -1, info )
332  lwork_zunmlq=dum(1)
333 * Compute total workspace needed
334  maxwrk = m + lwork_zgelqf
335  maxwrk = max( maxwrk, 3*m + m*m + lwork_zgebrd )
336  maxwrk = max( maxwrk, 3*m + m*m + lwork_zunmbr )
337  maxwrk = max( maxwrk, 3*m + m*m + lwork_zungbr )
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_zunmlq )
344  ELSE
345 *
346 * Path 2 - underdetermined
347 *
348 * Compute space needed for ZGEBRD
349  CALL zgebrd( m, n, a, lda, s, dum(1), dum(1),
350  \$ dum(1), dum(1), -1, info )
351  lwork_zgebrd=dum(1)
352 * Compute space needed for ZUNMBR
353  CALL zunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
354  \$ dum(1), b, ldb, dum(1), -1, info )
355  lwork_zunmbr=dum(1)
356 * Compute space needed for ZUNGBR
357  CALL zungbr( 'P', m, n, m, a, lda, dum(1),
358  \$ dum(1), -1, info )
359  lwork_zungbr=dum(1)
360  maxwrk = 2*m + lwork_zgebrd
361  maxwrk = max( maxwrk, 2*m + lwork_zunmbr )
362  maxwrk = max( maxwrk, 2*m + lwork_zungbr )
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( 'ZGELSS', -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 = dlamch( 'P' )
391  sfmin = dlamch( 'S' )
392  smlnum = sfmin / eps
393  bignum = one / smlnum
394  CALL dlabad( smlnum, bignum )
395 *
396 * Scale A if max element outside range [SMLNUM,BIGNUM]
397 *
398  anrm = zlange( '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 zlascl( '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 zlascl( '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 zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
417  CALL dlaset( '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 = zlange( '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 zlascl( '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 zlascl( '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 zgeqrf( 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 zunmqr( '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 zlaset( '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 zgebrd( 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 zunmbr( '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 zungbr( '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 zbdsqr( '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 zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
524  rank = rank + 1
525  ELSE
526  CALL zlaset( '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 zgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
536  \$ czero, work, ldb )
537  CALL zlacpy( '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 zgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
543  \$ ldb, czero, work, n )
544  CALL zlacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
545  20 CONTINUE
546  ELSE
547  CALL zgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
548  CALL zcopy( 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 zgelqf( 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 zlacpy( 'L', m, m, a, lda, work( il ), ldwork )
576  CALL zlaset( '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 zgebrd( 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 zunmbr( '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 zungbr( '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 zbdsqr( '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 zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
627  rank = rank + 1
628  ELSE
629  CALL zlaset( '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 zgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
640  \$ b, ldb, czero, work( iwork ), ldb )
641  CALL zlacpy( '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 zgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
647  \$ b( 1, i ), ldb, czero, work( iwork ), m )
648  CALL zlacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
649  \$ ldb )
650  40 CONTINUE
651  ELSE
652  CALL zgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
653  \$ 1, czero, work( iwork ), 1 )
654  CALL zcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
655  END IF
656 *
657 * Zero out below first M rows of B
658 *
659  CALL zlaset( '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 zunmlq( '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 zgebrd( 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 zunmbr( '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 zungbr( '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 zbdsqr( '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 zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
721  rank = rank + 1
722  ELSE
723  CALL zlaset( '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 zgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
733  \$ czero, work, ldb )
734  CALL zlacpy( '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 zgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
740  \$ ldb, czero, work, n )
741  CALL zlacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
742  60 CONTINUE
743  ELSE
744  CALL zgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
745  CALL zcopy( n, work, 1, b, 1 )
746  END IF
747  END IF
748 *
749 * Undo scaling
750 *
751  IF( iascl.EQ.1 ) THEN
752  CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
753  CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
754  \$ info )
755  ELSE IF( iascl.EQ.2 ) THEN
756  CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
757  CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
758  \$ info )
759  END IF
760  IF( ibscl.EQ.1 ) THEN
761  CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
762  ELSE IF( ibscl.EQ.2 ) THEN
763  CALL zlascl( '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 ZGELSS
770 *
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
Definition: zbdsqr.f:224
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
Definition: zunmlq.f:169
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:207
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
Definition: zungbr.f:159
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:198
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:137
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:86
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:141
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgelsx ( integer M, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, integer, dimension( * ) JPVT, double precision RCOND, integer RANK, complex*16, dimension( * ) WORK, double precision, dimension( * ) RWORK, integer INFO )

ZGELSX solves overdetermined or underdetermined systems for GE matrices

Purpose:
``` This routine is deprecated and has been replaced by routine ZGELSY.

ZGELSX 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*16 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*16 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 DOUBLE PRECISION 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*16 array, dimension (min(M,N) + max( N, 2*min(M,N)+NRHS )),``` [out] RWORK ` RWORK is DOUBLE PRECISION array, dimension (2*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Date
November 2011

Definition at line 186 of file zgelsx.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  DOUBLE PRECISION rcond
195 * ..
196 * .. Array Arguments ..
197  INTEGER jpvt( * )
198  DOUBLE PRECISION rwork( * )
199  COMPLEX*16 a( lda, * ), b( ldb, * ), work( * )
200 * ..
201 *
202 * =====================================================================
203 *
204 * .. Parameters ..
205  INTEGER imax, imin
206  parameter( imax = 1, imin = 2 )
207  DOUBLE PRECISION zero, one, done, ntdone
208  parameter( zero = 0.0d+0, one = 1.0d+0, done = zero,
209  \$ ntdone = one )
210  COMPLEX*16 czero, cone
211  parameter( czero = ( 0.0d+0, 0.0d+0 ),
212  \$ cone = ( 1.0d+0, 0.0d+0 ) )
213 * ..
214 * .. Local Scalars ..
215  INTEGER i, iascl, ibscl, ismax, ismin, j, k, mn
216  DOUBLE PRECISION anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
217  \$ smlnum
218  COMPLEX*16 c1, c2, s1, s2, t1, t2
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL xerbla, zgeqpf, zlaic1, zlascl, zlaset, zlatzm,
222  \$ ztrsm, ztzrqf, zunm2r
223 * ..
224 * .. External Functions ..
225  DOUBLE PRECISION dlamch, zlange
226  EXTERNAL dlamch, zlange
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, dconjg, 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( 'ZGELSX', -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 = dlamch( 'S' ) / dlamch( 'P' )
267  bignum = one / smlnum
268  CALL dlabad( smlnum, bignum )
269 *
270 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
271 *
272  anrm = zlange( '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 zlascl( '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 zlascl( '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 zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
291  rank = 0
292  GO TO 100
293  END IF
294 *
295  bnrm = zlange( '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 zlascl( '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 zlascl( '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 zgeqpf( 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 zlaset( '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 zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
338  \$ a( i, i ), sminpr, s1, c1 )
339  CALL zlaic1( 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 ztzrqf( 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 zunm2r( '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 ztrsm( '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 zlatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
390  \$ dconjg( 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 zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
428  CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
429  \$ info )
430  ELSE IF( iascl.EQ.2 ) THEN
431  CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
432  CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
433  \$ info )
434  END IF
435  IF( ibscl.EQ.1 ) THEN
436  CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
437  ELSE IF( ibscl.EQ.2 ) THEN
438  CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
439  END IF
440 *
441  100 CONTINUE
442 *
443  RETURN
444 *
445 * End of ZGELSX
446 *
subroutine zlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
ZLAIC1 applies one step of incremental condition estimation.
Definition: zlaic1.f:137
subroutine ztzrqf(M, N, A, LDA, TAU, INFO)
ZTZRQF
Definition: ztzrqf.f:140
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
ZUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: zunm2r.f:161
subroutine zlatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
ZLATZM
Definition: zlatzm.f:154
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine zgeqpf(M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
ZGEQPF
Definition: zgeqpf.f:150

Here is the call graph for this function:

 subroutine zgelsy ( integer M, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, integer, dimension( * ) JPVT, double precision RCOND, integer RANK, complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension( * ) RWORK, integer INFO )

ZGELSY solves overdetermined or underdetermined systems for GE matrices

Purpose:
``` ZGELSY 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*16 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*16 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 DOUBLE PRECISION 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*16 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 ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR, and ZUNMRZ. 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 DOUBLE PRECISION array, dimension (2*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
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 zgelsy.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  DOUBLE PRECISION rcond
221 * ..
222 * .. Array Arguments ..
223  INTEGER jpvt( * )
224  DOUBLE PRECISION rwork( * )
225  COMPLEX*16 a( lda, * ), b( ldb, * ), work( * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  INTEGER imax, imin
232  parameter( imax = 1, imin = 2 )
233  DOUBLE PRECISION zero, one
234  parameter( zero = 0.0d+0, one = 1.0d+0 )
235  COMPLEX*16 czero, cone
236  parameter( czero = ( 0.0d+0, 0.0d+0 ),
237  \$ cone = ( 1.0d+0, 0.0d+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  DOUBLE PRECISION anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
244  \$ smlnum, wsize
245  COMPLEX*16 c1, c2, s1, s2
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL dlabad, xerbla, zcopy, zgeqp3, zlaic1, zlascl,
250 * ..
251 * .. External Functions ..
252  INTEGER ilaenv
253  DOUBLE PRECISION dlamch, zlange
254  EXTERNAL ilaenv, dlamch, zlange
255 * ..
256 * .. Intrinsic Functions ..
257  INTRINSIC abs, dble, dcmplx, max, min
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, 'ZGEQRF', ' ', m, n, -1, -1 )
269  nb2 = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
270  nb3 = ilaenv( 1, 'ZUNMQR', ' ', m, n, nrhs, -1 )
271  nb4 = ilaenv( 1, 'ZUNMRQ', ' ', 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 ) = dcmplx( 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. .NOT.
287  \$ lquery ) THEN
288  info = -12
289  END IF
290 *
291  IF( info.NE.0 ) THEN
292  CALL xerbla( 'ZGELSY', -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 = dlamch( 'S' ) / dlamch( 'P' )
308  bignum = one / smlnum
309  CALL dlabad( smlnum, bignum )
310 *
311 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
312 *
313  anrm = zlange( '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 zlascl( '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 zlascl( '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 zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
332  rank = 0
333  GO TO 70
334  END IF
335 *
336  bnrm = zlange( '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 zlascl( '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 zlascl( '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 zgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356  \$ lwork-mn, rwork, info )
357  wsize = mn + dble( 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 zlaset( '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 zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
380  \$ a( i, i ), sminpr, s1, c1 )
381  CALL zlaic1( 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 ztzrzf( 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 zunmqr( '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+dble( 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 ztrsm( '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 zunmrz( '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 zcopy( 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 zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
457  CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
458  \$ info )
459  ELSE IF( iascl.EQ.2 ) THEN
460  CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
461  CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
462  \$ info )
463  END IF
464  IF( ibscl.EQ.1 ) THEN
465  CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
466  ELSE IF( ibscl.EQ.2 ) THEN
467  CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
468  END IF
469 *
470  70 CONTINUE
471  work( 1 ) = dcmplx( lwkopt )
472 *
473  RETURN
474 *
475 * End of ZGELSY
476 *
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
ZLAIC1 applies one step of incremental condition estimation.
Definition: zlaic1.f:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
ZGEQP3
Definition: zgeqp3.f:161
subroutine zunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRZ
Definition: zunmrz.f:189
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine ztzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZTZRZF
Definition: ztzrzf.f:153
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgesv ( integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, integer, dimension( * ) IPIV, complex*16, dimension( ldb, * ) B, integer LDB, integer INFO )

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

Purpose:
``` ZGESV 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*16 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*16 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.```
Date
November 2011

Definition at line 124 of file zgesv.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*16 a( lda, * ), b( ldb, * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. External Subroutines ..
141  EXTERNAL xerbla, zgetrf, zgetrs
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( 'ZGESV ', -info )
162  RETURN
163  END IF
164 *
165 * Compute the LU factorization of A.
166 *
167  CALL zgetrf( 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 zgetrs( 'No transpose', n, nrhs, a, lda, ipiv, b, ldb,
173  \$ info )
174  END IF
175  RETURN
176 *
177 * End of ZGESV
178 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgesvx ( character FACT, character TRANS, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldaf, * ) AF, integer LDAF, integer, dimension( * ) IPIV, character EQUED, double precision, dimension( * ) R, double precision, dimension( * ) C, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldx, * ) X, integer LDX, double precision RCOND, double precision, dimension( * ) FERR, double precision, dimension( * ) BERR, complex*16, dimension( * ) WORK, double precision, dimension( * ) RWORK, integer INFO )

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

Purpose:
``` ZGESVX 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*16 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*16 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 ZGETRF. 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 ZGETRF; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (2*N)` [out] RWORK ``` RWORK is DOUBLE PRECISION 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 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.```
Date
April 2012

Definition at line 352 of file zgesvx.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  DOUBLE PRECISION rcond
362 * ..
363 * .. Array Arguments ..
364  INTEGER ipiv( * )
365  DOUBLE PRECISION berr( * ), c( * ), ferr( * ), r( * ),
366  \$ rwork( * )
367  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
368  \$ work( * ), x( ldx, * )
369 * ..
370 *
371 * =====================================================================
372 *
373 * .. Parameters ..
374  DOUBLE PRECISION zero, one
375  parameter( zero = 0.0d+0, one = 1.0d+0 )
376 * ..
377 * .. Local Scalars ..
378  LOGICAL colequ, equil, nofact, notran, rowequ
379  CHARACTER norm
380  INTEGER i, infequ, j
381  DOUBLE PRECISION amax, anorm, bignum, colcnd, rcmax, rcmin,
382  \$ rowcnd, rpvgrw, smlnum
383 * ..
384 * .. External Functions ..
385  LOGICAL lsame
386  DOUBLE PRECISION dlamch, zlange, zlantr
387  EXTERNAL lsame, dlamch, zlange, zlantr
388 * ..
389 * .. External Subroutines ..
390  EXTERNAL xerbla, zgecon, zgeequ, zgerfs, zgetrf, zgetrs,
391  \$ zlacpy, zlaqge
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 = dlamch( '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( 'ZGESVX', -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 zgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
482  IF( infequ.EQ.0 ) THEN
483 *
484 * Equilibrate the matrix.
485 *
486  CALL zlaqge( 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 zlacpy( 'Full', n, n, a, lda, af, ldaf )
516  CALL zgetrf( 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 = zlantr( 'M', 'U', 'N', info, info, af, ldaf,
526  \$ rwork )
527  IF( rpvgrw.EQ.zero ) THEN
528  rpvgrw = one
529  ELSE
530  rpvgrw = zlange( '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 = zlange( norm, n, n, a, lda, rwork )
548  rpvgrw = zlantr( 'M', 'U', 'N', n, n, af, ldaf, rwork )
549  IF( rpvgrw.EQ.zero ) THEN
550  rpvgrw = one
551  ELSE
552  rpvgrw = zlange( 'M', n, n, a, lda, rwork ) / rpvgrw
553  END IF
554 *
555 * Compute the reciprocal of the condition number of A.
556 *
557  CALL zgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, info )
558 *
559 * Compute the solution matrix X.
560 *
561  CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
562  CALL zgetrs( 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 zgerfs( 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.dlamch( 'Epsilon' ) )
598  \$ info = n + 1
599 *
600  rwork( 1 ) = rpvgrw
601  RETURN
602 *
603 * End of ZGESVX
604 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function zlantr(NORM, UPLO, DIAG, M, N, A, LDA, WORK)
ZLANTR 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: zlantr.f:144
subroutine zgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZGECON
Definition: zgecon.f:126
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zlaqge(M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, EQUED)
ZLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ...
Definition: zlaqge.f:145
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine zgeequ(M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
ZGEEQU
Definition: zgeequ.f:142
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
subroutine zgerfs(TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZGERFS
Definition: zgerfs.f:188
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123

Here is the call graph for this function:

Here is the caller graph for this function:

 subroutine zgesvxx ( character FACT, character TRANS, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldaf, * ) AF, integer LDAF, integer, dimension( * ) IPIV, character EQUED, double precision, dimension( * ) R, double precision, dimension( * ) C, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldx , * ) X, integer LDX, double precision RCOND, double precision RPVGRW, double precision, dimension( * ) BERR, integer N_ERR_BNDS, double precision, dimension( nrhs, * ) ERR_BNDS_NORM, double precision, dimension( nrhs, * ) ERR_BNDS_COMP, integer NPARAMS, double precision, dimension( * ) PARAMS, complex*16, dimension( * ) WORK, double precision, dimension( * ) RWORK, integer INFO )

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

Purpose:
```    ZGESVXX uses the LU factorization to compute the solution to a
complex*16 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. ZGESVXX 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.

ZGESVXX 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
ZGESVXX 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 ZGESVXX would itself produce.```
Description:
```    The following steps are performed:

1. If FACT = 'E', double precision 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*16 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*16 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 ZGETRF. 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 ZGETRF; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 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.```
Date
April 2012

Definition at line 542 of file zgesvxx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: