OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
double
LAPACK » General Matrices » Linear Solve

Functions

subroutine dgelsx (m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, info)
  DGELSX solves overdetermined or underdetermined systems for GE matrices
subroutine dgels (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
  DGELS solves overdetermined or underdetermined systems for GE matrices
subroutine dgelsd (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info)
  DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
subroutine dgelss (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, info)
  DGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine dgelsy (m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)
  DGELSY solves overdetermined or underdetermined systems for GE matrices
subroutine dgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
  DGESV computes the solution to system of linear equations A * X = B for GE matrices
subroutine dgesvx (fact, trans, n, nrhs, a, lda, af, ldaf, ipiv, equed, r, c, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
  DGESVX computes the solution to system of linear equations A * X = B for GE matrices
subroutine dgesvxx (fact, trans, n, nrhs, a, lda, af, ldaf, ipiv, equed, r, c, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
  DGESVXX computes the solution to system of linear equations A * X = B for GE matrices
subroutine dgetsls (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
 DGETSLS
subroutine dsgesv (n, nrhs, a, lda, ipiv, b, ldb, x, ldx, work, swork, iter, info)
  DSGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement)

Detailed Description

This is the group of double solve driver functions for GE matrices

Function Documentation

◆ dgels()

subroutine dgels ( character trans,
integer m,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) work,
integer lwork,
integer info )

DGELS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> DGELS solves overdetermined or underdetermined real linear systems
!> involving an M-by-N matrix A, or its transpose, using a QR or LQ
!> factorization of A.  It is assumed that A has full rank.
!>
!> The following options are provided:
!>
!> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A*X ||.
!>
!> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
!>    an underdetermined system A * X = B.
!>
!> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
!>    an underdetermined system A**T * X = B.
!>
!> 4. If TRANS = 'T' and m < n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A**T * X ||.
!>
!> 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.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': the linear system involves A;
!>          = 'T': the linear system involves A**T.
!>
[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 the 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,
!>            if M >= N, A is overwritten by details of its QR
!>                       factorization as returned by DGEQRF;
!>            if M <  N, A is overwritten by details of its LQ
!>                       factorization as returned by DGELQF.
!>
[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 matrix B of right hand side vectors, stored
!>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
!>          if TRANS = 'T'.
!>          On exit, if INFO = 0, B is overwritten by the solution
!>          vectors, stored columnwise:
!>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
!>          squares solution vectors; the residual sum of squares for the
!>          solution in each column is given by the sum of squares of
!>          elements N+1 to M in that column;
!>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
!>          minimum norm solution vectors;
!>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
!>          minimum norm solution vectors;
!>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
!>          least squares solution vectors; the residual sum of squares
!>          for the solution in each column is given by the sum of
!>          squares of elements M+1 to N in that column.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= MAX(1,M,N).
!>
[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.
!>          LWORK >= max( 1, MN + max( MN, NRHS ) ).
!>          For optimal performance,
!>          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
!>          where MN = min(M,N) and NB is the optimum block size.
!>
!>          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
!>          > 0:  if INFO =  i, the i-th diagonal element of the
!>                triangular factor of A is zero, so that A does not have
!>                full rank; the least squares solution could not be
!>                computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 181 of file dgels.f.

183*
184* -- LAPACK driver routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER TRANS
190 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
191* ..
192* .. Array Arguments ..
193 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 DOUBLE PRECISION ZERO, ONE
200 parameter( zero = 0.0d0, one = 1.0d0 )
201* ..
202* .. Local Scalars ..
203 LOGICAL LQUERY, TPSD
204 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
205 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
206* ..
207* .. Local Arrays ..
208 DOUBLE PRECISION RWORK( 1 )
209* ..
210* .. External Functions ..
211 LOGICAL LSAME
212 INTEGER ILAENV
213 DOUBLE PRECISION DLAMCH, DLANGE
214 EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
215* ..
216* .. External Subroutines ..
217 EXTERNAL dgelqf, dgeqrf, dlascl, dlaset, dormlq, dormqr,
218 $ dtrtrs, xerbla
219* ..
220* .. Intrinsic Functions ..
221 INTRINSIC dble, max, min
222* ..
223* .. Executable Statements ..
224*
225* Test the input arguments.
226*
227 info = 0
228 mn = min( m, n )
229 lquery = ( lwork.EQ.-1 )
230 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'T' ) ) ) THEN
231 info = -1
232 ELSE IF( m.LT.0 ) THEN
233 info = -2
234 ELSE IF( n.LT.0 ) THEN
235 info = -3
236 ELSE IF( nrhs.LT.0 ) THEN
237 info = -4
238 ELSE IF( lda.LT.max( 1, m ) ) THEN
239 info = -6
240 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
241 info = -8
242 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND. .NOT.lquery )
243 $ THEN
244 info = -10
245 END IF
246*
247* Figure out optimal block size
248*
249 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
250*
251 tpsd = .true.
252 IF( lsame( trans, 'N' ) )
253 $ tpsd = .false.
254*
255 IF( m.GE.n ) THEN
256 nb = ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
257 IF( tpsd ) THEN
258 nb = max( nb, ilaenv( 1, 'DORMQR', 'LN', m, nrhs, n,
259 $ -1 ) )
260 ELSE
261 nb = max( nb, ilaenv( 1, 'DORMQR', 'LT', m, nrhs, n,
262 $ -1 ) )
263 END IF
264 ELSE
265 nb = ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
266 IF( tpsd ) THEN
267 nb = max( nb, ilaenv( 1, 'DORMLQ', 'LT', n, nrhs, m,
268 $ -1 ) )
269 ELSE
270 nb = max( nb, ilaenv( 1, 'DORMLQ', 'LN', n, nrhs, m,
271 $ -1 ) )
272 END IF
273 END IF
274*
275 wsize = max( 1, mn+max( mn, nrhs )*nb )
276 work( 1 ) = dble( wsize )
277*
278 END IF
279*
280 IF( info.NE.0 ) THEN
281 CALL xerbla( 'DGELS ', -info )
282 RETURN
283 ELSE IF( lquery ) THEN
284 RETURN
285 END IF
286*
287* Quick return if possible
288*
289 IF( min( m, n, nrhs ).EQ.0 ) THEN
290 CALL dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
291 RETURN
292 END IF
293*
294* Get machine parameters
295*
296 smlnum = dlamch( 'S' ) / dlamch( 'P' )
297 bignum = one / smlnum
298 CALL dlabad( smlnum, bignum )
299*
300* Scale A, B if max element outside range [SMLNUM,BIGNUM]
301*
302 anrm = dlange( 'M', m, n, a, lda, rwork )
303 iascl = 0
304 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
305*
306* Scale matrix norm up to SMLNUM
307*
308 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
309 iascl = 1
310 ELSE IF( anrm.GT.bignum ) THEN
311*
312* Scale matrix norm down to BIGNUM
313*
314 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
315 iascl = 2
316 ELSE IF( anrm.EQ.zero ) THEN
317*
318* Matrix all zero. Return zero solution.
319*
320 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
321 GO TO 50
322 END IF
323*
324 brow = m
325 IF( tpsd )
326 $ brow = n
327 bnrm = dlange( 'M', brow, nrhs, b, ldb, rwork )
328 ibscl = 0
329 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
330*
331* Scale matrix norm up to SMLNUM
332*
333 CALL dlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
334 $ info )
335 ibscl = 1
336 ELSE IF( bnrm.GT.bignum ) THEN
337*
338* Scale matrix norm down to BIGNUM
339*
340 CALL dlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
341 $ info )
342 ibscl = 2
343 END IF
344*
345 IF( m.GE.n ) THEN
346*
347* compute QR factorization of A
348*
349 CALL dgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
350 $ info )
351*
352* workspace at least N, optimally N*NB
353*
354 IF( .NOT.tpsd ) THEN
355*
356* Least-Squares Problem min || A * X - B ||
357*
358* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
359*
360 CALL dormqr( 'Left', 'Transpose', m, nrhs, n, a, lda,
361 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
362 $ info )
363*
364* workspace at least NRHS, optimally NRHS*NB
365*
366* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
367*
368 CALL dtrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
369 $ a, lda, b, ldb, info )
370*
371 IF( info.GT.0 ) THEN
372 RETURN
373 END IF
374*
375 scllen = n
376*
377 ELSE
378*
379* Underdetermined system of equations A**T * X = B
380*
381* B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
382*
383 CALL dtrtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
384 $ a, lda, b, ldb, info )
385*
386 IF( info.GT.0 ) THEN
387 RETURN
388 END IF
389*
390* B(N+1:M,1:NRHS) = ZERO
391*
392 DO 20 j = 1, nrhs
393 DO 10 i = n + 1, m
394 b( i, j ) = zero
395 10 CONTINUE
396 20 CONTINUE
397*
398* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
399*
400 CALL dormqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
401 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
402 $ info )
403*
404* workspace at least NRHS, optimally NRHS*NB
405*
406 scllen = m
407*
408 END IF
409*
410 ELSE
411*
412* Compute LQ factorization of A
413*
414 CALL dgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
415 $ info )
416*
417* workspace at least M, optimally M*NB.
418*
419 IF( .NOT.tpsd ) THEN
420*
421* underdetermined system of equations A * X = B
422*
423* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
424*
425 CALL dtrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
426 $ a, lda, b, ldb, info )
427*
428 IF( info.GT.0 ) THEN
429 RETURN
430 END IF
431*
432* B(M+1:N,1:NRHS) = 0
433*
434 DO 40 j = 1, nrhs
435 DO 30 i = m + 1, n
436 b( i, j ) = zero
437 30 CONTINUE
438 40 CONTINUE
439*
440* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
441*
442 CALL dormlq( 'Left', 'Transpose', n, nrhs, m, a, lda,
443 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
444 $ info )
445*
446* workspace at least NRHS, optimally NRHS*NB
447*
448 scllen = n
449*
450 ELSE
451*
452* overdetermined system min || A**T * X - B ||
453*
454* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
455*
456 CALL dormlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
457 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
458 $ info )
459*
460* workspace at least NRHS, optimally NRHS*NB
461*
462* B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
463*
464 CALL dtrtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
465 $ a, lda, b, ldb, info )
466*
467 IF( info.GT.0 ) THEN
468 RETURN
469 END IF
470*
471 scllen = m
472*
473 END IF
474*
475 END IF
476*
477* Undo scaling
478*
479 IF( iascl.EQ.1 ) THEN
480 CALL dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
481 $ info )
482 ELSE IF( iascl.EQ.2 ) THEN
483 CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
484 $ info )
485 END IF
486 IF( ibscl.EQ.1 ) THEN
487 CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
488 $ info )
489 ELSE IF( ibscl.EQ.2 ) THEN
490 CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
491 $ info )
492 END IF
493*
494 50 CONTINUE
495 work( 1 ) = dble( wsize )
496*
497 RETURN
498*
499* End of DGELS
500*
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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:143
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:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
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:114
dgelqf
subroutine dgelqf(m, n, a, lda, tau, work, lwork, info)
DGELQF
Definition dgelqf.f:143
dgeqrf
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
dormqr
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
dtrtrs
subroutine dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
DTRTRS
Definition dtrtrs.f:140
dormlq
subroutine dormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMLQ
Definition dormlq.f:167
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ dgelsd()

subroutine dgelsd ( integer m,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) s,
double precision rcond,
integer rank,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer info )

DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices

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

Purpose:
!>
!> DGELSD computes the minimum-norm solution to a real linear least
!> squares problem:
!>     minimize 2-norm(| b - A*x |)
!> using the singular value decomposition (SVD) 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 problem is solved in three steps:
!> (1) Reduce the coefficient matrix A to bidiagonal form with
!>     Householder transformations, reducing the original problem
!>     into a  (BLS)
!> (2) Solve the BLS using a divide and conquer approach.
!> (3) Apply back all the Householder transformations to solve
!>     the original least squares problem.
!>
!> The effective rank of A is determined by treating as zero those
!> singular values which are less than RCOND times the largest singular
!> value.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of A. M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of A. N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the 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 destroyed.
!>
[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, B is overwritten by 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,max(M,N)).
!>
[out]S
!>          S is DOUBLE PRECISION array, dimension (min(M,N))
!>          The singular values of A in decreasing order.
!>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
!>
[in]RCOND
!>          RCOND is DOUBLE PRECISION
!>          RCOND is used to determine the effective rank of A.
!>          Singular values S(i) <= RCOND*S(1) are treated as zero.
!>          If RCOND < 0, machine precision is used instead.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the number of singular values
!>          which are greater than RCOND*S(1).
!>
[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. LWORK must be at least 1.
!>          The exact minimum amount of workspace needed depends on M,
!>          N and NRHS. As long as LWORK is at least
!>              12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
!>          if M is greater than or equal to N or
!>              12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
!>          if M is less than N, the code will execute correctly.
!>          SMLSIZ is returned by ILAENV and is equal to the maximum
!>          size of the subproblems at the bottom of the computation
!>          tree (usually about 25), and
!>             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
!>          For good performance, LWORK should generally be larger.
!>
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN),
!>          where MINMN = MIN( M,N ).
!>          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  the algorithm for computing the SVD failed to converge;
!>                if INFO = i, i off-diagonal elements of an intermediate
!>                bidiagonal form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 207 of file dgelsd.f.

