LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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```

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 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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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 zgeqpf(M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
ZGEQPF
Definition: zgeqpf.f:148
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 ztzrqf(M, N, A, LDA, TAU, INFO)
ZTZRQF
Definition: ztzrqf.f:138
subroutine zlatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
ZLATZM
Definition: zlatzm.f:152
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
Here is the call graph for this function: