LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ 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

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```
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,
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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
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
subroutine zgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
ZGEQP3
Definition: zgeqp3.f:159
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:143
subroutine zlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
ZLAIC1 applies one step of incremental condition estimation.
Definition: zlaic1.f:135
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
subroutine zunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRZ
Definition: zunmrz.f:187
subroutine ztzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZTZRZF
Definition: ztzrzf.f:151
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: