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

Functions

subroutine sgelsx (m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, info)
  SGELSX solves overdetermined or underdetermined systems for GE matrices
subroutine sgels (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
  SGELS solves overdetermined or underdetermined systems for GE matrices
subroutine sgelsd (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info)
  SGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
subroutine sgelss (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, info)
  SGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine sgelsy (m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)
  SGELSY solves overdetermined or underdetermined systems for GE matrices
subroutine sgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
  SGESV computes the solution to system of linear equations A * X = B for GE matrices (simple driver)
subroutine sgesvx (fact, trans, n, nrhs, a, lda, af, ldaf, ipiv, equed, r, c, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
  SGESVX computes the solution to system of linear equations A * X = B for GE matrices
subroutine sgesvxx (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)
  SGESVXX computes the solution to system of linear equations A * X = B for GE matrices
subroutine sgetsls (trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
 SGETSLS

Detailed Description

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

Function Documentation

◆ sgels()

subroutine sgels ( character trans,
integer m,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) work,
integer lwork,
integer info )

SGELS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> SGELS 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 REAL 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 SGEQRF;
!>            if M <  N, A is overwritten by details of its LQ
!>                       factorization as returned by SGELQF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]B
!>          B is REAL 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 REAL 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 sgels.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 REAL A( LDA, * ), B( LDB, * ), WORK( * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO, ONE
200 parameter( zero = 0.0e0, one = 1.0e0 )
201* ..
202* .. Local Scalars ..
203 LOGICAL LQUERY, TPSD
204 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
205 REAL ANRM, BIGNUM, BNRM, SMLNUM
206* ..
207* .. Local Arrays ..
208 REAL RWORK( 1 )
209* ..
210* .. External Functions ..
211 LOGICAL LSAME
212 INTEGER ILAENV
213 REAL SLAMCH, SLANGE
214 EXTERNAL lsame, ilaenv, slamch, slange
215* ..
216* .. External Subroutines ..
217 EXTERNAL sgelqf, sgeqrf, slabad, slascl, slaset, sormlq,
218 $ sormqr, strtrs, xerbla
219* ..
220* .. Intrinsic Functions ..
221 INTRINSIC max, min, real
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.
243 $ .NOT.lquery ) 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, 'SGEQRF', ' ', m, n, -1, -1 )
257 IF( tpsd ) THEN
258 nb = max( nb, ilaenv( 1, 'SORMQR', 'LN', m, nrhs, n,
259 $ -1 ) )
260 ELSE
261 nb = max( nb, ilaenv( 1, 'SORMQR', 'LT', m, nrhs, n,
262 $ -1 ) )
263 END IF
264 ELSE
265 nb = ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
266 IF( tpsd ) THEN
267 nb = max( nb, ilaenv( 1, 'SORMLQ', 'LT', n, nrhs, m,
268 $ -1 ) )
269 ELSE
270 nb = max( nb, ilaenv( 1, 'SORMLQ', '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 ) = real( wsize )
277*
278 END IF
279*
280 IF( info.NE.0 ) THEN
281 CALL xerbla( 'SGELS ', -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 slaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
291 RETURN
292 END IF
293*
294* Get machine parameters
295*
296 smlnum = slamch( 'S' ) / slamch( 'P' )
297 bignum = one / smlnum
298 CALL slabad( smlnum, bignum )
299*
300* Scale A, B if max element outside range [SMLNUM,BIGNUM]
301*
302 anrm = slange( '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 slascl( '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 slascl( '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 slaset( '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 = slange( '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 slascl( '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 slascl( '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 sgeqrf( 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 sormqr( '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 strtrs( '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 strtrs( '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 sormqr( '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 sgelqf( 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 strtrs( '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 sormlq( '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 sormlq( '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 strtrs( '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 slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
481 $ info )
482 ELSE IF( iascl.EQ.2 ) THEN
483 CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
484 $ info )
485 END IF
486 IF( ibscl.EQ.1 ) THEN
487 CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
488 $ info )
489 ELSE IF( ibscl.EQ.2 ) THEN
490 CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
491 $ info )
492 END IF
493*
494 50 CONTINUE
495 work( 1 ) = real( wsize )
496*
497 RETURN
498*
499* End of SGELS
500*
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:143
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
sgelqf
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:143
strtrs
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
sormlq
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
Definition sormlq.f:168
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ sgelsd()

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

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

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

Purpose:
!>
!> SGELSD 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 REAL 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 REAL 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 REAL 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 REAL
!>          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 REAL 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 array WORK and the
!>          minimum size of the array IWORK, and returns these values as
!>          the first entries of the WORK and IWORK arrays, 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 208 of file sgelsd.f.

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

◆ sgelss()

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

SGELSS solves overdetermined or underdetermined systems for GE matrices

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

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

◆ sgelsx()

subroutine sgelsx ( integer m,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
real rcond,
integer rank,
real, dimension( * ) work,
integer info )

SGELSX solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine SGELSY.
!>
!> SGELSX 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 REAL 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 REAL 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 REAL
!>          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 REAL 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 sgelsx.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 REAL RCOND
186* ..
187* .. Array Arguments ..
188 INTEGER JPVT( * )
189 REAL A( LDA, * ), B( LDB, * ), WORK( * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 INTEGER IMAX, IMIN
196 parameter( imax = 1, imin = 2 )
197 REAL ZERO, ONE, DONE, NTDONE
198 parameter( zero = 0.0e0, one = 1.0e0, done = zero,
199 $ ntdone = one )
200* ..
201* .. Local Scalars ..
202 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
203 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
204 $ SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2
205* ..
206* .. External Functions ..
207 REAL SLAMCH, SLANGE
208 EXTERNAL slamch, slange
209* ..
210* .. External Subroutines ..
211 EXTERNAL sgeqpf, slabad, slaic1, slascl, slaset, slatzm,
212 $ sorm2r, strsm, stzrqf, 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( 'SGELSX', -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 = slamch( 'S' ) / slamch( 'P' )
253 bignum = one / smlnum
254 CALL slabad( smlnum, bignum )
255*
256* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
257*
258 anrm = slange( '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 slascl( '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 slascl( '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 slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
277 rank = 0
278 GO TO 100
279 END IF
280*
281 bnrm = slange( '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 slascl( '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 slascl( '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 sgeqpf( 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 slaset( '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 slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
323 $ a( i, i ), sminpr, s1, c1 )
324 CALL slaic1( 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 stzrqf( 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 sorm2r( '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 strsm( '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 slatzm( '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 slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
413 CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
414 $ info )
415 ELSE IF( iascl.EQ.2 ) THEN
416 CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
417 CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
418 $ info )
419 END IF
420 IF( ibscl.EQ.1 ) THEN
421 CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
422 ELSE IF( ibscl.EQ.2 ) THEN
423 CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
424 END IF
425*
426 100 CONTINUE
427*
428 RETURN
429*
430* End of SGELSX
431*
sgeqpf
subroutine sgeqpf(m, n, a, lda, jpvt, tau, work, info)
SGEQPF
Definition sgeqpf.f:142
slaic1
subroutine slaic1(job, j, x, sest, w, gamma, sestpr, s, c)
SLAIC1 applies one step of incremental condition estimation.
Definition slaic1.f:134
stzrqf
subroutine stzrqf(m, n, a, lda, tau, info)
STZRQF
Definition stzrqf.f:138
slatzm
subroutine slatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
SLATZM
Definition slatzm.f:151
sorm2r
subroutine sorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition sorm2r.f:159
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181

◆ sgelsy()

subroutine sgelsy ( integer m,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
real rcond,
integer rank,
real, dimension( * ) work,
integer lwork,
integer info )

SGELSY solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> SGELSY 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 REAL 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 REAL 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 REAL
!>          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 REAL 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 SGEQP3, STZRZF, STZRQF, SORMQR,
!>          and SORMRZ.
!>
!>          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 sgelsy.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 REAL RCOND
212* ..
213* .. Array Arguments ..
214 INTEGER JPVT( * )
215 REAL A( LDA, * ), B( LDB, * ), WORK( * )
216* ..
217*
218* =====================================================================
219*
220* .. Parameters ..
221 INTEGER IMAX, IMIN
222 parameter( imax = 1, imin = 2 )
223 REAL ZERO, ONE
224 parameter( zero = 0.0e+0, one = 1.0e+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 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
231 $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
232* ..
233* .. External Functions ..
234 INTEGER ILAENV
235 REAL SLAMCH, SLANGE
236 EXTERNAL ilaenv, slamch, slange
237* ..
238* .. External Subroutines ..
239 EXTERNAL scopy, sgeqp3, slabad, slaic1, slascl, slaset,
240 $ sormqr, sormrz, strsm, stzrzf, 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, 'SGEQRF', ' ', m, n, -1, -1 )
275 nb2 = ilaenv( 1, 'SGERQF', ' ', m, n, -1, -1 )
276 nb3 = ilaenv( 1, 'SORMQR', ' ', m, n, nrhs, -1 )
277 nb4 = ilaenv( 1, 'SORMRQ', ' ', 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( 'SGELSY', -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 = slamch( 'S' ) / slamch( 'P' )
307 bignum = one / smlnum
308 CALL slabad( smlnum, bignum )
309*
310* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
311*
312 anrm = slange( '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 slascl( '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 slascl( '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 slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
331 rank = 0
332 GO TO 70
333 END IF
334*
335 bnrm = slange( '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 slascl( '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 slascl( '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 sgeqp3( 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 slaset( '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 slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
379 $ a( i, i ), sminpr, s1, c1 )
380 CALL slaic1( 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 stzrzf( 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 sormqr( '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 strsm( '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 sormrz( '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 scopy( 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 slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
456 CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
457 $ info )
458 ELSE IF( iascl.EQ.2 ) THEN
459 CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
460 CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
461 $ info )
462 END IF
463 IF( ibscl.EQ.1 ) THEN
464 CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
465 ELSE IF( ibscl.EQ.2 ) THEN
466 CALL slascl( '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 SGELSY
475*
sgeqp3
subroutine sgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
SGEQP3
Definition sgeqp3.f:151
stzrzf
subroutine stzrzf(m, n, a, lda, tau, work, lwork, info)
STZRZF
Definition stzrzf.f:151
sormrz
subroutine sormrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
SORMRZ
Definition sormrz.f:187

◆ sgesv()

subroutine sgesv ( integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
real, dimension( ldb, * ) b,
integer ldb,
integer info )

SGESV computes the solution to system of linear equations A * X = B for GE matrices (simple driver)

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

Purpose:
!>
!> SGESV 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 REAL 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 REAL 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 sgesv.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 REAL A( LDA, * ), B( LDB, * )
133* ..
134*
135* =====================================================================
136*
137* .. External Subroutines ..
138 EXTERNAL sgetrf, sgetrs, 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( 'SGESV ', -info )
159 RETURN
160 END IF
161*
162* Compute the LU factorization of A.
163*
164 CALL sgetrf( 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 sgetrs( 'No transpose', n, nrhs, a, lda, ipiv, b, ldb,
170 $ info )
171 END IF
172 RETURN
173*
174* End of SGESV
175*
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

◆ sgesvx()

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

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

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

Purpose:
!>
!> SGESVX 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 REAL 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 REAL 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 SGETRF.  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 SGETRF; 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL
!>          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 REAL 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 REAL 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 REAL 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 sgesvx.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 REAL RCOND
358* ..
359* .. Array Arguments ..
360 INTEGER IPIV( * ), IWORK( * )
361 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
362 $ BERR( * ), C( * ), FERR( * ), R( * ),
363 $ WORK( * ), X( LDX, * )
364* ..
365*
366* =====================================================================
367*
368* .. Parameters ..
369 REAL ZERO, ONE
370 parameter( zero = 0.0e+0, one = 1.0e+0 )
371* ..
372* .. Local Scalars ..
373 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
374 CHARACTER NORM
375 INTEGER I, INFEQU, J
376 REAL AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
377 $ ROWCND, RPVGRW, SMLNUM
378* ..
379* .. External Functions ..
380 LOGICAL LSAME
381 REAL SLAMCH, SLANGE, SLANTR
382 EXTERNAL lsame, slamch, slange, slantr
383* ..
384* .. External Subroutines ..
385 EXTERNAL sgecon, sgeequ, sgerfs, sgetrf, sgetrs, slacpy,
386 $ slaqge, 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 = slamch( '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( 'SGESVX', -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 sgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
477 IF( infequ.EQ.0 ) THEN
478*
479* Equilibrate the matrix.
480*
481 CALL slaqge( 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 slacpy( 'Full', n, n, a, lda, af, ldaf )
511 CALL sgetrf( 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 = slantr( 'M', 'U', 'N', info, info, af, ldaf,
521 $ work )
522 IF( rpvgrw.EQ.zero ) THEN
523 rpvgrw = one
524 ELSE
525 rpvgrw = slange( '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 = slange( norm, n, n, a, lda, work )
542 rpvgrw = slantr( 'M', 'U', 'N', n, n, af, ldaf, work )
543 IF( rpvgrw.EQ.zero ) THEN
544 rpvgrw = one
545 ELSE
546 rpvgrw = slange( 'M', n, n, a, lda, work ) / rpvgrw
547 END IF
548*
549* Compute the reciprocal of the condition number of A.
550*
551 CALL sgecon( norm, n, af, ldaf, anorm, rcond, work, iwork, info )
552*
553* Compute the solution matrix X.
554*
555 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
556 CALL sgetrs( 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 sgerfs( 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* Set INFO = N+1 if the matrix is singular to working precision.
590*
591 IF( rcond.LT.slamch( 'Epsilon' ) )
592 $ info = n + 1
593*
594 work( 1 ) = rpvgrw
595 RETURN
596*
597* End of SGESVX
598*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
slaqge
subroutine slaqge(m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
SLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
Definition slaqge.f:142
sgecon
subroutine sgecon(norm, n, a, lda, anorm, rcond, work, iwork, info)
SGECON
Definition sgecon.f:124
sgeequ
subroutine sgeequ(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
SGEEQU
Definition sgeequ.f:139
sgerfs
subroutine sgerfs(trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SGERFS
Definition sgerfs.f:185
slantr
real function slantr(norm, uplo, diag, m, n, a, lda, work)
SLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slantr.f:141

◆ sgesvxx()

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

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

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

Purpose:
!>
!>    SGESVXX 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.
!>
!>    If requested, both normwise and maximum componentwise error bounds
!>    are returned. SGESVXX 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.
!>
!>    SGESVXX 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
!>    SGESVXX 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 SGESVXX would itself produce.
!>
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 (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 REAL 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 REAL 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 SGETRF.  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 SGETRF; 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL
!>     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 REAL
!>     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 SGESVX, this quantity is
!>     returned in WORK(1).
!>
[out]BERR
!>          BERR is REAL 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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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.0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the double-precision refinement algorithm,
!>                    possibly with doubled-single computations if the
!>                    compilation environment does not support DOUBLE
!>                    PRECISION.
!>              (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 REAL 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 538 of file sgesvxx.f.

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

◆ sgetsls()

subroutine sgetsls ( character trans,
integer m,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) work,
integer lwork,
integer info )

SGETSLS

Purpose:
!>
!> SGETSLS 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 REAL 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 SGEQR or SGELQ.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]B
!>          B is REAL 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) REAL 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 sgetsls.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 REAL A( LDA, * ), B( LDB, * ), WORK( * )
173*
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 REAL ZERO, ONE
180 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ( 1 )
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 REAL SLAMCH, SLANGE
192 EXTERNAL lsame, slabad, slamch, slange
193* ..
194* .. External Subroutines ..
195 EXTERNAL sgeqr, sgemqr, slascl, slaset,
196 $ strtrs, xerbla, sgelq, sgemlq
197* ..
198* .. Intrinsic Functions ..
199 INTRINSIC real, 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 sgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
231 tszo = int( tq( 1 ) )
232 lwo = int( workq( 1 ) )
233 CALL sgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
234 $ tszo, b, ldb, workq, -1, info2 )
235 lwo = max( lwo, int( workq( 1 ) ) )
236 CALL sgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
237 tszm = int( tq( 1 ) )
238 lwm = int( workq( 1 ) )
239 CALL sgemqr( '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 sgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
246 tszo = int( tq( 1 ) )
247 lwo = int( workq( 1 ) )
248 CALL sgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
249 $ tszo, b, ldb, workq, -1, info2 )
250 lwo = max( lwo, int( workq( 1 ) ) )
251 CALL sgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
252 tszm = int( tq( 1 ) )
253 lwm = int( workq( 1 ) )
254 CALL sgemlq( '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 ) = real( wsizeo )
266*
267 END IF
268*
269 IF( info.NE.0 ) THEN
270 CALL xerbla( 'SGETSLS', -info )
271 RETURN
272 END IF
273 IF( lquery ) THEN
274 IF( lwork.EQ.-2 ) work( 1 ) = real( 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 slaset( 'FULL', max( m, n ), nrhs, zero, zero,
289 $ b, ldb )
290 RETURN
291 END IF
292*
293* Get machine parameters
294*
295 smlnum = slamch( 'S' ) / slamch( 'P' )
296 bignum = one / smlnum
297 CALL slabad( smlnum, bignum )
298*
299* Scale A, B if max element outside range [SMLNUM,BIGNUM]
300*
301 anrm = slange( '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 slascl( '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 slascl( '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 slaset( '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 = slange( '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 slascl( '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 slascl( '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 sgeqr( 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 sgemqr( '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 strtrs( '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 strtrs( '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 sgemqr( '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 sgelq( 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 strtrs( '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 sgemlq( '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 sgemlq( '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 strtrs( '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 slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
471 $ info )
472 ELSE IF( iascl.EQ.2 ) THEN
473 CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
474 $ info )
475 END IF
476 IF( ibscl.EQ.1 ) THEN
477 CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
478 $ info )
479 ELSE IF( ibscl.EQ.2 ) THEN
480 CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
481 $ info )
482 END IF
483*
484 50 CONTINUE
485 work( 1 ) = real( tszo + lwo )
486 RETURN
487*
488* End of SGETSLS
489*
sgelq
subroutine sgelq(m, n, a, lda, t, tsize, work, lwork, info)
SGELQ
Definition sgelq.f:172
sgemlq
subroutine sgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
SGEMLQ
Definition sgemlq.f:170
sgemqr
subroutine sgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
SGEMQR
Definition sgemqr.f:172
sgeqr
subroutine sgeqr(m, n, a, lda, t, tsize, work, lwork, info)
SGEQR
Definition sgeqr.f:174