209*
210* -- LAPACK driver routine --
211* -- LAPACK is a software package provided by Univ. of Tennessee, --
212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214* .. Scalar Arguments ..
215 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
216 DOUBLE PRECISION RCOND
217* ..
218* .. Array Arguments ..
219 INTEGER IWORK( * )
220 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
221* ..
222*
223* =====================================================================
224*
225* .. Parameters ..
226 DOUBLE PRECISION ZERO, ONE, TWO
227 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
228* ..
229* .. Local Scalars ..
230 LOGICAL LQUERY
231 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
232 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
233 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
234 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
235* ..
236* .. External Subroutines ..
237 EXTERNAL dgebrd, dgelqf, dgeqrf, dlabad, dlacpy, dlalsd,
238 $ dlascl, dlaset, dormbr, dormlq, dormqr, xerbla
239* ..
240* .. External Functions ..
241 INTEGER ILAENV
242 DOUBLE PRECISION DLAMCH, DLANGE
243 EXTERNAL ilaenv, dlamch, dlange
244* ..
245* .. Intrinsic Functions ..
246 INTRINSIC dble, int, log, max, min
247* ..
248* .. Executable Statements ..
249*
250* Test the input arguments.
251*
252 info = 0
253 minmn = min( m, n )
254 maxmn = max( m, n )
255 mnthr = ilaenv( 6, 'DGELSD', ' ', m, n, nrhs, -1 )
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, maxmn ) ) THEN
266 info = -7
267 END IF
268*
269 smlsiz = ilaenv( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
270*
271* Compute workspace.
272* (Note: Comments in the code beginning "Workspace:" describe the
273* minimal amount of workspace needed at that point in the code,
274* as well as the preferred amount for good performance.
275* NB refers to the optimal block size for the immediately
276* following subroutine, as returned by ILAENV.)
277*
278 minwrk = 1
279 liwork = 1
280 minmn = max( 1, minmn )
281 nlvl = max( int( log( dble( minmn ) / dble( smlsiz+1 ) ) /
282 $ log( two ) ) + 1, 0 )
283*
284 IF( info.EQ.0 ) THEN
285 maxwrk = 0
286 liwork = 3*minmn*nlvl + 11*minmn
287 mm = m
288 IF( m.GE.n .AND. m.GE.mnthr ) THEN
289*
290* Path 1a - overdetermined, with many more rows than columns.
291*
292 mm = n
293 maxwrk = max( maxwrk, n+n*ilaenv( 1, 'DGEQRF', ' ', m, n,
294 $ -1, -1 ) )
295 maxwrk = max( maxwrk, n+nrhs*
296 $ ilaenv( 1, 'DORMQR', 'LT', m, nrhs, n, -1 ) )
297 END IF
298 IF( m.GE.n ) THEN
299*
300* Path 1 - overdetermined or exactly determined.
301*
302 maxwrk = max( maxwrk, 3*n+( mm+n )*
303 $ ilaenv( 1, 'DGEBRD', ' ', mm, n, -1, -1 ) )
304 maxwrk = max( maxwrk, 3*n+nrhs*
305 $ ilaenv( 1, 'DORMBR', 'QLT', mm, nrhs, n, -1 ) )
306 maxwrk = max( maxwrk, 3*n+( n-1 )*
307 $ ilaenv( 1, 'DORMBR', 'PLN', n, nrhs, n, -1 ) )
308 wlalsd = 9*n+2*n*smlsiz+8*n*nlvl+n*nrhs+(smlsiz+1)**2
309 maxwrk = max( maxwrk, 3*n+wlalsd )
310 minwrk = max( 3*n+mm, 3*n+nrhs, 3*n+wlalsd )
311 END IF
312 IF( n.GT.m ) THEN
313 wlalsd = 9*m+2*m*smlsiz+8*m*nlvl+m*nrhs+(smlsiz+1)**2
314 IF( n.GE.mnthr ) THEN
315*
316* Path 2a - underdetermined, with many more columns
317* than rows.
318*
319 maxwrk = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
320 maxwrk = max( maxwrk, m*m+4*m+2*m*
321 $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
322 maxwrk = max( maxwrk, m*m+4*m+nrhs*
323 $ ilaenv( 1, 'DORMBR', 'QLT', m, nrhs, m, -1 ) )
324 maxwrk = max( maxwrk, m*m+4*m+( m-1 )*
325 $ ilaenv( 1, 'DORMBR', 'PLN', m, nrhs, m, -1 ) )
326 IF( nrhs.GT.1 ) THEN
327 maxwrk = max( maxwrk, m*m+m+m*nrhs )
328 ELSE
329 maxwrk = max( maxwrk, m*m+2*m )
330 END IF
331 maxwrk = max( maxwrk, m+nrhs*
332 $ ilaenv( 1, 'DORMLQ', 'LT', n, nrhs, m, -1 ) )
333 maxwrk = max( maxwrk, m*m+4*m+wlalsd )
334! XXX: Ensure the Path 2a case below is triggered. The workspace
335! calculation should use queries for all routines eventually.
336 maxwrk = max( maxwrk,
337 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
338 ELSE
339*
340* Path 2 - remaining underdetermined cases.
341*
342 maxwrk = 3*m + ( n+m )*ilaenv( 1, 'DGEBRD', ' ', m, n,
343 $ -1, -1 )
344 maxwrk = max( maxwrk, 3*m+nrhs*
345 $ ilaenv( 1, 'DORMBR', 'QLT', m, nrhs, n, -1 ) )
346 maxwrk = max( maxwrk, 3*m+m*
347 $ ilaenv( 1, 'DORMBR', 'PLN', n, nrhs, m, -1 ) )
348 maxwrk = max( maxwrk, 3*m+wlalsd )
349 END IF
350 minwrk = max( 3*m+nrhs, 3*m+m, 3*m+wlalsd )
351 END IF
352 minwrk = min( minwrk, maxwrk )
353 work( 1 ) = maxwrk
354 iwork( 1 ) = liwork
355
356 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
357 info = -12
358 END IF
359 END IF
360*
361 IF( info.NE.0 ) THEN
362 CALL xerbla( 'DGELSD', -info )
363 RETURN
364 ELSE IF( lquery ) THEN
365 GO TO 10
366 END IF
367*
368* Quick return if possible.
369*
370 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
371 rank = 0
372 RETURN
373 END IF
374*
375* Get machine parameters.
376*
377 eps = dlamch( 'P' )
378 sfmin = dlamch( 'S' )
379 smlnum = sfmin / eps
380 bignum = one / smlnum
381 CALL dlabad( smlnum, bignum )
382*
383* Scale A if max entry outside range [SMLNUM,BIGNUM].
384*
385 anrm = dlange( 'M', m, n, a, lda, work )
386 iascl = 0
387 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
388*
389* Scale matrix norm up to SMLNUM.
390*
391 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
392 iascl = 1
393 ELSE IF( anrm.GT.bignum ) THEN
394*
395* Scale matrix norm down to BIGNUM.
396*
397 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
398 iascl = 2
399 ELSE IF( anrm.EQ.zero ) THEN
400*
401* Matrix all zero. Return zero solution.
402*
403 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
404 CALL dlaset( 'F', minmn, 1, zero, zero, s, 1 )
405 rank = 0
406 GO TO 10
407 END IF
408*
409* Scale B if max entry outside range [SMLNUM,BIGNUM].
410*
411 bnrm = dlange( 'M', m, nrhs, b, ldb, work )
412 ibscl = 0
413 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
414*
415* Scale matrix norm up to SMLNUM.
416*
417 CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
418 ibscl = 1
419 ELSE IF( bnrm.GT.bignum ) THEN
420*
421* Scale matrix norm down to BIGNUM.
422*
423 CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
424 ibscl = 2
425 END IF
426*
427* If M < N make sure certain entries of B are zero.
428*
429 IF( m.LT.n )
430 $ CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
431*
432* Overdetermined case.
433*
434 IF( m.GE.n ) THEN
435*
436* Path 1 - overdetermined or exactly determined.
437*
438 mm = m
439 IF( m.GE.mnthr ) THEN
440*
441* Path 1a - overdetermined, with many more rows than columns.
442*
443 mm = n
444 itau = 1
445 nwork = itau + n
446*
447* Compute A=Q*R.
448* (Workspace: need 2*N, prefer N+N*NB)
449*
450 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
451 $ lwork-nwork+1, info )
452*
453* Multiply B by transpose(Q).
454* (Workspace: need N+NRHS, prefer N+NRHS*NB)
455*
456 CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
457 $ ldb, work( nwork ), lwork-nwork+1, info )
458*
459* Zero out below R.
460*
461 IF( n.GT.1 ) THEN
462 CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
463 END IF
464 END IF
465*
466 ie = 1
467 itauq = ie + n
468 itaup = itauq + n
469 nwork = itaup + n
470*
471* Bidiagonalize R in A.
472* (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
473*
474 CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
475 $ work( itaup ), work( nwork ), lwork-nwork+1,
476 $ info )
477*
478* Multiply B by transpose of left bidiagonalizing vectors of R.
479* (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
480*
481 CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
482 $ b, ldb, work( nwork ), lwork-nwork+1, info )
483*
484* Solve the bidiagonal least squares problem.
485*
486 CALL dlalsd( 'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
487 $ rcond, rank, work( nwork ), iwork, info )
488 IF( info.NE.0 ) THEN
489 GO TO 10
490 END IF
491*
492* Multiply B by right bidiagonalizing vectors of R.
493*
494 CALL dormbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
495 $ b, ldb, work( nwork ), lwork-nwork+1, info )
496*
497 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
498 $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
499*
500* Path 2a - underdetermined, with many more columns than rows
501* and sufficient workspace for an efficient algorithm.
502*
503 ldwork = m
504 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
505 $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
506 itau = 1
507 nwork = m + 1
508*
509* Compute A=L*Q.
510* (Workspace: need 2*M, prefer M+M*NB)
511*
512 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
513 $ lwork-nwork+1, info )
514 il = nwork
515*
516* Copy L to WORK(IL), zeroing out above its diagonal.
517*
518 CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwork )
519 CALL dlaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
520 $ ldwork )
521 ie = il + ldwork*m
522 itauq = ie + m
523 itaup = itauq + m
524 nwork = itaup + m
525*
526* Bidiagonalize L in WORK(IL).
527* (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
528*
529 CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
530 $ work( itauq ), work( itaup ), work( nwork ),
531 $ lwork-nwork+1, info )
532*
533* Multiply B by transpose of left bidiagonalizing vectors of L.
534* (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
535*
536 CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
537 $ work( itauq ), b, ldb, work( nwork ),
538 $ lwork-nwork+1, info )
539*
540* Solve the bidiagonal least squares problem.
541*
542 CALL dlalsd( 'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
543 $ rcond, rank, work( nwork ), iwork, info )
544 IF( info.NE.0 ) THEN
545 GO TO 10
546 END IF
547*
548* Multiply B by right bidiagonalizing vectors of L.
549*
550 CALL dormbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
551 $ work( itaup ), b, ldb, work( nwork ),
552 $ lwork-nwork+1, info )
553*
554* Zero out below first M rows of B.
555*
556 CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
557 nwork = itau + m
558*
559* Multiply transpose(Q) by B.
560* (Workspace: need M+NRHS, prefer M+NRHS*NB)
561*
562 CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
563 $ ldb, work( nwork ), lwork-nwork+1, info )
564*
565 ELSE
566*
567* Path 2 - remaining underdetermined cases.
568*
569 ie = 1
570 itauq = ie + m
571 itaup = itauq + m
572 nwork = itaup + m
573*
574* Bidiagonalize A.
575* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
576*
577 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
578 $ work( itaup ), work( nwork ), lwork-nwork+1,
579 $ info )
580*
581* Multiply B by transpose of left bidiagonalizing vectors.
582* (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
583*
584 CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
585 $ b, ldb, work( nwork ), lwork-nwork+1, info )
586*
587* Solve the bidiagonal least squares problem.
588*
589 CALL dlalsd( 'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
590 $ rcond, rank, work( nwork ), iwork, info )
591 IF( info.NE.0 ) THEN
592 GO TO 10
593 END IF
594*
595* Multiply B by right bidiagonalizing vectors of A.
596*
597 CALL dormbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
598 $ b, ldb, work( nwork ), lwork-nwork+1, info )
599*
600 END IF
601*
602* Undo scaling.
603*
604 IF( iascl.EQ.1 ) THEN
605 CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
606 CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
607 $ info )
608 ELSE IF( iascl.EQ.2 ) THEN
609 CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
610 CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
611 $ info )
612 END IF
613 IF( ibscl.EQ.1 ) THEN
614 CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
615 ELSE IF( ibscl.EQ.2 ) THEN
616 CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
617 END IF
618*
619 10 CONTINUE
620 work( 1 ) = maxwrk
621 iwork( 1 ) = liwork
622 RETURN
623*
624* End of DGELSD
625*
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dgebrd
subroutine dgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
DGEBRD
Definition dgebrd.f:205
dormbr
subroutine dormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMBR
Definition dormbr.f:195
dlalsd
subroutine dlalsd(uplo, smlsiz, n, nrhs, d, e, b, ldb, rcond, rank, work, iwork, info)
DLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition dlalsd.f:179

◆ dgelss()

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

DGELSS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> DGELSS computes the minimum norm solution to a real linear least
!> squares problem:
!>
!> Minimize 2-norm(| b - A*x |).
!>
!> using the singular value decomposition (SVD) 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 effective rank of A is determined by treating as zero those
!> singular values which are less than RCOND times the largest singular
!> value.
!>
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 the 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, the first min(m,n) rows of A are overwritten with
!>          its right singular vectors, stored rowwise.
!>
[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, B is overwritten by 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,max(M,N)).
!>
[out]S
!>          S is DOUBLE PRECISION array, dimension (min(M,N))
!>          The singular values of A in decreasing order.
!>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
!>
[in]RCOND
!>          RCOND is DOUBLE PRECISION
!>          RCOND is used to determine the effective rank of A.
!>          Singular values S(i) <= RCOND*S(1) are treated as zero.
!>          If RCOND < 0, machine precision is used instead.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the number of singular values
!>          which are greater than RCOND*S(1).
!>
[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. LWORK >= 1, and also:
!>          LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
!>          For good performance, LWORK should generally be larger.
!>
!>          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.
!>          > 0:  the algorithm for computing the SVD failed to converge;
!>                if INFO = i, i off-diagonal elements of an intermediate
!>                bidiagonal form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 170 of file dgelss.f.

172*
173* -- LAPACK driver routine --
174* -- LAPACK is a software package provided by Univ. of Tennessee, --
175* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176*
177* .. Scalar Arguments ..
178 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
179 DOUBLE PRECISION RCOND
180* ..
181* .. Array Arguments ..
182 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
183* ..
184*
185* =====================================================================
186*
187* .. Parameters ..
188 DOUBLE PRECISION ZERO, ONE
189 parameter( zero = 0.0d+0, one = 1.0d+0 )
190* ..
191* .. Local Scalars ..
192 LOGICAL LQUERY
193 INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
194 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
195 $ MAXWRK, MINMN, MINWRK, MM, MNTHR
196 INTEGER LWORK_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD,
197 $ LWORK_DORMBR, LWORK_DORGBR, LWORK_DORMLQ,
198 $ LWORK_DGELQF
199 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
200* ..
201* .. Local Arrays ..
202 DOUBLE PRECISION DUM( 1 )
203* ..
204* .. External Subroutines ..
205 EXTERNAL dbdsqr, dcopy, dgebrd, dgelqf, dgemm, dgemv,
206 $ dgeqrf, dlabad, dlacpy, dlascl, dlaset, dorgbr,
207 $ dormbr, dormlq, dormqr, drscl, xerbla
208* ..
209* .. External Functions ..
210 INTEGER ILAENV
211 DOUBLE PRECISION DLAMCH, DLANGE
212 EXTERNAL ilaenv, dlamch, dlange
213* ..
214* .. Intrinsic Functions ..
215 INTRINSIC max, min
216* ..
217* .. Executable Statements ..
218*
219* Test the input arguments
220*
221 info = 0
222 minmn = min( m, n )
223 maxmn = max( m, n )
224 lquery = ( lwork.EQ.-1 )
225 IF( m.LT.0 ) THEN
226 info = -1
227 ELSE IF( n.LT.0 ) THEN
228 info = -2
229 ELSE IF( nrhs.LT.0 ) THEN
230 info = -3
231 ELSE IF( lda.LT.max( 1, m ) ) THEN
232 info = -5
233 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
234 info = -7
235 END IF
236*
237* Compute workspace
238* (Note: Comments in the code beginning "Workspace:" describe the
239* minimal amount of workspace needed at that point in the code,
240* as well as the preferred amount for good performance.
241* NB refers to the optimal block size for the immediately
242* following subroutine, as returned by ILAENV.)
243*
244 IF( info.EQ.0 ) THEN
245 minwrk = 1
246 maxwrk = 1
247 IF( minmn.GT.0 ) THEN
248 mm = m
249 mnthr = ilaenv( 6, 'DGELSS', ' ', m, n, nrhs, -1 )
250 IF( m.GE.n .AND. m.GE.mnthr ) THEN
251*
252* Path 1a - overdetermined, with many more rows than
253* columns
254*
255* Compute space needed for DGEQRF
256 CALL dgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
257 lwork_dgeqrf=dum(1)
258* Compute space needed for DORMQR
259 CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, dum(1), b,
260 $ ldb, dum(1), -1, info )
261 lwork_dormqr=dum(1)
262 mm = n
263 maxwrk = max( maxwrk, n + lwork_dgeqrf )
264 maxwrk = max( maxwrk, n + lwork_dormqr )
265 END IF
266 IF( m.GE.n ) THEN
267*
268* Path 1 - overdetermined or exactly determined
269*
270* Compute workspace needed for DBDSQR
271*
272 bdspac = max( 1, 5*n )
273* Compute space needed for DGEBRD
274 CALL dgebrd( mm, n, a, lda, s, dum(1), dum(1),
275 $ dum(1), dum(1), -1, info )
276 lwork_dgebrd=dum(1)
277* Compute space needed for DORMBR
278 CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, dum(1),
279 $ b, ldb, dum(1), -1, info )
280 lwork_dormbr=dum(1)
281* Compute space needed for DORGBR
282 CALL dorgbr( 'P', n, n, n, a, lda, dum(1),
283 $ dum(1), -1, info )
284 lwork_dorgbr=dum(1)
285* Compute total workspace needed
286 maxwrk = max( maxwrk, 3*n + lwork_dgebrd )
287 maxwrk = max( maxwrk, 3*n + lwork_dormbr )
288 maxwrk = max( maxwrk, 3*n + lwork_dorgbr )
289 maxwrk = max( maxwrk, bdspac )
290 maxwrk = max( maxwrk, n*nrhs )
291 minwrk = max( 3*n + mm, 3*n + nrhs, bdspac )
292 maxwrk = max( minwrk, maxwrk )
293 END IF
294 IF( n.GT.m ) THEN
295*
296* Compute workspace needed for DBDSQR
297*
298 bdspac = max( 1, 5*m )
299 minwrk = max( 3*m+nrhs, 3*m+n, bdspac )
300 IF( n.GE.mnthr ) THEN
301*
302* Path 2a - underdetermined, with many more columns
303* than rows
304*
305* Compute space needed for DGELQF
306 CALL dgelqf( m, n, a, lda, dum(1), dum(1),
307 $ -1, info )
308 lwork_dgelqf=dum(1)
309* Compute space needed for DGEBRD
310 CALL dgebrd( m, m, a, lda, s, dum(1), dum(1),
311 $ dum(1), dum(1), -1, info )
312 lwork_dgebrd=dum(1)
313* Compute space needed for DORMBR
314 CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda,
315 $ dum(1), b, ldb, dum(1), -1, info )
316 lwork_dormbr=dum(1)
317* Compute space needed for DORGBR
318 CALL dorgbr( 'P', m, m, m, a, lda, dum(1),
319 $ dum(1), -1, info )
320 lwork_dorgbr=dum(1)
321* Compute space needed for DORMLQ
322 CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, dum(1),
323 $ b, ldb, dum(1), -1, info )
324 lwork_dormlq=dum(1)
325* Compute total workspace needed
326 maxwrk = m + lwork_dgelqf
327 maxwrk = max( maxwrk, m*m + 4*m + lwork_dgebrd )
328 maxwrk = max( maxwrk, m*m + 4*m + lwork_dormbr )
329 maxwrk = max( maxwrk, m*m + 4*m + lwork_dorgbr )
330 maxwrk = max( maxwrk, m*m + m + bdspac )
331 IF( nrhs.GT.1 ) THEN
332 maxwrk = max( maxwrk, m*m + m + m*nrhs )
333 ELSE
334 maxwrk = max( maxwrk, m*m + 2*m )
335 END IF
336 maxwrk = max( maxwrk, m + lwork_dormlq )
337 ELSE
338*
339* Path 2 - underdetermined
340*
341* Compute space needed for DGEBRD
342 CALL dgebrd( m, n, a, lda, s, dum(1), dum(1),
343 $ dum(1), dum(1), -1, info )
344 lwork_dgebrd=dum(1)
345* Compute space needed for DORMBR
346 CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, a, lda,
347 $ dum(1), b, ldb, dum(1), -1, info )
348 lwork_dormbr=dum(1)
349* Compute space needed for DORGBR
350 CALL dorgbr( 'P', m, n, m, a, lda, dum(1),
351 $ dum(1), -1, info )
352 lwork_dorgbr=dum(1)
353 maxwrk = 3*m + lwork_dgebrd
354 maxwrk = max( maxwrk, 3*m + lwork_dormbr )
355 maxwrk = max( maxwrk, 3*m + lwork_dorgbr )
356 maxwrk = max( maxwrk, bdspac )
357 maxwrk = max( maxwrk, n*nrhs )
358 END IF
359 END IF
360 maxwrk = max( minwrk, maxwrk )
361 END IF
362 work( 1 ) = maxwrk
363*
364 IF( lwork.LT.minwrk .AND. .NOT.lquery )
365 $ info = -12
366 END IF
367*
368 IF( info.NE.0 ) THEN
369 CALL xerbla( 'DGELSS', -info )
370 RETURN
371 ELSE IF( lquery ) THEN
372 RETURN
373 END IF
374*
375* Quick return if possible
376*
377 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
378 rank = 0
379 RETURN
380 END IF
381*
382* Get machine parameters
383*
384 eps = dlamch( 'P' )
385 sfmin = dlamch( 'S' )
386 smlnum = sfmin / eps
387 bignum = one / smlnum
388 CALL dlabad( smlnum, bignum )
389*
390* Scale A if max element outside range [SMLNUM,BIGNUM]
391*
392 anrm = dlange( 'M', m, n, a, lda, work )
393 iascl = 0
394 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
395*
396* Scale matrix norm up to SMLNUM
397*
398 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
399 iascl = 1
400 ELSE IF( anrm.GT.bignum ) THEN
401*
402* Scale matrix norm down to BIGNUM
403*
404 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
405 iascl = 2
406 ELSE IF( anrm.EQ.zero ) THEN
407*
408* Matrix all zero. Return zero solution.
409*
410 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
411 CALL dlaset( 'F', minmn, 1, zero, zero, s, minmn )
412 rank = 0
413 GO TO 70
414 END IF
415*
416* Scale B if max element outside range [SMLNUM,BIGNUM]
417*
418 bnrm = dlange( 'M', m, nrhs, b, ldb, work )
419 ibscl = 0
420 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
421*
422* Scale matrix norm up to SMLNUM
423*
424 CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
425 ibscl = 1
426 ELSE IF( bnrm.GT.bignum ) THEN
427*
428* Scale matrix norm down to BIGNUM
429*
430 CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
431 ibscl = 2
432 END IF
433*
434* Overdetermined case
435*
436 IF( m.GE.n ) THEN
437*
438* Path 1 - overdetermined or exactly determined
439*
440 mm = m
441 IF( m.GE.mnthr ) THEN
442*
443* Path 1a - overdetermined, with many more rows than columns
444*
445 mm = n
446 itau = 1
447 iwork = itau + n
448*
449* Compute A=Q*R
450* (Workspace: need 2*N, prefer N+N*NB)
451*
452 CALL dgeqrf( m, n, a, lda, work( itau ), work( iwork ),
453 $ lwork-iwork+1, info )
454*
455* Multiply B by transpose(Q)
456* (Workspace: need N+NRHS, prefer N+NRHS*NB)
457*
458 CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
459 $ ldb, work( iwork ), lwork-iwork+1, info )
460*
461* Zero out below R
462*
463 IF( n.GT.1 )
464 $ CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
465 END IF
466*
467 ie = 1
468 itauq = ie + n
469 itaup = itauq + n
470 iwork = itaup + n
471*
472* Bidiagonalize R in A
473* (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
474*
475 CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
476 $ work( itaup ), work( iwork ), lwork-iwork+1,
477 $ info )
478*
479* Multiply B by transpose of left bidiagonalizing vectors of R
480* (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
481*
482 CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
483 $ b, ldb, work( iwork ), lwork-iwork+1, info )
484*
485* Generate right bidiagonalizing vectors of R in A
486* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
487*
488 CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
489 $ work( iwork ), lwork-iwork+1, info )
490 iwork = ie + n
491*
492* Perform bidiagonal QR iteration
493* multiply B by transpose of left singular vectors
494* compute right singular vectors in A
495* (Workspace: need BDSPAC)
496*
497 CALL dbdsqr( 'U', n, n, 0, nrhs, s, work( ie ), a, lda, dum,
498 $ 1, b, ldb, work( iwork ), info )
499 IF( info.NE.0 )
500 $ GO TO 70
501*
502* Multiply B by reciprocals of singular values
503*
504 thr = max( rcond*s( 1 ), sfmin )
505 IF( rcond.LT.zero )
506 $ thr = max( eps*s( 1 ), sfmin )
507 rank = 0
508 DO 10 i = 1, n
509 IF( s( i ).GT.thr ) THEN
510 CALL drscl( nrhs, s( i ), b( i, 1 ), ldb )
511 rank = rank + 1
512 ELSE
513 CALL dlaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
514 END IF
515 10 CONTINUE
516*
517* Multiply B by right singular vectors
518* (Workspace: need N, prefer N*NRHS)
519*
520 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
521 CALL dgemm( 'T', 'N', n, nrhs, n, one, a, lda, b, ldb, zero,
522 $ work, ldb )
523 CALL dlacpy( 'G', n, nrhs, work, ldb, b, ldb )
524 ELSE IF( nrhs.GT.1 ) THEN
525 chunk = lwork / n
526 DO 20 i = 1, nrhs, chunk
527 bl = min( nrhs-i+1, chunk )
528 CALL dgemm( 'T', 'N', n, bl, n, one, a, lda, b( 1, i ),
529 $ ldb, zero, work, n )
530 CALL dlacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
531 20 CONTINUE
532 ELSE
533 CALL dgemv( 'T', n, n, one, a, lda, b, 1, zero, work, 1 )
534 CALL dcopy( n, work, 1, b, 1 )
535 END IF
536*
537 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
538 $ max( m, 2*m-4, nrhs, n-3*m ) ) THEN
539*
540* Path 2a - underdetermined, with many more columns than rows
541* and sufficient workspace for an efficient algorithm
542*
543 ldwork = m
544 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
545 $ m*lda+m+m*nrhs ) )ldwork = lda
546 itau = 1
547 iwork = m + 1
548*
549* Compute A=L*Q
550* (Workspace: need 2*M, prefer M+M*NB)
551*
552 CALL dgelqf( m, n, a, lda, work( itau ), work( iwork ),
553 $ lwork-iwork+1, info )
554 il = iwork
555*
556* Copy L to WORK(IL), zeroing out above it
557*
558 CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwork )
559 CALL dlaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
560 $ ldwork )
561 ie = il + ldwork*m
562 itauq = ie + m
563 itaup = itauq + m
564 iwork = itaup + m
565*
566* Bidiagonalize L in WORK(IL)
567* (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
568*
569 CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
570 $ work( itauq ), work( itaup ), work( iwork ),
571 $ lwork-iwork+1, info )
572*
573* Multiply B by transpose of left bidiagonalizing vectors of L
574* (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
575*
576 CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
577 $ work( itauq ), b, ldb, work( iwork ),
578 $ lwork-iwork+1, info )
579*
580* Generate right bidiagonalizing vectors of R in WORK(IL)
581* (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
582*
583 CALL dorgbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
584 $ work( iwork ), lwork-iwork+1, info )
585 iwork = ie + m
586*
587* Perform bidiagonal QR iteration,
588* computing right singular vectors of L in WORK(IL) and
589* multiplying B by transpose of left singular vectors
590* (Workspace: need M*M+M+BDSPAC)
591*
592 CALL dbdsqr( 'U', m, m, 0, nrhs, s, work( ie ), work( il ),
593 $ ldwork, a, lda, b, ldb, work( iwork ), info )
594 IF( info.NE.0 )
595 $ GO TO 70
596*
597* Multiply B by reciprocals of singular values
598*
599 thr = max( rcond*s( 1 ), sfmin )
600 IF( rcond.LT.zero )
601 $ thr = max( eps*s( 1 ), sfmin )
602 rank = 0
603 DO 30 i = 1, m
604 IF( s( i ).GT.thr ) THEN
605 CALL drscl( nrhs, s( i ), b( i, 1 ), ldb )
606 rank = rank + 1
607 ELSE
608 CALL dlaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
609 END IF
610 30 CONTINUE
611 iwork = ie
612*
613* Multiply B by right singular vectors of L in WORK(IL)
614* (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
615*
616 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
617 CALL dgemm( 'T', 'N', m, nrhs, m, one, work( il ), ldwork,
618 $ b, ldb, zero, work( iwork ), ldb )
619 CALL dlacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
620 ELSE IF( nrhs.GT.1 ) THEN
621 chunk = ( lwork-iwork+1 ) / m
622 DO 40 i = 1, nrhs, chunk
623 bl = min( nrhs-i+1, chunk )
624 CALL dgemm( 'T', 'N', m, bl, m, one, work( il ), ldwork,
625 $ b( 1, i ), ldb, zero, work( iwork ), m )
626 CALL dlacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
627 $ ldb )
628 40 CONTINUE
629 ELSE
630 CALL dgemv( 'T', m, m, one, work( il ), ldwork, b( 1, 1 ),
631 $ 1, zero, work( iwork ), 1 )
632 CALL dcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
633 END IF
634*
635* Zero out below first M rows of B
636*
637 CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
638 iwork = itau + m
639*
640* Multiply transpose(Q) by B
641* (Workspace: need M+NRHS, prefer M+NRHS*NB)
642*
643 CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
644 $ ldb, work( iwork ), lwork-iwork+1, info )
645*
646 ELSE
647*
648* Path 2 - remaining underdetermined cases
649*
650 ie = 1
651 itauq = ie + m
652 itaup = itauq + m
653 iwork = itaup + m
654*
655* Bidiagonalize A
656* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
657*
658 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
659 $ work( itaup ), work( iwork ), lwork-iwork+1,
660 $ info )
661*
662* Multiply B by transpose of left bidiagonalizing vectors
663* (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
664*
665 CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
666 $ b, ldb, work( iwork ), lwork-iwork+1, info )
667*
668* Generate right bidiagonalizing vectors in A
669* (Workspace: need 4*M, prefer 3*M+M*NB)
670*
671 CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
672 $ work( iwork ), lwork-iwork+1, info )
673 iwork = ie + m
674*
675* Perform bidiagonal QR iteration,
676* computing right singular vectors of A in A and
677* multiplying B by transpose of left singular vectors
678* (Workspace: need BDSPAC)
679*
680 CALL dbdsqr( 'L', m, n, 0, nrhs, s, work( ie ), a, lda, dum,
681 $ 1, b, ldb, work( iwork ), info )
682 IF( info.NE.0 )
683 $ GO TO 70
684*
685* Multiply B by reciprocals of singular values
686*
687 thr = max( rcond*s( 1 ), sfmin )
688 IF( rcond.LT.zero )
689 $ thr = max( eps*s( 1 ), sfmin )
690 rank = 0
691 DO 50 i = 1, m
692 IF( s( i ).GT.thr ) THEN
693 CALL drscl( nrhs, s( i ), b( i, 1 ), ldb )
694 rank = rank + 1
695 ELSE
696 CALL dlaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
697 END IF
698 50 CONTINUE
699*
700* Multiply B by right singular vectors of A
701* (Workspace: need N, prefer N*NRHS)
702*
703 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
704 CALL dgemm( 'T', 'N', n, nrhs, m, one, a, lda, b, ldb, zero,
705 $ work, ldb )
706 CALL dlacpy( 'F', n, nrhs, work, ldb, b, ldb )
707 ELSE IF( nrhs.GT.1 ) THEN
708 chunk = lwork / n
709 DO 60 i = 1, nrhs, chunk
710 bl = min( nrhs-i+1, chunk )
711 CALL dgemm( 'T', 'N', n, bl, m, one, a, lda, b( 1, i ),
712 $ ldb, zero, work, n )
713 CALL dlacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
714 60 CONTINUE
715 ELSE
716 CALL dgemv( 'T', m, n, one, a, lda, b, 1, zero, work, 1 )
717 CALL dcopy( n, work, 1, b, 1 )
718 END IF
719 END IF
720*
721* Undo scaling
722*
723 IF( iascl.EQ.1 ) THEN
724 CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
725 CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
726 $ info )
727 ELSE IF( iascl.EQ.2 ) THEN
728 CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
729 CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
730 $ info )
731 END IF
732 IF( ibscl.EQ.1 ) THEN
733 CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
734 ELSE IF( ibscl.EQ.2 ) THEN
735 CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
736 END IF
737*
738 70 CONTINUE
739 work( 1 ) = maxwrk
740 RETURN
741*
742* End of DGELSS
743*
dbdsqr
subroutine dbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DBDSQR
Definition dbdsqr.f:241
dorgbr
subroutine dorgbr(vect, m, n, k, a, lda, tau, work, lwork, info)
DORGBR
Definition dorgbr.f:157
drscl
subroutine drscl(n, sa, sx, incx)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition drscl.f:84
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187

◆ 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*
dgeqpf
subroutine dgeqpf(m, n, a, lda, jpvt, tau, work, info)
DGEQPF
Definition dgeqpf.f:142
dlaic1
subroutine dlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
DLAIC1 applies one step of incremental condition estimation.
Definition dlaic1.f:134
dtzrqf
subroutine dtzrqf(m, n, a, lda, tau, info)
DTZRQF
Definition dtzrqf.f:138
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:159
dlatzm
subroutine dlatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
DLATZM
Definition dlatzm.f:151
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181

◆ 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.
!>
[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.
!>
[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.

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

◆ dgesv()

subroutine dgesv ( integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
double precision, dimension( ldb, * ) b,
integer ldb,
integer info )

DGESV computes the solution to system of linear equations A * X = B for GE matrices

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

Purpose:
!>
!> DGESV computes the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!> The LU decomposition with partial pivoting and row interchanges is
!> used to factor A as
!>    A = P * L * U,
!> where P is a permutation matrix, L is unit lower triangular, and U is
!> upper triangular.  The factored form of A is then used to solve the
!> system of equations A * X = B.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the N-by-N coefficient matrix A.
!>          On exit, the factors L and U from the factorization
!>          A = P*L*U; the unit diagonal elements of L are not stored.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices that define the permutation matrix P;
!>          row i of the matrix was interchanged with row IPIV(i).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS matrix of right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
!>                has been completed, but the factor U is exactly
!>                singular, so the solution could not be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 121 of file dgesv.f.

122*
123* -- LAPACK driver routine --
124* -- LAPACK is a software package provided by Univ. of Tennessee, --
125* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126*
127* .. Scalar Arguments ..
128 INTEGER INFO, LDA, LDB, N, NRHS
129* ..
130* .. Array Arguments ..
131 INTEGER IPIV( * )
132 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
133* ..
134*
135* =====================================================================
136*
137* .. External Subroutines ..
138 EXTERNAL dgetrf, dgetrs, xerbla
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC max
142* ..
143* .. Executable Statements ..
144*
145* Test the input parameters.
146*
147 info = 0
148 IF( n.LT.0 ) THEN
149 info = -1
150 ELSE IF( nrhs.LT.0 ) THEN
151 info = -2
152 ELSE IF( lda.LT.max( 1, n ) ) THEN
153 info = -4
154 ELSE IF( ldb.LT.max( 1, n ) ) THEN
155 info = -7
156 END IF
157 IF( info.NE.0 ) THEN
158 CALL xerbla( 'DGESV ', -info )
159 RETURN
160 END IF
161*
162* Compute the LU factorization of A.
163*
164 CALL dgetrf( n, n, a, lda, ipiv, info )
165 IF( info.EQ.0 ) THEN
166*
167* Solve the system A*X = B, overwriting B with X.
168*
169 CALL dgetrs( 'No transpose', n, nrhs, a, lda, ipiv, b, ldb,
170 $ info )
171 END IF
172 RETURN
173*
174* End of DGESV
175*
dgetrs
subroutine dgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
DGETRS
Definition dgetrs.f:121
dgetrf
subroutine dgetrf(m, n, a, lda, ipiv, info)
DGETRF
Definition dgetrf.f:108

◆ dgesvx()

subroutine dgesvx ( character fact,
character trans,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldaf, * ) af,
integer ldaf,
integer, dimension( * ) ipiv,
character equed,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldx, * ) x,
integer ldx,
double precision rcond,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

DGESVX computes the solution to system of linear equations A * X = B for GE matrices

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

Purpose:
!>
!> DGESVX uses the LU factorization to compute the solution to a real
!> system of linear equations
!>    A * X = B,
!> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>       TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
!>       TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
!>       TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
!>    or diag(C)*B (if TRANS = 'T' or 'C').
!>
!> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
!>    matrix A (after equilibration if FACT = 'E') as
!>       A = P * L * U,
!>    where P is a permutation matrix, L is a unit lower triangular
!>    matrix, and U is upper triangular.
!>
!> 3. If some U(i,i)=0, so that U is exactly singular, then the routine
!>    returns with INFO = i. Otherwise, the factored form of A is used
!>    to estimate the condition number of the matrix A.  If the
!>    reciprocal of the condition number is less than machine precision,
!>    INFO = N+1 is returned as a warning, but the routine still goes on
!>    to solve for X and compute error bounds as described below.
!>
!> 4. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 5. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
!> 6. If equilibration was used, the matrix X is premultiplied by
!>    diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
!>    that it solves the original system before equilibration.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AF and IPIV contain the factored form of A.
!>                  If EQUED is not 'N', the matrix A has been
!>                  equilibrated with scaling factors given by R and C.
!>                  A, AF, and IPIV are not modified.
!>          = 'N':  The matrix A will be copied to AF and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AF and factored.
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations:
!>          = 'N':  A * X = B     (No transpose)
!>          = 'T':  A**T * X = B  (Transpose)
!>          = 'C':  A**H * X = B  (Transpose)
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the N-by-N matrix A.  If FACT = 'F' and EQUED is
!>          not 'N', then A must have been equilibrated by the scaling
!>          factors in R and/or C.  A is not modified if FACT = 'F' or
!>          'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>          On exit, if EQUED .ne. 'N', A is scaled as follows:
!>          EQUED = 'R':  A := diag(R) * A
!>          EQUED = 'C':  A := A * diag(C)
!>          EQUED = 'B':  A := diag(R) * A * diag(C).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]AF
!>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
!>          If FACT = 'F', then AF is an input argument and on entry
!>          contains the factors L and U from the factorization
!>          A = P*L*U as computed by DGETRF.  If EQUED .ne. 'N', then
!>          AF is the factored form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AF is an output argument and on exit
!>          returns the factors L and U from the factorization A = P*L*U
!>          of the original matrix A.
!>
!>          If FACT = 'E', then AF is an output argument and on exit
!>          returns the factors L and U from the factorization A = P*L*U
!>          of the equilibrated matrix A (see the description of A for
!>          the form of the equilibrated matrix).
!>
[in]LDAF
!>          LDAF is INTEGER
!>          The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          If FACT = 'F', then IPIV is an input argument and on entry
!>          contains the pivot indices from the factorization A = P*L*U
!>          as computed by DGETRF; row i of the matrix was interchanged
!>          with row IPIV(i).
!>
!>          If FACT = 'N', then IPIV is an output argument and on exit
!>          contains the pivot indices from the factorization A = P*L*U
!>          of the original matrix A.
!>
!>          If FACT = 'E', then IPIV is an output argument and on exit
!>          contains the pivot indices from the factorization A = P*L*U
!>          of the equilibrated matrix A.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'R':  Row equilibration, i.e., A has been premultiplied by
!>                  diag(R).
!>          = 'C':  Column equilibration, i.e., A has been postmultiplied
!>                  by diag(C).
!>          = 'B':  Both row and column equilibration, i.e., A has been
!>                  replaced by diag(R) * A * diag(C).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]R
!>          R is DOUBLE PRECISION array, dimension (N)
!>          The row scale factors for A.  If EQUED = 'R' or 'B', A is
!>          multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
!>          is not accessed.  R is an input argument if FACT = 'F';
!>          otherwise, R is an output argument.  If FACT = 'F' and
!>          EQUED = 'R' or 'B', each element of R must be positive.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>          The column scale factors for A.  If EQUED = 'C' or 'B', A is
!>          multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
!>          is not accessed.  C is an input argument if FACT = 'F';
!>          otherwise, C is an output argument.  If FACT = 'F' and
!>          EQUED = 'C' or 'B', each element of C must be positive.
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit,
!>          if EQUED = 'N', B is not modified;
!>          if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
!>          diag(R)*B;
!>          if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
!>          overwritten by diag(C)*B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
!>          to the original system of equations.  Note that A and B are
!>          modified on exit if EQUED .ne. 'N', and the solution to the
!>          equilibrated system is inv(diag(C))*X if TRANS = 'N' and
!>          EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
!>          and EQUED = 'R' or 'B'.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (4*N)
!>          On exit, WORK(1) contains the reciprocal pivot growth
!>          factor norm(A)/norm(U). The  norm is
!>          used. If WORK(1) is much less than 1, then the stability
!>          of the LU factorization of the (equilibrated) matrix A
!>          could be poor. This also means that the solution X, condition
!>          estimator RCOND, and forward error bound FERR could be
!>          unreliable. If factorization fails with 0<INFO<=N, then
!>          WORK(1) contains the reciprocal pivot growth factor for the
!>          leading INFO columns of A.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  U(i,i) is exactly zero.  The factorization has
!>                       been completed, but the factor U is exactly
!>                       singular, so the solution and error bounds
!>                       could not be computed. RCOND = 0 is returned.
!>                = N+1: U is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 346 of file dgesvx.f.

349*
350* -- LAPACK driver routine --
351* -- LAPACK is a software package provided by Univ. of Tennessee, --
352* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
353*
354* .. Scalar Arguments ..
355 CHARACTER EQUED, FACT, TRANS
356 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
357 DOUBLE PRECISION RCOND
358* ..
359* .. Array Arguments ..
360 INTEGER IPIV( * ), IWORK( * )
361 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
362 $ BERR( * ), C( * ), FERR( * ), R( * ),
363 $ WORK( * ), X( LDX, * )
364* ..
365*
366* =====================================================================
367*
368* .. Parameters ..
369 DOUBLE PRECISION ZERO, ONE
370 parameter( zero = 0.0d+0, one = 1.0d+0 )
371* ..
372* .. Local Scalars ..
373 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
374 CHARACTER NORM
375 INTEGER I, INFEQU, J
376 DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
377 $ ROWCND, RPVGRW, SMLNUM
378* ..
379* .. External Functions ..
380 LOGICAL LSAME
381 DOUBLE PRECISION DLAMCH, DLANGE, DLANTR
382 EXTERNAL lsame, dlamch, dlange, dlantr
383* ..
384* .. External Subroutines ..
385 EXTERNAL dgecon, dgeequ, dgerfs, dgetrf, dgetrs, dlacpy,
386 $ dlaqge, xerbla
387* ..
388* .. Intrinsic Functions ..
389 INTRINSIC max, min
390* ..
391* .. Executable Statements ..
392*
393 info = 0
394 nofact = lsame( fact, 'N' )
395 equil = lsame( fact, 'E' )
396 notran = lsame( trans, 'N' )
397 IF( nofact .OR. equil ) THEN
398 equed = 'N'
399 rowequ = .false.
400 colequ = .false.
401 ELSE
402 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
403 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
404 smlnum = dlamch( 'Safe minimum' )
405 bignum = one / smlnum
406 END IF
407*
408* Test the input parameters.
409*
410 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
411 $ THEN
412 info = -1
413 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
414 $ lsame( trans, 'C' ) ) THEN
415 info = -2
416 ELSE IF( n.LT.0 ) THEN
417 info = -3
418 ELSE IF( nrhs.LT.0 ) THEN
419 info = -4
420 ELSE IF( lda.LT.max( 1, n ) ) THEN
421 info = -6
422 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
423 info = -8
424 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
425 $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
426 info = -10
427 ELSE
428 IF( rowequ ) THEN
429 rcmin = bignum
430 rcmax = zero
431 DO 10 j = 1, n
432 rcmin = min( rcmin, r( j ) )
433 rcmax = max( rcmax, r( j ) )
434 10 CONTINUE
435 IF( rcmin.LE.zero ) THEN
436 info = -11
437 ELSE IF( n.GT.0 ) THEN
438 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
439 ELSE
440 rowcnd = one
441 END IF
442 END IF
443 IF( colequ .AND. info.EQ.0 ) THEN
444 rcmin = bignum
445 rcmax = zero
446 DO 20 j = 1, n
447 rcmin = min( rcmin, c( j ) )
448 rcmax = max( rcmax, c( j ) )
449 20 CONTINUE
450 IF( rcmin.LE.zero ) THEN
451 info = -12
452 ELSE IF( n.GT.0 ) THEN
453 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
454 ELSE
455 colcnd = one
456 END IF
457 END IF
458 IF( info.EQ.0 ) THEN
459 IF( ldb.LT.max( 1, n ) ) THEN
460 info = -14
461 ELSE IF( ldx.LT.max( 1, n ) ) THEN
462 info = -16
463 END IF
464 END IF
465 END IF
466*
467 IF( info.NE.0 ) THEN
468 CALL xerbla( 'DGESVX', -info )
469 RETURN
470 END IF
471*
472 IF( equil ) THEN
473*
474* Compute row and column scalings to equilibrate the matrix A.
475*
476 CALL dgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
477 IF( infequ.EQ.0 ) THEN
478*
479* Equilibrate the matrix.
480*
481 CALL dlaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
482 $ equed )
483 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
484 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
485 END IF
486 END IF
487*
488* Scale the right hand side.
489*
490 IF( notran ) THEN
491 IF( rowequ ) THEN
492 DO 40 j = 1, nrhs
493 DO 30 i = 1, n
494 b( i, j ) = r( i )*b( i, j )
495 30 CONTINUE
496 40 CONTINUE
497 END IF
498 ELSE IF( colequ ) THEN
499 DO 60 j = 1, nrhs
500 DO 50 i = 1, n
501 b( i, j ) = c( i )*b( i, j )
502 50 CONTINUE
503 60 CONTINUE
504 END IF
505*
506 IF( nofact .OR. equil ) THEN
507*
508* Compute the LU factorization of A.
509*
510 CALL dlacpy( 'Full', n, n, a, lda, af, ldaf )
511 CALL dgetrf( n, n, af, ldaf, ipiv, info )
512*
513* Return if INFO is non-zero.
514*
515 IF( info.GT.0 ) THEN
516*
517* Compute the reciprocal pivot growth factor of the
518* leading rank-deficient INFO columns of A.
519*
520 rpvgrw = dlantr( 'M', 'U', 'N', info, info, af, ldaf,
521 $ work )
522 IF( rpvgrw.EQ.zero ) THEN
523 rpvgrw = one
524 ELSE
525 rpvgrw = dlange( 'M', n, info, a, lda, work ) / rpvgrw
526 END IF
527 work( 1 ) = rpvgrw
528 rcond = zero
529 RETURN
530 END IF
531 END IF
532*
533* Compute the norm of the matrix A and the
534* reciprocal pivot growth factor RPVGRW.
535*
536 IF( notran ) THEN
537 norm = '1'
538 ELSE
539 norm = 'I'
540 END IF
541 anorm = dlange( norm, n, n, a, lda, work )
542 rpvgrw = dlantr( 'M', 'U', 'N', n, n, af, ldaf, work )
543 IF( rpvgrw.EQ.zero ) THEN
544 rpvgrw = one
545 ELSE
546 rpvgrw = dlange( 'M', n, n, a, lda, work ) / rpvgrw
547 END IF
548*
549* Compute the reciprocal of the condition number of A.
550*
551 CALL dgecon( norm, n, af, ldaf, anorm, rcond, work, iwork, info )
552*
553* Compute the solution matrix X.
554*
555 CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
556 CALL dgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
557*
558* Use iterative refinement to improve the computed solution and
559* compute error bounds and backward error estimates for it.
560*
561 CALL dgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
562 $ ldx, ferr, berr, work, iwork, info )
563*
564* Transform the solution matrix X to a solution of the original
565* system.
566*
567 IF( notran ) THEN
568 IF( colequ ) THEN
569 DO 80 j = 1, nrhs
570 DO 70 i = 1, n
571 x( i, j ) = c( i )*x( i, j )
572 70 CONTINUE
573 80 CONTINUE
574 DO 90 j = 1, nrhs
575 ferr( j ) = ferr( j ) / colcnd
576 90 CONTINUE
577 END IF
578 ELSE IF( rowequ ) THEN
579 DO 110 j = 1, nrhs
580 DO 100 i = 1, n
581 x( i, j ) = r( i )*x( i, j )
582 100 CONTINUE
583 110 CONTINUE
584 DO 120 j = 1, nrhs
585 ferr( j ) = ferr( j ) / rowcnd
586 120 CONTINUE
587 END IF
588*
589 work( 1 ) = rpvgrw
590*
591* Set INFO = N+1 if the matrix is singular to working precision.
592*
593 IF( rcond.LT.dlamch( 'Epsilon' ) )
594 $ info = n + 1
595 RETURN
596*
597* End of DGESVX
598*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
dlaqge
subroutine dlaqge(m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
DLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
Definition dlaqge.f:142
dgecon
subroutine dgecon(norm, n, a, lda, anorm, rcond, work, iwork, info)
DGECON
Definition dgecon.f:124
dgeequ
subroutine dgeequ(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
DGEEQU
Definition dgeequ.f:139
dgerfs
subroutine dgerfs(trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
DGERFS
Definition dgerfs.f:185
dlantr
double precision function dlantr(norm, uplo, diag, m, n, a, lda, work)
DLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlantr.f:141

◆ dgesvxx()

subroutine dgesvxx ( character fact,
character trans,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldaf, * ) af,
integer ldaf,
integer, dimension( * ) ipiv,
character equed,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldx , * ) x,
integer ldx,
double precision rcond,
double precision rpvgrw,
double precision, dimension( * ) berr,
integer n_err_bnds,
double precision, dimension( nrhs, * ) err_bnds_norm,
double precision, dimension( nrhs, * ) err_bnds_comp,
integer nparams,
double precision, dimension( * ) params,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

DGESVXX computes the solution to system of linear equations A * X = B for GE matrices

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

Purpose:
!>
!>    DGESVXX uses the LU factorization to compute the solution to a
!>    double precision system of linear equations  A * X = B,  where A is an
!>    N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!>    If requested, both normwise and maximum componentwise error bounds
!>    are returned. DGESVXX will return a solution with a tiny
!>    guaranteed error (O(eps) where eps is the working machine
!>    precision) unless the matrix is very ill-conditioned, in which
!>    case a warning is returned. Relevant condition numbers also are
!>    calculated and returned.
!>
!>    DGESVXX accepts user-provided factorizations and equilibration
!>    factors; see the definitions of the FACT and EQUED options.
!>    Solving with refinement and using a factorization from a previous
!>    DGESVXX call will also produce a solution with either O(eps)
!>    errors or warnings, but we cannot make that claim for general
!>    user-provided factorizations and equilibration factors if they
!>    differ from what DGESVXX would itself produce.
!>
Description:
!>
!>    The following steps are performed:
!>
!>    1. If FACT = 'E', double precision scaling factors are computed to equilibrate
!>    the system:
!>
!>      TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
!>      TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
!>      TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
!>
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
!>    or diag(C)*B (if TRANS = 'T' or 'C').
!>
!>    2. If FACT = 'N' or 'E', the LU decomposition is used to factor
!>    the matrix A (after equilibration if FACT = 'E') as
!>
!>      A = P * L * U,
!>
!>    where P is a permutation matrix, L is a unit lower triangular
!>    matrix, and U is upper triangular.
!>
!>    3. If some U(i,i)=0, so that U is exactly singular, then the
!>    routine returns with INFO = i. Otherwise, the factored form of A
!>    is used to estimate the condition number of the matrix A (see
!>    argument RCOND). If the reciprocal of the condition number is less
!>    than machine precision, the routine still goes on to solve for X
!>    and compute error bounds as described below.
!>
!>    4. The system of equations is solved for X using the factored form
!>    of A.
!>
!>    5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
!>    the routine will use iterative refinement to try to get a small
!>    error and error bounds.  Refinement calculates the residual to at
!>    least twice the working precision.
!>
!>    6. If equilibration was used, the matrix X is premultiplied by
!>    diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
!>    that it solves the original system before equilibration.
!> 
!>     Some optional parameters are bundled in the PARAMS array.  These
!>     settings determine how refinement is performed, but often the
!>     defaults are acceptable.  If the defaults are acceptable, users
!>     can pass NPARAMS = 0 which prevents the source code from accessing
!>     the PARAMS argument.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>     Specifies whether or not the factored form of the matrix A is
!>     supplied on entry, and if not, whether the matrix A should be
!>     equilibrated before it is factored.
!>       = 'F':  On entry, AF and IPIV contain the factored form of A.
!>               If EQUED is not 'N', the matrix A has been
!>               equilibrated with scaling factors given by R and C.
!>               A, AF, and IPIV are not modified.
!>       = 'N':  The matrix A will be copied to AF and factored.
!>       = 'E':  The matrix A will be equilibrated if necessary, then
!>               copied to AF and factored.
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>     Specifies the form of the system of equations:
!>       = 'N':  A * X = B     (No transpose)
!>       = 'T':  A**T * X = B  (Transpose)
!>       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
!>
[in]N
!>          N is INTEGER
!>     The number of linear equations, i.e., the order of the
!>     matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>     The number of right hand sides, i.e., the number of columns
!>     of the matrices B and X.  NRHS >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>     On entry, the N-by-N matrix A.  If FACT = 'F' and EQUED is
!>     not 'N', then A must have been equilibrated by the scaling
!>     factors in R and/or C.  A is not modified if FACT = 'F' or
!>     'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>     On exit, if EQUED .ne. 'N', A is scaled as follows:
!>     EQUED = 'R':  A := diag(R) * A
!>     EQUED = 'C':  A := A * diag(C)
!>     EQUED = 'B':  A := diag(R) * A * diag(C).
!>
[in]LDA
!>          LDA is INTEGER
!>     The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]AF
!>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
!>     If FACT = 'F', then AF is an input argument and on entry
!>     contains the factors L and U from the factorization
!>     A = P*L*U as computed by DGETRF.  If EQUED .ne. 'N', then
!>     AF is the factored form of the equilibrated matrix A.
!>
!>     If FACT = 'N', then AF is an output argument and on exit
!>     returns the factors L and U from the factorization A = P*L*U
!>     of the original matrix A.
!>
!>     If FACT = 'E', then AF is an output argument and on exit
!>     returns the factors L and U from the factorization A = P*L*U
!>     of the equilibrated matrix A (see the description of A for
!>     the form of the equilibrated matrix).
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>     If FACT = 'F', then IPIV is an input argument and on entry
!>     contains the pivot indices from the factorization A = P*L*U
!>     as computed by DGETRF; row i of the matrix was interchanged
!>     with row IPIV(i).
!>
!>     If FACT = 'N', then IPIV is an output argument and on exit
!>     contains the pivot indices from the factorization A = P*L*U
!>     of the original matrix A.
!>
!>     If FACT = 'E', then IPIV is an output argument and on exit
!>     contains the pivot indices from the factorization A = P*L*U
!>     of the equilibrated matrix A.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>     Specifies the form of equilibration that was done.
!>       = 'N':  No equilibration (always true if FACT = 'N').
!>       = 'R':  Row equilibration, i.e., A has been premultiplied by
!>               diag(R).
!>       = 'C':  Column equilibration, i.e., A has been postmultiplied
!>               by diag(C).
!>       = 'B':  Both row and column equilibration, i.e., A has been
!>               replaced by diag(R) * A * diag(C).
!>     EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>     output argument.
!>
[in,out]R
!>          R is DOUBLE PRECISION array, dimension (N)
!>     The row scale factors for A.  If EQUED = 'R' or 'B', A is
!>     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
!>     is not accessed.  R is an input argument if FACT = 'F';
!>     otherwise, R is an output argument.  If FACT = 'F' and
!>     EQUED = 'R' or 'B', each element of R must be positive.
!>     If R is output, each element of R is a power of the radix.
!>     If R is input, each element of R should be a power of the radix
!>     to ensure a reliable solution and error estimates. Scaling by
!>     powers of the radix does not cause rounding errors unless the
!>     result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>     The column scale factors for A.  If EQUED = 'C' or 'B', A is
!>     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
!>     is not accessed.  C is an input argument if FACT = 'F';
!>     otherwise, C is an output argument.  If FACT = 'F' and
!>     EQUED = 'C' or 'B', each element of C must be positive.
!>     If C is output, each element of C is a power of the radix.
!>     If C is input, each element of C should be a power of the radix
!>     to ensure a reliable solution and error estimates. Scaling by
!>     powers of the radix does not cause rounding errors unless the
!>     result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>     On entry, the N-by-NRHS right hand side matrix B.
!>     On exit,
!>     if EQUED = 'N', B is not modified;
!>     if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
!>        diag(R)*B;
!>     if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
!>        overwritten by diag(C)*B.
!>
[in]LDB
!>          LDB is INTEGER
!>     The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
!>     If INFO = 0, the N-by-NRHS solution matrix X to the original
!>     system of equations.  Note that A and B are modified on exit
!>     if EQUED .ne. 'N', and the solution to the equilibrated system is
!>     inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
!>     inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
!>
[in]LDX
!>          LDX is INTEGER
!>     The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[out]RPVGRW
!>          RPVGRW is DOUBLE PRECISION
!>     Reciprocal pivot growth.  On exit, this contains the reciprocal
!>     pivot growth factor norm(A)/norm(U). The 
!>     norm is used.  If this is much less than 1, then the stability of
!>     the LU factorization of the (equilibrated) matrix A could be poor.
!>     This also means that the solution X, estimated condition numbers,
!>     and error bounds could be unreliable. If factorization fails with
!>     0<INFO<=N, then this contains the reciprocal pivot growth factor
!>     for the leading INFO columns of A.  In DGESVX, this quantity is
!>     returned in WORK(1).
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>     Componentwise relative backward error.  This is the
!>     componentwise relative backward error of each solution vector X(j)
!>     (i.e., the smallest relative change in any element of A or B that
!>     makes X(j) an exact solution).
!>
[in]N_ERR_BNDS
!>          N_ERR_BNDS is INTEGER
!>     Number of error bounds to return for each right hand side
!>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
!>     ERR_BNDS_COMP below.
!>
[out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * dlamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * dlamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * dlamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * dlamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.0D+0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the extra-precise refinement algorithm.
!>              (other values are reserved for future use)
!>
!>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
!>            computations allowed for refinement.
!>         Default: 10
!>         Aggressive: Set to 100 to permit convergence using approximate
!>                     factorizations or factorizations other than LU. If
!>                     the factorization uses a technique other than
!>                     Gaussian elimination, the guarantees in
!>                     err_bnds_norm and err_bnds_comp may no longer be
!>                     trustworthy.
!>
!>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
!>            will attempt to find a solution with small componentwise
!>            relative error in the double-precision algorithm.  Positive
!>            is true, 0.0 is false.
!>         Default: 1.0 (attempt componentwise convergence)
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (4*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit. The solution to every right-hand side is
!>         guaranteed.
!>       < 0:  If INFO = -i, the i-th argument had an illegal value
!>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
!>         has been completed, but the factor U is exactly singular, so
!>         the solution and error bounds could not be computed. RCOND = 0
!>         is returned.
!>       = N+J: The solution corresponding to the Jth right-hand side is
!>         not guaranteed. The solutions corresponding to other right-
!>         hand sides K with K > J may not be guaranteed as well, but
!>         only the first such right-hand side is reported. If a small
!>         componentwise error is not requested (PARAMS(3) = 0.0) then
!>         the Jth right-hand side is the first with a normwise error
!>         bound that is not guaranteed (the smallest J such
!>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
!>         the Jth right-hand side is the first with either a normwise or
!>         componentwise error bound that is not guaranteed (the smallest
!>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
!>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
!>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
!>         about all of the right-hand sides check ERR_BNDS_NORM or
!>         ERR_BNDS_COMP.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 535 of file dgesvxx.f.

540*
541* -- LAPACK driver routine --
542* -- LAPACK is a software package provided by Univ. of Tennessee, --
543* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
544*
545* .. Scalar Arguments ..
546 CHARACTER EQUED, FACT, TRANS
547 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
548 $ N_ERR_BNDS
549 DOUBLE PRECISION RCOND, RPVGRW
550* ..
551* .. Array Arguments ..
552 INTEGER IPIV( * ), IWORK( * )
553 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
554 $ X( LDX , * ),WORK( * )
555 DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
556 $ ERR_BNDS_NORM( NRHS, * ),
557 $ ERR_BNDS_COMP( NRHS, * )
558* ..
559*
560* =====================================================================
561*
562* .. Parameters ..
563 DOUBLE PRECISION ZERO, ONE
564 parameter( zero = 0.0d+0, one = 1.0d+0 )
565 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
566 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
567 INTEGER CMP_ERR_I, PIV_GROWTH_I
568 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
569 $ berr_i = 3 )
570 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
571 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
572 $ piv_growth_i = 9 )
573* ..
574* .. Local Scalars ..
575 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
576 INTEGER INFEQU, J
577 DOUBLE PRECISION AMAX, BIGNUM, COLCND, RCMAX, RCMIN, ROWCND,
578 $ SMLNUM
579* ..
580* .. External Functions ..
581 EXTERNAL lsame, dlamch, dla_gerpvgrw
582 LOGICAL LSAME
583 DOUBLE PRECISION DLAMCH, DLA_GERPVGRW
584* ..
585* .. External Subroutines ..
586 EXTERNAL dgeequb, dgetrf, dgetrs, dlacpy, dlaqge,
587 $ xerbla, dlascl2, dgerfsx
588* ..
589* .. Intrinsic Functions ..
590 INTRINSIC max, min
591* ..
592* .. Executable Statements ..
593*
594 info = 0
595 nofact = lsame( fact, 'N' )
596 equil = lsame( fact, 'E' )
597 notran = lsame( trans, 'N' )
598 smlnum = dlamch( 'Safe minimum' )
599 bignum = one / smlnum
600 IF( nofact .OR. equil ) THEN
601 equed = 'N'
602 rowequ = .false.
603 colequ = .false.
604 ELSE
605 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
606 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
607 END IF
608*
609* Default is failure. If an input parameter is wrong or
610* factorization fails, make everything look horrible. Only the
611* pivot growth is set here, the rest is initialized in DGERFSX.
612*
613 rpvgrw = zero
614*
615* Test the input parameters. PARAMS is not tested until DGERFSX.
616*
617 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
618 $ lsame( fact, 'F' ) ) THEN
619 info = -1
620 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
621 $ lsame( trans, 'C' ) ) THEN
622 info = -2
623 ELSE IF( n.LT.0 ) THEN
624 info = -3
625 ELSE IF( nrhs.LT.0 ) THEN
626 info = -4
627 ELSE IF( lda.LT.max( 1, n ) ) THEN
628 info = -6
629 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
630 info = -8
631 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
632 $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
633 info = -10
634 ELSE
635 IF( rowequ ) THEN
636 rcmin = bignum
637 rcmax = zero
638 DO 10 j = 1, n
639 rcmin = min( rcmin, r( j ) )
640 rcmax = max( rcmax, r( j ) )
641 10 CONTINUE
642 IF( rcmin.LE.zero ) THEN
643 info = -11
644 ELSE IF( n.GT.0 ) THEN
645 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
646 ELSE
647 rowcnd = one
648 END IF
649 END IF
650 IF( colequ .AND. info.EQ.0 ) THEN
651 rcmin = bignum
652 rcmax = zero
653 DO 20 j = 1, n
654 rcmin = min( rcmin, c( j ) )
655 rcmax = max( rcmax, c( j ) )
656 20 CONTINUE
657 IF( rcmin.LE.zero ) THEN
658 info = -12
659 ELSE IF( n.GT.0 ) THEN
660 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
661 ELSE
662 colcnd = one
663 END IF
664 END IF
665 IF( info.EQ.0 ) THEN
666 IF( ldb.LT.max( 1, n ) ) THEN
667 info = -14
668 ELSE IF( ldx.LT.max( 1, n ) ) THEN
669 info = -16
670 END IF
671 END IF
672 END IF
673*
674 IF( info.NE.0 ) THEN
675 CALL xerbla( 'DGESVXX', -info )
676 RETURN
677 END IF
678*
679 IF( equil ) THEN
680*
681* Compute row and column scalings to equilibrate the matrix A.
682*
683 CALL dgeequb( n, n, a, lda, r, c, rowcnd, colcnd, amax,
684 $ infequ )
685 IF( infequ.EQ.0 ) THEN
686*
687* Equilibrate the matrix.
688*
689 CALL dlaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
690 $ equed )
691 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
692 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
693 END IF
694*
695* If the scaling factors are not applied, set them to 1.0.
696*
697 IF ( .NOT.rowequ ) THEN
698 DO j = 1, n
699 r( j ) = 1.0d+0
700 END DO
701 END IF
702 IF ( .NOT.colequ ) THEN
703 DO j = 1, n
704 c( j ) = 1.0d+0
705 END DO
706 END IF
707 END IF
708*
709* Scale the right-hand side.
710*
711 IF( notran ) THEN
712 IF( rowequ ) CALL dlascl2( n, nrhs, r, b, ldb )
713 ELSE
714 IF( colequ ) CALL dlascl2( n, nrhs, c, b, ldb )
715 END IF
716*
717 IF( nofact .OR. equil ) THEN
718*
719* Compute the LU factorization of A.
720*
721 CALL dlacpy( 'Full', n, n, a, lda, af, ldaf )
722 CALL dgetrf( n, n, af, ldaf, ipiv, info )
723*
724* Return if INFO is non-zero.
725*
726 IF( info.GT.0 ) THEN
727*
728* Pivot in column INFO is exactly 0
729* Compute the reciprocal pivot growth factor of the
730* leading rank-deficient INFO columns of A.
731*
732 rpvgrw = dla_gerpvgrw( n, info, a, lda, af, ldaf )
733 RETURN
734 END IF
735 END IF
736*
737* Compute the reciprocal pivot growth factor RPVGRW.
738*
739 rpvgrw = dla_gerpvgrw( n, n, a, lda, af, ldaf )
740*
741* Compute the solution matrix X.
742*
743 CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
744 CALL dgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
745*
746* Use iterative refinement to improve the computed solution and
747* compute error bounds and backward error estimates for it.
748*
749 CALL dgerfsx( trans, equed, n, nrhs, a, lda, af, ldaf,
750 $ ipiv, r, c, b, ldb, x, ldx, rcond, berr,
751 $ n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params,
752 $ work, iwork, info )
753*
754* Scale solutions.
755*
756 IF ( colequ .AND. notran ) THEN
757 CALL dlascl2 ( n, nrhs, c, x, ldx )
758 ELSE IF ( rowequ .AND. .NOT.notran ) THEN
759 CALL dlascl2 ( n, nrhs, r, x, ldx )
760 END IF
761*
762 RETURN
763*
764* End of DGESVXX
765
dla_gerpvgrw
double precision function dla_gerpvgrw(n, ncols, a, lda, af, ldaf)
DLA_GERPVGRW
Definition dla_gerpvgrw.f:100
dgerfsx
subroutine dgerfsx(trans, equed, n, nrhs, a, lda, af, ldaf, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DGERFSX
Definition dgerfsx.f:414
dgeequb
subroutine dgeequb(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
DGEEQUB
Definition dgeequb.f:146
dlascl2
subroutine dlascl2(m, n, d, x, ldx)
DLASCL2 performs diagonal scaling on a vector.
Definition dlascl2.f:90

◆ dgetsls()

subroutine dgetsls ( character trans,
integer m,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) work,
integer lwork,
integer info )

DGETSLS

Purpose:
!>
!> DGETSLS solves overdetermined or underdetermined real linear systems
!> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
!> factorization of A.  It is assumed that A has full rank.
!>
!>
!>
!> The following options are provided:
!>
!> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A*X ||.
!>
!> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
!>    an underdetermined system A * X = B.
!>
!> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
!>    an undetermined system A**T * X = B.
!>
!> 4. If TRANS = 'T' and m < n:  find the least squares solution of
!>    an overdetermined system, i.e., solve the least squares problem
!>                 minimize || B - A**T * X ||.
!>
!> 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.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': the linear system involves A;
!>          = 'T': the linear system involves A**T.
!>
[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 the 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 is overwritten by details of its QR or LQ
!>          factorization as returned by DGEQR or DGELQ.
!>
[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 matrix B of right hand side vectors, stored
!>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
!>          if TRANS = 'T'.
!>          On exit, if INFO = 0, B is overwritten by the solution
!>          vectors, stored columnwise:
!>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
!>          squares solution vectors.
!>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
!>          minimum norm solution vectors;
!>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
!>          minimum norm solution vectors;
!>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
!>          least squares solution vectors.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= MAX(1,M,N).
!>
[out]WORK
!>          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
!>          or optimal, if query was assumed) LWORK.
!>          See LWORK for details.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If LWORK = -1 or -2, then a workspace query is assumed.
!>          If LWORK = -1, the routine calculates optimal size of WORK for the
!>          optimal performance and returns this value in WORK(1).
!>          If LWORK = -2, the routine calculates minimal size of WORK and 
!>          returns this value in WORK(1).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO =  i, the i-th diagonal element of the
!>                triangular factor of A is zero, so that A does not have
!>                full rank; the least squares solution could not be
!>                computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 160 of file dgetsls.f.

162*
163* -- LAPACK driver routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments ..
168 CHARACTER TRANS
169 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
170* ..
171* .. Array Arguments ..
172 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
173*
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 DOUBLE PRECISION ZERO, ONE
180 parameter( zero = 0.0d0, one = 1.0d0 )
181* ..
182* .. Local Scalars ..
183 LOGICAL LQUERY, TRAN
184 INTEGER I, IASCL, IBSCL, J, MAXMN, BROW,
185 $ SCLLEN, TSZO, TSZM, LWO, LWM, LW1, LW2,
186 $ WSIZEO, WSIZEM, INFO2
187 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ( 1 )
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 DOUBLE PRECISION DLAMCH, DLANGE
192 EXTERNAL lsame, dlabad, dlamch, dlange
193* ..
194* .. External Subroutines ..
195 EXTERNAL dgeqr, dgemqr, dlascl, dlaset,
196 $ dtrtrs, xerbla, dgelq, dgemlq
197* ..
198* .. Intrinsic Functions ..
199 INTRINSIC dble, max, min, int
200* ..
201* .. Executable Statements ..
202*
203* Test the input arguments.
204*
205 info = 0
206 maxmn = max( m, n )
207 tran = lsame( trans, 'T' )
208*
209 lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
210 IF( .NOT.( lsame( trans, 'N' ) .OR.
211 $ lsame( trans, 'T' ) ) ) THEN
212 info = -1
213 ELSE IF( m.LT.0 ) THEN
214 info = -2
215 ELSE IF( n.LT.0 ) THEN
216 info = -3
217 ELSE IF( nrhs.LT.0 ) THEN
218 info = -4
219 ELSE IF( lda.LT.max( 1, m ) ) THEN
220 info = -6
221 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
222 info = -8
223 END IF
224*
225 IF( info.EQ.0 ) THEN
226*
227* Determine the optimum and minimum LWORK
228*
229 IF( m.GE.n ) THEN
230 CALL dgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
231 tszo = int( tq( 1 ) )
232 lwo = int( workq( 1 ) )
233 CALL dgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
234 $ tszo, b, ldb, workq, -1, info2 )
235 lwo = max( lwo, int( workq( 1 ) ) )
236 CALL dgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
237 tszm = int( tq( 1 ) )
238 lwm = int( workq( 1 ) )
239 CALL dgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
240 $ tszm, b, ldb, workq, -1, info2 )
241 lwm = max( lwm, int( workq( 1 ) ) )
242 wsizeo = tszo + lwo
243 wsizem = tszm + lwm
244 ELSE
245 CALL dgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
246 tszo = int( tq( 1 ) )
247 lwo = int( workq( 1 ) )
248 CALL dgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
249 $ tszo, b, ldb, workq, -1, info2 )
250 lwo = max( lwo, int( workq( 1 ) ) )
251 CALL dgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
252 tszm = int( tq( 1 ) )
253 lwm = int( workq( 1 ) )
254 CALL dgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
255 $ tszm, b, ldb, workq, -1, info2 )
256 lwm = max( lwm, int( workq( 1 ) ) )
257 wsizeo = tszo + lwo
258 wsizem = tszm + lwm
259 END IF
260*
261 IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
262 info = -10
263 END IF
264*
265 work( 1 ) = dble( wsizeo )
266*
267 END IF
268*
269 IF( info.NE.0 ) THEN
270 CALL xerbla( 'DGETSLS', -info )
271 RETURN
272 END IF
273 IF( lquery ) THEN
274 IF( lwork.EQ.-2 ) work( 1 ) = dble( wsizem )
275 RETURN
276 END IF
277 IF( lwork.LT.wsizeo ) THEN
278 lw1 = tszm
279 lw2 = lwm
280 ELSE
281 lw1 = tszo
282 lw2 = lwo
283 END IF
284*
285* Quick return if possible
286*
287 IF( min( m, n, nrhs ).EQ.0 ) THEN
288 CALL dlaset( 'FULL', max( m, n ), nrhs, zero, zero,
289 $ b, ldb )
290 RETURN
291 END IF
292*
293* Get machine parameters
294*
295 smlnum = dlamch( 'S' ) / dlamch( 'P' )
296 bignum = one / smlnum
297 CALL dlabad( smlnum, bignum )
298*
299* Scale A, B if max element outside range [SMLNUM,BIGNUM]
300*
301 anrm = dlange( 'M', m, n, a, lda, work )
302 iascl = 0
303 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
304*
305* Scale matrix norm up to SMLNUM
306*
307 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
308 iascl = 1
309 ELSE IF( anrm.GT.bignum ) THEN
310*
311* Scale matrix norm down to BIGNUM
312*
313 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
314 iascl = 2
315 ELSE IF( anrm.EQ.zero ) THEN
316*
317* Matrix all zero. Return zero solution.
318*
319 CALL dlaset( 'F', maxmn, nrhs, zero, zero, b, ldb )
320 GO TO 50
321 END IF
322*
323 brow = m
324 IF ( tran ) THEN
325 brow = n
326 END IF
327 bnrm = dlange( 'M', brow, nrhs, b, ldb, work )
328 ibscl = 0
329 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
330*
331* Scale matrix norm up to SMLNUM
332*
333 CALL dlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
334 $ info )
335 ibscl = 1
336 ELSE IF( bnrm.GT.bignum ) THEN
337*
338* Scale matrix norm down to BIGNUM
339*
340 CALL dlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
341 $ info )
342 ibscl = 2
343 END IF
344*
345 IF ( m.GE.n ) THEN
346*
347* compute QR factorization of A
348*
349 CALL dgeqr( m, n, a, lda, work( lw2+1 ), lw1,
350 $ work( 1 ), lw2, info )
351 IF ( .NOT.tran ) THEN
352*
353* Least-Squares Problem min || A * X - B ||
354*
355* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
356*
357 CALL dgemqr( 'L' , 'T', m, nrhs, n, a, lda,
358 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
359 $ info )
360*
361* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
362*
363 CALL dtrtrs( 'U', 'N', 'N', n, nrhs,
364 $ a, lda, b, ldb, info )
365 IF( info.GT.0 ) THEN
366 RETURN
367 END IF
368 scllen = n
369 ELSE
370*
371* Overdetermined system of equations A**T * X = B
372*
373* B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
374*
375 CALL dtrtrs( 'U', 'T', 'N', n, nrhs,
376 $ a, lda, b, ldb, info )
377*
378 IF( info.GT.0 ) THEN
379 RETURN
380 END IF
381*
382* B(N+1:M,1:NRHS) = ZERO
383*
384 DO 20 j = 1, nrhs
385 DO 10 i = n + 1, m
386 b( i, j ) = zero
387 10 CONTINUE
388 20 CONTINUE
389*
390* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
391*
392 CALL dgemqr( 'L', 'N', m, nrhs, n, a, lda,
393 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
394 $ info )
395*
396 scllen = m
397*
398 END IF
399*
400 ELSE
401*
402* Compute LQ factorization of A
403*
404 CALL dgelq( m, n, a, lda, work( lw2+1 ), lw1,
405 $ work( 1 ), lw2, info )
406*
407* workspace at least M, optimally M*NB.
408*
409 IF( .NOT.tran ) THEN
410*
411* underdetermined system of equations A * X = B
412*
413* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
414*
415 CALL dtrtrs( 'L', 'N', 'N', m, nrhs,
416 $ a, lda, b, ldb, info )
417*
418 IF( info.GT.0 ) THEN
419 RETURN
420 END IF
421*
422* B(M+1:N,1:NRHS) = 0
423*
424 DO 40 j = 1, nrhs
425 DO 30 i = m + 1, n
426 b( i, j ) = zero
427 30 CONTINUE
428 40 CONTINUE
429*
430* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
431*
432 CALL dgemlq( 'L', 'T', n, nrhs, m, a, lda,
433 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
434 $ info )
435*
436* workspace at least NRHS, optimally NRHS*NB
437*
438 scllen = n
439*
440 ELSE
441*
442* overdetermined system min || A**T * X - B ||
443*
444* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
445*
446 CALL dgemlq( 'L', 'N', n, nrhs, m, a, lda,
447 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
448 $ info )
449*
450* workspace at least NRHS, optimally NRHS*NB
451*
452* B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
453*
454 CALL dtrtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
455 $ a, lda, b, ldb, info )
456*
457 IF( info.GT.0 ) THEN
458 RETURN
459 END IF
460*
461 scllen = m
462*
463 END IF
464*
465 END IF
466*
467* Undo scaling
468*
469 IF( iascl.EQ.1 ) THEN
470 CALL dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
471 $ info )
472 ELSE IF( iascl.EQ.2 ) THEN
473 CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
474 $ info )
475 END IF
476 IF( ibscl.EQ.1 ) THEN
477 CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
478 $ info )
479 ELSE IF( ibscl.EQ.2 ) THEN
480 CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
481 $ info )
482 END IF
483*
484 50 CONTINUE
485 work( 1 ) = dble( tszo + lwo )
486 RETURN
487*
488* End of DGETSLS
489*
dgelq
subroutine dgelq(m, n, a, lda, t, tsize, work, lwork, info)
DGELQ
Definition dgelq.f:172
dgemlq
subroutine dgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
DGEMLQ
Definition dgemlq.f:171
dgemqr
subroutine dgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
DGEMQR
Definition dgemqr.f:172
dgeqr
subroutine dgeqr(m, n, a, lda, t, tsize, work, lwork, info)
DGEQR
Definition dgeqr.f:174

◆ dsgesv()

subroutine dsgesv ( integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldx, * ) x,
integer ldx,
double precision, dimension( n, * ) work,
real, dimension( * ) swork,
integer iter,
integer info )

DSGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement)

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

Purpose:
!>
!> DSGESV computes the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!>
!> DSGESV first attempts to factorize the matrix in SINGLE PRECISION
!> and use this factorization within an iterative refinement procedure
!> to produce a solution with DOUBLE PRECISION normwise backward error
!> quality (see below). If the approach fails the method switches to a
!> DOUBLE PRECISION factorization and solve.
!>
!> The iterative refinement is not going to be a winning strategy if
!> the ratio SINGLE PRECISION performance over DOUBLE PRECISION
!> performance is too small. A reasonable strategy should take the
!> number of right-hand sides and the size of the matrix into account.
!> This might be done with a call to ILAENV in the future. Up to now, we
!> always try iterative refinement.
!>
!> The iterative refinement process is stopped if
!>     ITER > ITERMAX
!> or for all the RHS we have:
!>     RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
!> where
!>     o ITER is the number of the current iteration in the iterative
!>       refinement process
!>     o RNRM is the infinity-norm of the residual
!>     o XNRM is the infinity-norm of the solution
!>     o ANRM is the infinity-operator-norm of the matrix A
!>     o EPS is the machine epsilon returned by DLAMCH('Epsilon')
!> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
!> respectively.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array,
!>          dimension (LDA,N)
!>          On entry, the N-by-N coefficient matrix A.
!>          On exit, if iterative refinement has been successfully used
!>          (INFO = 0 and ITER >= 0, see description below), then A is
!>          unchanged, if double precision factorization has been used
!>          (INFO = 0 and ITER < 0, see description below), then the
!>          array A contains the factors L and U from the factorization
!>          A = P*L*U; the unit diagonal elements of L are not stored.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices that define the permutation matrix P;
!>          row i of the matrix was interchanged with row IPIV(i).
!>          Corresponds either to the single precision factorization
!>          (if INFO = 0 and ITER >= 0) or the double precision
!>          factorization (if INFO = 0 and ITER < 0).
!>
[in]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>          The N-by-NRHS right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
!>          If INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (N,NRHS)
!>          This array is used to hold the residual vectors.
!>
[out]SWORK
!>          SWORK is REAL array, dimension (N*(N+NRHS))
!>          This array is used to use the single precision matrix and the
!>          right-hand sides or solutions in single precision.
!>
[out]ITER
!>          ITER is INTEGER
!>          < 0: iterative refinement has failed, double precision
!>               factorization has been performed
!>               -1 : the routine fell back to full precision for
!>                    implementation- or machine-specific reasons
!>               -2 : narrowing the precision induced an overflow,
!>                    the routine fell back to full precision
!>               -3 : failure of SGETRF
!>               -31: stop the iterative refinement after the 30th
!>                    iterations
!>          > 0: iterative refinement has been successfully used.
!>               Returns the number of iterations
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, U(i,i) computed in DOUBLE PRECISION is
!>                exactly zero.  The factorization has been completed,
!>                but the factor U is exactly singular, so the solution
!>                could not be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 193 of file dsgesv.f.

195*
196* -- LAPACK driver routine --
197* -- LAPACK is a software package provided by Univ. of Tennessee, --
198* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199*
200* .. Scalar Arguments ..
201 INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
202* ..
203* .. Array Arguments ..
204 INTEGER IPIV( * )
205 REAL SWORK( * )
206 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
207 $ X( LDX, * )
208* ..
209*
210* =====================================================================
211*
212* .. Parameters ..
213 LOGICAL DOITREF
214 parameter( doitref = .true. )
215*
216 INTEGER ITERMAX
217 parameter( itermax = 30 )
218*
219 DOUBLE PRECISION BWDMAX
220 parameter( bwdmax = 1.0e+00 )
221*
222 DOUBLE PRECISION NEGONE, ONE
223 parameter( negone = -1.0d+0, one = 1.0d+0 )
224*
225* .. Local Scalars ..
226 INTEGER I, IITER, PTSA, PTSX
227 DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
228*
229* .. External Subroutines ..
230 EXTERNAL daxpy, dgemm, dlacpy, dlag2s, dgetrf, dgetrs,
231 $ sgetrf, sgetrs, slag2d, xerbla
232* ..
233* .. External Functions ..
234 INTEGER IDAMAX
235 DOUBLE PRECISION DLAMCH, DLANGE
236 EXTERNAL idamax, dlamch, dlange
237* ..
238* .. Intrinsic Functions ..
239 INTRINSIC abs, dble, max, sqrt
240* ..
241* .. Executable Statements ..
242*
243 info = 0
244 iter = 0
245*
246* Test the input parameters.
247*
248 IF( n.LT.0 ) THEN
249 info = -1
250 ELSE IF( nrhs.LT.0 ) THEN
251 info = -2
252 ELSE IF( lda.LT.max( 1, n ) ) THEN
253 info = -4
254 ELSE IF( ldb.LT.max( 1, n ) ) THEN
255 info = -7
256 ELSE IF( ldx.LT.max( 1, n ) ) THEN
257 info = -9
258 END IF
259 IF( info.NE.0 ) THEN
260 CALL xerbla( 'DSGESV', -info )
261 RETURN
262 END IF
263*
264* Quick return if (N.EQ.0).
265*
266 IF( n.EQ.0 )
267 $ RETURN
268*
269* Skip single precision iterative refinement if a priori slower
270* than double precision factorization.
271*
272 IF( .NOT.doitref ) THEN
273 iter = -1
274 GO TO 40
275 END IF
276*
277* Compute some constants.
278*
279 anrm = dlange( 'I', n, n, a, lda, work )
280 eps = dlamch( 'Epsilon' )
281 cte = anrm*eps*sqrt( dble( n ) )*bwdmax
282*
283* Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
284*
285 ptsa = 1
286 ptsx = ptsa + n*n
287*
288* Convert B from double precision to single precision and store the
289* result in SX.
290*
291 CALL dlag2s( n, nrhs, b, ldb, swork( ptsx ), n, info )
292*
293 IF( info.NE.0 ) THEN
294 iter = -2
295 GO TO 40
296 END IF
297*
298* Convert A from double precision to single precision and store the
299* result in SA.
300*
301 CALL dlag2s( n, n, a, lda, swork( ptsa ), n, info )
302*
303 IF( info.NE.0 ) THEN
304 iter = -2
305 GO TO 40
306 END IF
307*
308* Compute the LU factorization of SA.
309*
310 CALL sgetrf( n, n, swork( ptsa ), n, ipiv, info )
311*
312 IF( info.NE.0 ) THEN
313 iter = -3
314 GO TO 40
315 END IF
316*
317* Solve the system SA*SX = SB.
318*
319 CALL sgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
320 $ swork( ptsx ), n, info )
321*
322* Convert SX back to double precision
323*
324 CALL slag2d( n, nrhs, swork( ptsx ), n, x, ldx, info )
325*
326* Compute R = B - AX (R is WORK).
327*
328 CALL dlacpy( 'All', n, nrhs, b, ldb, work, n )
329*
330 CALL dgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone, a,
331 $ lda, x, ldx, one, work, n )
332*
333* Check whether the NRHS normwise backward errors satisfy the
334* stopping criterion. If yes, set ITER=0 and return.
335*
336 DO i = 1, nrhs
337 xnrm = abs( x( idamax( n, x( 1, i ), 1 ), i ) )
338 rnrm = abs( work( idamax( n, work( 1, i ), 1 ), i ) )
339 IF( rnrm.GT.xnrm*cte )
340 $ GO TO 10
341 END DO
342*
343* If we are here, the NRHS normwise backward errors satisfy the
344* stopping criterion. We are good to exit.
345*
346 iter = 0
347 RETURN
348*
349 10 CONTINUE
350*
351 DO 30 iiter = 1, itermax
352*
353* Convert R (in WORK) from double precision to single precision
354* and store the result in SX.
355*
356 CALL dlag2s( n, nrhs, work, n, swork( ptsx ), n, info )
357*
358 IF( info.NE.0 ) THEN
359 iter = -2
360 GO TO 40
361 END IF
362*
363* Solve the system SA*SX = SR.
364*
365 CALL sgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
366 $ swork( ptsx ), n, info )
367*
368* Convert SX back to double precision and update the current
369* iterate.
370*
371 CALL slag2d( n, nrhs, swork( ptsx ), n, work, n, info )
372*
373 DO i = 1, nrhs
374 CALL daxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
375 END DO
376*
377* Compute R = B - AX (R is WORK).
378*
379 CALL dlacpy( 'All', n, nrhs, b, ldb, work, n )
380*
381 CALL dgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone,
382 $ a, lda, x, ldx, one, work, n )
383*
384* Check whether the NRHS normwise backward errors satisfy the
385* stopping criterion. If yes, set ITER=IITER>0 and return.
386*
387 DO i = 1, nrhs
388 xnrm = abs( x( idamax( n, x( 1, i ), 1 ), i ) )
389 rnrm = abs( work( idamax( n, work( 1, i ), 1 ), i ) )
390 IF( rnrm.GT.xnrm*cte )
391 $ GO TO 20
392 END DO
393*
394* If we are here, the NRHS normwise backward errors satisfy the
395* stopping criterion, we are good to exit.
396*
397 iter = iiter
398*
399 RETURN
400*
401 20 CONTINUE
402*
403 30 CONTINUE
404*
405* If we are at this place of the code, this is because we have
406* performed ITER=ITERMAX iterations and never satisfied the
407* stopping criterion, set up the ITER flag accordingly and follow up
408* on double precision routine.
409*
410 iter = -itermax - 1
411*
412 40 CONTINUE
413*
414* Single-precision iterative refinement failed to converge to a
415* satisfactory solution, so we resort to double precision.
416*
417 CALL dgetrf( n, n, a, lda, ipiv, info )
418*
419 IF( info.NE.0 )
420 $ RETURN
421*
422 CALL dlacpy( 'All', n, nrhs, b, ldb, x, ldx )
423 CALL dgetrs( 'No transpose', n, nrhs, a, lda, ipiv, x, ldx,
424 $ info )
425*
426 RETURN
427*
428* End of DSGESV
429*
slag2d
subroutine slag2d(m, n, sa, ldsa, a, lda, info)
SLAG2D converts a single precision matrix to a double precision matrix.
Definition slag2d.f:104
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlag2s
subroutine dlag2s(m, n, a, lda, sa, ldsa, info)
DLAG2S converts a double precision matrix to a single precision matrix.
Definition dlag2s.f:108
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
sgetrf
subroutine sgetrf(m, n, a, lda, ipiv, info)
SGETRF
Definition sgetrf.f:108
sgetrs
subroutine sgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
SGETRS
Definition sgetrs.f:121