LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.
          If M = 0 or N = 0, B is not referenced.
[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.
          If NRHS = 0, RANK = 0 on output.
[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 210 of file zgelsy.f.

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