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

◆ 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,
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*
266* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
267*
268 anrm = zlange( 'M', m, n, a, lda, rwork )
269 iascl = 0
270 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
271*
272* Scale matrix norm up to SMLNUM
273*
274 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
275 iascl = 1
276 ELSE IF( anrm.GT.bignum ) THEN
277*
278* Scale matrix norm down to BIGNUM
279*
280 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
281 iascl = 2
282 ELSE IF( anrm.EQ.zero ) THEN
283*
284* Matrix all zero. Return zero solution.
285*
286 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
287 rank = 0
288 GO TO 100
289 END IF
290*
291 bnrm = zlange( 'M', m, nrhs, b, ldb, rwork )
292 ibscl = 0
293 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
294*
295* Scale matrix norm up to SMLNUM
296*
297 CALL zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
298 ibscl = 1
299 ELSE IF( bnrm.GT.bignum ) THEN
300*
301* Scale matrix norm down to BIGNUM
302*
303 CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
304 ibscl = 2
305 END IF
306*
307* Compute QR factorization with column pivoting of A:
308* A * P = Q * R
309*
310 CALL zgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), rwork,
311 $ info )
312*
313* complex workspace MN+N. Real workspace 2*N. Details of Householder
314* rotations stored in WORK(1:MN).
315*
316* Determine RANK using incremental condition estimation
317*
318 work( ismin ) = cone
319 work( ismax ) = cone
320 smax = abs( a( 1, 1 ) )
321 smin = smax
322 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
323 rank = 0
324 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
325 GO TO 100
326 ELSE
327 rank = 1
328 END IF
329*
330 10 CONTINUE
331 IF( rank.LT.mn ) THEN
332 i = rank + 1
333 CALL zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
334 $ a( i, i ), sminpr, s1, c1 )
335 CALL zlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
336 $ a( i, i ), smaxpr, s2, c2 )
337*
338 IF( smaxpr*rcond.LE.sminpr ) THEN
339 DO 20 i = 1, rank
340 work( ismin+i-1 ) = s1*work( ismin+i-1 )
341 work( ismax+i-1 ) = s2*work( ismax+i-1 )
342 20 CONTINUE
343 work( ismin+rank ) = c1
344 work( ismax+rank ) = c2
345 smin = sminpr
346 smax = smaxpr
347 rank = rank + 1
348 GO TO 10
349 END IF
350 END IF
351*
352* Logically partition R = [ R11 R12 ]
353* [ 0 R22 ]
354* where R11 = R(1:RANK,1:RANK)
355*
356* [R11,R12] = [ T11, 0 ] * Y
357*
358 IF( rank.LT.n )
359 $ CALL ztzrqf( rank, n, a, lda, work( mn+1 ), info )
360*
361* Details of Householder rotations stored in WORK(MN+1:2*MN)
362*
363* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
364*
365 CALL zunm2r( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
366 $ work( 1 ), b, ldb, work( 2*mn+1 ), info )
367*
368* workspace NRHS
369*
370* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
371*
372 CALL ztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
373 $ nrhs, cone, a, lda, b, ldb )
374*
375 DO 40 i = rank + 1, n
376 DO 30 j = 1, nrhs
377 b( i, j ) = czero
378 30 CONTINUE
379 40 CONTINUE
380*
381* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
382*
383 IF( rank.LT.n ) THEN
384 DO 50 i = 1, rank
385 CALL zlatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
386 $ dconjg( work( mn+i ) ), b( i, 1 ),
387 $ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
388 50 CONTINUE
389 END IF
390*
391* workspace NRHS
392*
393* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
394*
395 DO 90 j = 1, nrhs
396 DO 60 i = 1, n
397 work( 2*mn+i ) = ntdone
398 60 CONTINUE
399 DO 80 i = 1, n
400 IF( work( 2*mn+i ).EQ.ntdone ) THEN
401 IF( jpvt( i ).NE.i ) THEN
402 k = i
403 t1 = b( k, j )
404 t2 = b( jpvt( k ), j )
405 70 CONTINUE
406 b( jpvt( k ), j ) = t1
407 work( 2*mn+k ) = done
408 t1 = t2
409 k = jpvt( k )
410 t2 = b( jpvt( k ), j )
411 IF( jpvt( k ).NE.i )
412 $ GO TO 70
413 b( i, j ) = t1
414 work( 2*mn+k ) = done
415 END IF
416 END IF
417 80 CONTINUE
418 90 CONTINUE
419*
420* Undo scaling
421*
422 IF( iascl.EQ.1 ) THEN
423 CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
424 CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
425 $ info )
426 ELSE IF( iascl.EQ.2 ) THEN
427 CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
428 CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
429 $ info )
430 END IF
431 IF( ibscl.EQ.1 ) THEN
432 CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
433 ELSE IF( ibscl.EQ.2 ) THEN
434 CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
435 END IF
436*
437 100 CONTINUE
438*
439 RETURN
440*
441* End of ZGELSX
442*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 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
subroutine zgeqpf(m, n, a, lda, jpvt, tau, work, rwork, info)
ZGEQPF
Definition zgeqpf.f:148
subroutine zlatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
ZLATZM
Definition zlatzm.f:152
subroutine ztzrqf(m, n, a, lda, tau, info)
ZTZRQF
Definition ztzrqf.f:138
Here is the call graph for this function: