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

ZGELSX solves overdetermined or underdetermined systems for GE matrices

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

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

 ZGELSX computes the minimum-norm solution to a complex linear least
 squares problem:
     minimize || A * X - B ||
 using a complete orthogonal factorization of A.  A is an M-by-N
 matrix which may be rank-deficient.

 Several right hand side vectors b and solution vectors x can be
 handled in a single call; they are stored as the columns of the
 M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 matrix X.

 The routine first computes a QR factorization with column pivoting:
     A * P = Q * [ R11 R12 ]
                 [  0  R22 ]
 with R11 defined as the largest leading submatrix whose estimated
 condition number is less than 1/RCOND.  The order of R11, RANK,
 is the effective rank of A.

 Then, R22 is considered to be negligible, and R12 is annihilated
 by unitary transformations from the right, arriving at the
 complete orthogonal factorization:
    A * P = Q * [ T11 0 ] * Z
                [  0  0 ]
 The minimum-norm solution is then
    X = P * Z**H [ inv(T11)*Q1**H*B ]
                 [        0         ]
 where Q1 consists of the first RANK columns of Q.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of
          columns of matrices B and X. NRHS >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, A has been overwritten by details of its
          complete orthogonal factorization.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          On entry, the M-by-NRHS right hand side matrix B.
          On exit, the N-by-NRHS solution matrix X.
          If m >= n and RANK = n, the residual sum-of-squares for
          the solution in the i-th column is given by the sum of
          squares of elements N+1:M in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,M,N).
[in,out]JPVT
          JPVT is INTEGER array, dimension (N)
          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
          initial column, otherwise it is a free column.  Before
          the QR factorization of A, all initial columns are
          permuted to the leading positions; only the remaining
          free columns are moved as a result of column pivoting
          during the factorization.
          On exit, if JPVT(i) = k, then the i-th column of A*P
          was the k-th column of A.
[in]RCOND
          RCOND is DOUBLE PRECISION
          RCOND is used to determine the effective rank of A, which
          is defined as the order of the largest leading triangular
          submatrix R11 in the QR factorization with pivoting of A,
          whose estimated condition number < 1/RCOND.
[out]RANK
          RANK is INTEGER
          The effective rank of A, i.e., the order of the submatrix
          R11.  This is the same as the order of the submatrix T11
          in the complete orthogonal factorization of A.
[out]WORK
          WORK is COMPLEX*16 array, dimension
                      (min(M,N) + max( N, 2*min(M,N)+NRHS )),
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 186 of file zgelsx.f.

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

Here is the call graph for this function: