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

◆ dgelsx()

subroutine dgelsx ( integer m,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
double precision rcond,
integer rank,
double precision, dimension( * ) work,
integer info )

DGELSX solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine DGELSY.
!>
!> DGELSX computes the minimum-norm solution to a real 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 orthogonal 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**T [ inv(T11)*Q1**T*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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
!>                      (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
!>
[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 174 of file dgelsx.f.

176*
177* -- LAPACK driver routine --
178* -- LAPACK is a software package provided by Univ. of Tennessee, --
179* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180*
181* .. Scalar Arguments ..
182 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
183 DOUBLE PRECISION RCOND
184* ..
185* .. Array Arguments ..
186 INTEGER JPVT( * )
187 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
188* ..
189*
190* =====================================================================
191*
192* .. Parameters ..
193 INTEGER IMAX, IMIN
194 parameter( imax = 1, imin = 2 )
195 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
196 parameter( zero = 0.0d0, one = 1.0d0, done = zero,
197 $ ntdone = one )
198* ..
199* .. Local Scalars ..
200 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
201 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
202 $ SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2
203* ..
204* .. External Functions ..
205 DOUBLE PRECISION DLAMCH, DLANGE
206 EXTERNAL dlamch, dlange
207* ..
208* .. External Subroutines ..
209 EXTERNAL dgeqpf, dlaic1, dlascl, dlaset, dlatzm, dorm2r,
210 $ dtrsm, dtzrqf, xerbla
211* ..
212* .. Intrinsic Functions ..
213 INTRINSIC abs, max, min
214* ..
215* .. Executable Statements ..
216*
217 mn = min( m, n )
218 ismin = mn + 1
219 ismax = 2*mn + 1
220*
221* Test the input arguments.
222*
223 info = 0
224 IF( m.LT.0 ) THEN
225 info = -1
226 ELSE IF( n.LT.0 ) THEN
227 info = -2
228 ELSE IF( nrhs.LT.0 ) THEN
229 info = -3
230 ELSE IF( lda.LT.max( 1, m ) ) THEN
231 info = -5
232 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
233 info = -7
234 END IF
235*
236 IF( info.NE.0 ) THEN
237 CALL xerbla( 'DGELSX', -info )
238 RETURN
239 END IF
240*
241* Quick return if possible
242*
243 IF( min( m, n, nrhs ).EQ.0 ) THEN
244 rank = 0
245 RETURN
246 END IF
247*
248* Get machine parameters
249*
250 smlnum = dlamch( 'S' ) / dlamch( 'P' )
251 bignum = one / smlnum
252*
253* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
254*
255 anrm = dlange( 'M', m, n, a, lda, work )
256 iascl = 0
257 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
258*
259* Scale matrix norm up to SMLNUM
260*
261 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
262 iascl = 1
263 ELSE IF( anrm.GT.bignum ) THEN
264*
265* Scale matrix norm down to BIGNUM
266*
267 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
268 iascl = 2
269 ELSE IF( anrm.EQ.zero ) THEN
270*
271* Matrix all zero. Return zero solution.
272*
273 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
274 rank = 0
275 GO TO 100
276 END IF
277*
278 bnrm = dlange( 'M', m, nrhs, b, ldb, work )
279 ibscl = 0
280 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
281*
282* Scale matrix norm up to SMLNUM
283*
284 CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
285 $ info )
286 ibscl = 1
287 ELSE IF( bnrm.GT.bignum ) THEN
288*
289* Scale matrix norm down to BIGNUM
290*
291 CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
292 $ info )
293 ibscl = 2
294 END IF
295*
296* Compute QR factorization with column pivoting of A:
297* A * P = Q * R
298*
299 CALL dgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
300 $ info )
301*
302* workspace 3*N. Details of Householder rotations stored
303* in WORK(1:MN).
304*
305* Determine RANK using incremental condition estimation
306*
307 work( ismin ) = one
308 work( ismax ) = one
309 smax = abs( a( 1, 1 ) )
310 smin = smax
311 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
312 rank = 0
313 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
314 GO TO 100
315 ELSE
316 rank = 1
317 END IF
318*
319 10 CONTINUE
320 IF( rank.LT.mn ) THEN
321 i = rank + 1
322 CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
323 $ a( i, i ), sminpr, s1, c1 )
324 CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
325 $ a( i, i ), smaxpr, s2, c2 )
326*
327 IF( smaxpr*rcond.LE.sminpr ) THEN
328 DO 20 i = 1, rank
329 work( ismin+i-1 ) = s1*work( ismin+i-1 )
330 work( ismax+i-1 ) = s2*work( ismax+i-1 )
331 20 CONTINUE
332 work( ismin+rank ) = c1
333 work( ismax+rank ) = c2
334 smin = sminpr
335 smax = smaxpr
336 rank = rank + 1
337 GO TO 10
338 END IF
339 END IF
340*
341* Logically partition R = [ R11 R12 ]
342* [ 0 R22 ]
343* where R11 = R(1:RANK,1:RANK)
344*
345* [R11,R12] = [ T11, 0 ] * Y
346*
347 IF( rank.LT.n )
348 $ CALL dtzrqf( rank, n, a, lda, work( mn+1 ), info )
349*
350* Details of Householder rotations stored in WORK(MN+1:2*MN)
351*
352* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
353*
354 CALL dorm2r( 'Left', 'Transpose', m, nrhs, mn, a, lda,
355 $ work( 1 ), b, ldb, work( 2*mn+1 ), info )
356*
357* workspace NRHS
358*
359* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
360*
361 CALL dtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
362 $ nrhs, one, a, lda, b, ldb )
363*
364 DO 40 i = rank + 1, n
365 DO 30 j = 1, nrhs
366 b( i, j ) = zero
367 30 CONTINUE
368 40 CONTINUE
369*
370* B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
371*
372 IF( rank.LT.n ) THEN
373 DO 50 i = 1, rank
374 CALL dlatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
375 $ work( mn+i ), b( i, 1 ), b( rank+1, 1 ), ldb,
376 $ work( 2*mn+1 ) )
377 50 CONTINUE
378 END IF
379*
380* workspace NRHS
381*
382* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
383*
384 DO 90 j = 1, nrhs
385 DO 60 i = 1, n
386 work( 2*mn+i ) = ntdone
387 60 CONTINUE
388 DO 80 i = 1, n
389 IF( work( 2*mn+i ).EQ.ntdone ) THEN
390 IF( jpvt( i ).NE.i ) THEN
391 k = i
392 t1 = b( k, j )
393 t2 = b( jpvt( k ), j )
394 70 CONTINUE
395 b( jpvt( k ), j ) = t1
396 work( 2*mn+k ) = done
397 t1 = t2
398 k = jpvt( k )
399 t2 = b( jpvt( k ), j )
400 IF( jpvt( k ).NE.i )
401 $ GO TO 70
402 b( i, j ) = t1
403 work( 2*mn+k ) = done
404 END IF
405 END IF
406 80 CONTINUE
407 90 CONTINUE
408*
409* Undo scaling
410*
411 IF( iascl.EQ.1 ) THEN
412 CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
413 $ info )
414 CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
415 $ info )
416 ELSE IF( iascl.EQ.2 ) THEN
417 CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
418 $ info )
419 CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
420 $ info )
421 END IF
422 IF( ibscl.EQ.1 ) THEN
423 CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
424 $ info )
425 ELSE IF( ibscl.EQ.2 ) THEN
426 CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
427 $ info )
428 END IF
429*
430 100 CONTINUE
431*
432 RETURN
433*
434* End of DGELSX
435*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgeqpf
subroutine dgeqpf(m, n, a, lda, jpvt, tau, work, info)
DGEQPF
Definition dgeqpf.f:140
dlatzm
subroutine dlatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
DLATZM
Definition dlatzm.f:150
dtzrqf
subroutine dtzrqf(m, n, a, lda, tau, info)
DTZRQF
Definition dtzrqf.f:136
dlaic1
subroutine dlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
DLAIC1 applies one step of incremental condition estimation.
Definition dlaic1.f:132
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:112
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:142
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
dorm2r
subroutine dorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition dorm2r.f:156
Here is the call graph for this function: