LAPACK  3.8.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.
Date
December 2016

Definition at line 180 of file dgelsx.f.

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