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

◆ dgelsy()

subroutine dgelsy ( 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 lwork,
integer info )

DGELSY solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> DGELSY 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.
!>
!> This routine is basically identical to the original xGELSX except
!> three differences:
!>   o The call to the subroutine xGEQPF has been substituted by the
!>     the call to the subroutine xGEQP3. This subroutine is a Blas-3
!>     version of the QR factorization with column pivoting.
!>   o Matrix B (the right hand side) is updated with Blas-3.
!>   o The permutation of matrix B (the right hand side) is faster and
!>     more simple.
!>
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 = 0 or N = 0, B is not referenced.
!>
[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 permuted
!>          to the front of AP, otherwise column i is a free column.
!>          On exit, if JPVT(i) = k, then the i-th column of AP
!>          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.
!>          If NRHS = 0, RANK = 0 on output.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          The unblocked strategy requires that:
!>             LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
!>          where MN = min( M, N ).
!>          The block algorithm requires that:
!>             LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
!>          where NB is an upper bound on the blocksize returned
!>          by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
!>          and DORMRZ.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[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.
Contributors:
A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

Definition at line 202 of file dgelsy.f.

205*
206* -- LAPACK driver routine --
207* -- LAPACK is a software package provided by Univ. of Tennessee, --
208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210* .. Scalar Arguments ..
211 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
212 DOUBLE PRECISION RCOND
213* ..
214* .. Array Arguments ..
215 INTEGER JPVT( * )
216 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 INTEGER IMAX, IMIN
223 parameter( imax = 1, imin = 2 )
224 DOUBLE PRECISION ZERO, ONE
225 parameter( zero = 0.0d+0, one = 1.0d+0 )
226* ..
227* .. Local Scalars ..
228 LOGICAL LQUERY
229 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
230 $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4
231 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
232 $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
233* ..
234* .. External Functions ..
235 INTEGER ILAENV
236 DOUBLE PRECISION DLAMCH, DLANGE
237 EXTERNAL ilaenv, dlamch, dlange
238* ..
239* .. External Subroutines ..
240 EXTERNAL dcopy, dgeqp3, dlaic1, dlascl,
241 $ dlaset,
242 $ dormqr, dormrz, dtrsm, dtzrzf, xerbla
243* ..
244* .. Intrinsic Functions ..
245 INTRINSIC abs, max, min
246* ..
247* .. Executable Statements ..
248*
249 mn = min( m, n )
250 ismin = mn + 1
251 ismax = 2*mn + 1
252*
253* Test the input arguments.
254*
255 info = 0
256 lquery = ( lwork.EQ.-1 )
257 IF( m.LT.0 ) THEN
258 info = -1
259 ELSE IF( n.LT.0 ) THEN
260 info = -2
261 ELSE IF( nrhs.LT.0 ) THEN
262 info = -3
263 ELSE IF( lda.LT.max( 1, m ) ) THEN
264 info = -5
265 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
266 info = -7
267 END IF
268*
269* Figure out optimal block size
270*
271 IF( info.EQ.0 ) THEN
272 IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
273 lwkmin = 1
274 lwkopt = 1
275 ELSE
276 nb1 = ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
277 nb2 = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
278 nb3 = ilaenv( 1, 'DORMQR', ' ', m, n, nrhs, -1 )
279 nb4 = ilaenv( 1, 'DORMRQ', ' ', m, n, nrhs, -1 )
280 nb = max( nb1, nb2, nb3, nb4 )
281 lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
282 lwkopt = max( lwkmin,
283 $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
284 END IF
285 work( 1 ) = lwkopt
286*
287 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
288 info = -12
289 END IF
290 END IF
291*
292 IF( info.NE.0 ) THEN
293 CALL xerbla( 'DGELSY', -info )
294 RETURN
295 ELSE IF( lquery ) THEN
296 RETURN
297 END IF
298*
299* Quick return if possible
300*
301 IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
302 rank = 0
303 RETURN
304 END IF
305*
306* Get machine parameters
307*
308 smlnum = dlamch( 'S' ) / dlamch( 'P' )
309 bignum = one / smlnum
310*
311* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
312*
313 anrm = dlange( 'M', m, n, a, lda, work )
314 iascl = 0
315 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
316*
317* Scale matrix norm up to SMLNUM
318*
319 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
320 iascl = 1
321 ELSE IF( anrm.GT.bignum ) THEN
322*
323* Scale matrix norm down to BIGNUM
324*
325 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
326 iascl = 2
327 ELSE IF( anrm.EQ.zero ) THEN
328*
329* Matrix all zero. Return zero solution.
330*
331 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
332 rank = 0
333 GO TO 70
334 END IF
335*
336 bnrm = dlange( 'M', m, nrhs, b, ldb, work )
337 ibscl = 0
338 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
339*
340* Scale matrix norm up to SMLNUM
341*
342 CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
343 $ info )
344 ibscl = 1
345 ELSE IF( bnrm.GT.bignum ) THEN
346*
347* Scale matrix norm down to BIGNUM
348*
349 CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
350 $ info )
351 ibscl = 2
352 END IF
353*
354* Compute QR factorization with column pivoting of A:
355* A * P = Q * R
356*
357 CALL dgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
358 $ lwork-mn, info )
359 wsize = mn + work( mn+1 )
360*
361* workspace: MN+2*N+NB*(N+1).
362* Details of Householder rotations stored in WORK(1:MN).
363*
364* Determine RANK using incremental condition estimation
365*
366 work( ismin ) = one
367 work( ismax ) = one
368 smax = abs( a( 1, 1 ) )
369 smin = smax
370 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
371 rank = 0
372 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
373 GO TO 70
374 ELSE
375 rank = 1
376 END IF
377*
378 10 CONTINUE
379 IF( rank.LT.mn ) THEN
380 i = rank + 1
381 CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
382 $ a( i, i ), sminpr, s1, c1 )
383 CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
384 $ a( i, i ), smaxpr, s2, c2 )
385*
386 IF( smaxpr*rcond.LE.sminpr ) THEN
387 DO 20 i = 1, rank
388 work( ismin+i-1 ) = s1*work( ismin+i-1 )
389 work( ismax+i-1 ) = s2*work( ismax+i-1 )
390 20 CONTINUE
391 work( ismin+rank ) = c1
392 work( ismax+rank ) = c2
393 smin = sminpr
394 smax = smaxpr
395 rank = rank + 1
396 GO TO 10
397 END IF
398 END IF
399*
400* workspace: 3*MN.
401*
402* Logically partition R = [ R11 R12 ]
403* [ 0 R22 ]
404* where R11 = R(1:RANK,1:RANK)
405*
406* [R11,R12] = [ T11, 0 ] * Y
407*
408 IF( rank.LT.n )
409 $ CALL dtzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
410 $ lwork-2*mn, info )
411*
412* workspace: 2*MN.
413* Details of Householder rotations stored in WORK(MN+1:2*MN)
414*
415* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
416*
417 CALL dormqr( 'Left', 'Transpose', m, nrhs, mn, a, lda,
418 $ work( 1 ),
419 $ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
420 wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
421*
422* workspace: 2*MN+NB*NRHS.
423*
424* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
425*
426 CALL dtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
427 $ nrhs, one, a, lda, b, ldb )
428*
429 DO 40 j = 1, nrhs
430 DO 30 i = rank + 1, n
431 b( i, j ) = zero
432 30 CONTINUE
433 40 CONTINUE
434*
435* B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
436*
437 IF( rank.LT.n ) THEN
438 CALL dormrz( 'Left', 'Transpose', n, nrhs, rank, n-rank, a,
439 $ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
440 $ lwork-2*mn, info )
441 END IF
442*
443* workspace: 2*MN+NRHS.
444*
445* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
446*
447 DO 60 j = 1, nrhs
448 DO 50 i = 1, n
449 work( jpvt( i ) ) = b( i, j )
450 50 CONTINUE
451 CALL dcopy( n, work( 1 ), 1, b( 1, j ), 1 )
452 60 CONTINUE
453*
454* workspace: N.
455*
456* Undo scaling
457*
458 IF( iascl.EQ.1 ) THEN
459 CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
460 $ info )
461 CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
462 $ info )
463 ELSE IF( iascl.EQ.2 ) THEN
464 CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
465 $ info )
466 CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
467 $ info )
468 END IF
469 IF( ibscl.EQ.1 ) THEN
470 CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
471 $ info )
472 ELSE IF( ibscl.EQ.2 ) THEN
473 CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
474 $ info )
475 END IF
476*
477 70 CONTINUE
478 work( 1 ) = lwkopt
479*
480 RETURN
481*
482* End of DGELSY
483*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgeqp3
subroutine dgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
DGEQP3
Definition dgeqp3.f:149
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
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
dtzrzf
subroutine dtzrzf(m, n, a, lda, tau, work, lwork, info)
DTZRZF
Definition dtzrzf.f:149
dormqr
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:165
dormrz
subroutine dormrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
DORMRZ
Definition dormrz.f:186
Here is the call graph for this function:
Here is the caller graph for this function: