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

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
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)
subroutine zgels (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
  ZGELS solves overdetermined or underdetermined systems for GE matrices
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
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
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
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)
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
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
subroutine zgetsls (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
 ZGETSLS

Detailed Description

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

Function Documentation

◆ zcgesv()

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 = 0 and ITER >= 0, see description below), then A is
!>          unchanged, if double precision factorization has been used
!>          (INFO = 0 and ITER < 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 = 0 and ITER >= 0) or the double precision
!>          factorization (if INFO = 0 and ITER < 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 successfully 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.

Definition at line 199 of file zcgesv.f.

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

◆ zgels()

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 underdetermined system A**H * X = B.
!>
!> 4. If TRANS = 'C' and m < n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A**H * X ||.
!>
!> Several right hand side vectors b and solution vectors x can be
!> handled in a single call; they are stored as the columns of the
!> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
!> matrix X.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': the linear system involves A;
!>          = 'C': the linear system involves A**H.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of
!>          columns of the matrices B and X. NRHS >= 0.
!>
[in,out]A
!>          A is COMPLEX*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.

Definition at line 180 of file zgels.f.

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

◆ zgelsd()

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 transformations, reducing the original problem
!>     into a  (BLS)
!> (2) Solve the BLS using a divide and conquer approach.
!> (3) Apply back all the Householder transformations to solve
!>     the original least squares problem.
!>
!> The effective rank of A is determined by treating as zero those
!> singular values which are less than RCOND times the largest singular
!> value.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A. M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A. N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X. NRHS >= 0.
!>
[in,out]A
!>          A is COMPLEX*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.
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 223 of file zgelsd.f.

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

◆ zgelss()

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.

Definition at line 176 of file zgelss.f.

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

◆ zgelsx()

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.

Definition at line 182 of file zgelsx.f.

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

◆ zgelsy()

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.
Contributors:
A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

Definition at line 208 of file zgelsy.f.

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

◆ zgesv()

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.

Definition at line 121 of file zgesv.f.

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

◆ zgesvx()

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

Definition at line 347 of file zgesvx.f.

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

◆ zgesvxx()

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 
!>     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  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * 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 . 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 < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * 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 . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.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.

Definition at line 535 of file zgesvxx.f.

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

◆ zgetsls()

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

ZGETSLS

Purpose:
!>
!> ZGETSLS solves overdetermined or underdetermined complex linear systems
!> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
!> factorization of A.  It is assumed that A has full rank.
!>
!>
!>
!> The following options are provided:
!>
!> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A*X ||.
!>
!> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
!>    an underdetermined system A * X = B.
!>
!> 3. If TRANS = 'C' and m >= n:  find the minimum norm solution of
!>    an undetermined system A**T * X = B.
!>
!> 4. If TRANS = 'C' and m < n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A**T * X ||.
!>
!> Several right hand side vectors b and solution vectors x can be
!> handled in a single call; they are stored as the columns of the
!> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
!> matrix X.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': the linear system involves A;
!>          = 'C': the linear system involves A**H.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of
!>          columns of the matrices B and X. NRHS >=0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit,
!>          A is overwritten by details of its QR or LQ
!>          factorization as returned by ZGEQR or ZGELQ.
!>
[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.
!>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
!>          minimum norm solution vectors;
!>          if TRANS = 'C' and m >= n, rows 1 to M of B contain the
!>          minimum norm solution vectors;
!>          if TRANS = 'C' and m < n, rows 1 to M of B contain the
!>          least squares solution vectors.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= MAX(1,M,N).
!>
[out]WORK
!>          (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
!>          or optimal, if query was assumed) LWORK.
!>          See LWORK for details.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If LWORK = -1 or -2, then a workspace query is assumed.
!>          If LWORK = -1, the routine calculates optimal size of WORK for the
!>          optimal performance and returns this value in WORK(1).
!>          If LWORK = -2, the routine calculates minimal size of WORK and 
!>          returns this value in WORK(1).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO =  i, the i-th diagonal element of the
!>                triangular factor of A is zero, so that A does not have
!>                full rank; the least squares solution could not be
!>                computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 160 of file zgetsls.f.

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