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

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)

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 227 of file 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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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

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

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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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

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

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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Contributors:
A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

Definition at line 212 of file 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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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)

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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

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

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<INFO<=N, then
          RWORK(1) contains the reciprocal pivot growth factor for the
          leading INFO columns of A.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, and i is
                <= N:  U(i,i) is exactly zero.  The factorization has
                       been completed, but the factor U is exactly
                       singular, so the solution and error bounds
                       could not be computed. RCOND = 0 is returned.
                = N+1: U is nonsingular, but RCOND is less than machine
                       precision, meaning that the matrix is singular
                       to working precision.  Nevertheless, the
                       solution and error bounds are computed because
                       there are a number of situations where the
                       computed solution can be more accurate than the
                       value of RCOND would suggest.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 352 of file 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

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

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<INFO<=N, then this contains the reciprocal pivot growth factor
     for the leading INFO columns of A.  In ZGESVX, this quantity is
     returned in WORK(1).
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
     Componentwise relative backward error.  This is the
     componentwise relative backward error of each solution vector X(j)
     (i.e., the smallest relative change in any element of A or B that
     makes X(j) an exact solution).
[in]N_ERR_BNDS
          N_ERR_BNDS is INTEGER
     Number of error bounds to return for each right hand side
     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
     ERR_BNDS_COMP below.
[out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
     the first (:,N_ERR_BNDS) entries are returned.

     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_COMP(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*(A*diag(x)), where x is the solution for the
              current right-hand side and S scales each row of
              A*diag(x) by a power of the radix so all absolute row
              sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[in]NPARAMS
          NPARAMS is INTEGER
     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
     PARAMS array is never referenced and default values are used.
[in,out]PARAMS
          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
     that entry will be filled with default value used for that
     parameter.  Only positions up to NPARAMS are accessed; defaults
     are used for higher-numbered parameters.

       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
            refinement or not.
         Default: 1.0D+0
            = 0.0 : No refinement is performed, and no error bounds are
                    computed.
            = 1.0 : Use the extra-precise refinement algorithm.
              (other values are reserved for future use)

       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
            computations allowed for refinement.
         Default: 10
         Aggressive: Set to 100 to permit convergence using approximate
                     factorizations or factorizations other than LU. If
                     the factorization uses a technique other than
                     Gaussian elimination, the guarantees in
                     err_bnds_norm and err_bnds_comp may no longer be
                     trustworthy.

       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
            will attempt to find a solution with small componentwise
            relative error in the double-precision algorithm.  Positive
            is true, 0.0 is false.
         Default: 1.0 (attempt componentwise convergence)
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit. The solution to every right-hand side is
         guaranteed.
       < 0:  If INFO = -i, the i-th argument had an illegal value
       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
         has been completed, but the factor U is exactly singular, so
         the solution and error bounds could not be computed. RCOND = 0
         is returned.
       = N+J: The solution corresponding to the Jth right-hand side is
         not guaranteed. The solutions corresponding to other right-
         hand sides K with K > J may not be guaranteed as well, but
         only the first such right-hand side is reported. If a small
         componentwise error is not requested (PARAMS(3) = 0.0) then
         the Jth right-hand side is the first with a normwise error
         bound that is not guaranteed (the smallest J such
         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
         the Jth right-hand side is the first with either a normwise or
         componentwise error bound that is not guaranteed (the smallest
         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
         about all of the right-hand sides check ERR_BNDS_NORM or
         ERR_BNDS_COMP.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 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: