LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 176 of file dgelsx.f.

178 *
179 * -- LAPACK driver routine --
180 * -- LAPACK is a software package provided by Univ. of Tennessee, --
181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182 *
183 * .. Scalar Arguments ..
184  INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
185  DOUBLE PRECISION RCOND
186 * ..
187 * .. Array Arguments ..
188  INTEGER JPVT( * )
189  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  INTEGER IMAX, IMIN
196  parameter( imax = 1, imin = 2 )
197  DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
198  parameter( zero = 0.0d0, one = 1.0d0, done = zero,
199  $ ntdone = one )
200 * ..
201 * .. Local Scalars ..
202  INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
203  DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
204  $ SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2
205 * ..
206 * .. External Functions ..
207  DOUBLE PRECISION DLAMCH, DLANGE
208  EXTERNAL dlamch, dlange
209 * ..
210 * .. External Subroutines ..
211  EXTERNAL dgeqpf, dlaic1, dlascl, dlaset, dlatzm, dorm2r,
212  $ dtrsm, dtzrqf, xerbla
213 * ..
214 * .. Intrinsic Functions ..
215  INTRINSIC abs, max, min
216 * ..
217 * .. Executable Statements ..
218 *
219  mn = min( m, n )
220  ismin = mn + 1
221  ismax = 2*mn + 1
222 *
223 * Test the input arguments.
224 *
225  info = 0
226  IF( m.LT.0 ) THEN
227  info = -1
228  ELSE IF( n.LT.0 ) THEN
229  info = -2
230  ELSE IF( nrhs.LT.0 ) THEN
231  info = -3
232  ELSE IF( lda.LT.max( 1, m ) ) THEN
233  info = -5
234  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
235  info = -7
236  END IF
237 *
238  IF( info.NE.0 ) THEN
239  CALL xerbla( 'DGELSX', -info )
240  RETURN
241  END IF
242 *
243 * Quick return if possible
244 *
245  IF( min( m, n, nrhs ).EQ.0 ) THEN
246  rank = 0
247  RETURN
248  END IF
249 *
250 * Get machine parameters
251 *
252  smlnum = dlamch( 'S' ) / dlamch( 'P' )
253  bignum = one / smlnum
254  CALL dlabad( smlnum, bignum )
255 *
256 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
257 *
258  anrm = dlange( 'M', m, n, a, lda, work )
259  iascl = 0
260  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
261 *
262 * Scale matrix norm up to SMLNUM
263 *
264  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
265  iascl = 1
266  ELSE IF( anrm.GT.bignum ) THEN
267 *
268 * Scale matrix norm down to BIGNUM
269 *
270  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
271  iascl = 2
272  ELSE IF( anrm.EQ.zero ) THEN
273 *
274 * Matrix all zero. Return zero solution.
275 *
276  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
277  rank = 0
278  GO TO 100
279  END IF
280 *
281  bnrm = dlange( 'M', m, nrhs, b, ldb, work )
282  ibscl = 0
283  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
284 *
285 * Scale matrix norm up to SMLNUM
286 *
287  CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
288  ibscl = 1
289  ELSE IF( bnrm.GT.bignum ) THEN
290 *
291 * Scale matrix norm down to BIGNUM
292 *
293  CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
294  ibscl = 2
295  END IF
296 *
297 * Compute QR factorization with column pivoting of A:
298 * A * P = Q * R
299 *
300  CALL dgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), 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, work( 1 ),
355  $ 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, info )
413  CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
414  $ info )
415  ELSE IF( iascl.EQ.2 ) THEN
416  CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
417  CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
418  $ info )
419  END IF
420  IF( ibscl.EQ.1 ) THEN
421  CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
422  ELSE IF( ibscl.EQ.2 ) THEN
423  CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
424  END IF
425 *
426  100 CONTINUE
427 *
428  RETURN
429 *
430 * End of DGELSX
431 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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:143
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:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:181
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:114
subroutine dgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
DGEQPF
Definition: dgeqpf.f:142
subroutine dlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
DLAIC1 applies one step of incremental condition estimation.
Definition: dlaic1.f:134
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:159
subroutine dtzrqf(M, N, A, LDA, TAU, INFO)
DTZRQF
Definition: dtzrqf.f:138
subroutine dlatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
DLATZM
Definition: dlatzm.f:151
Here is the call graph for this function: