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

◆ 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 (MAX(1,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 344 of file sgesvx.f.

348*
349* -- LAPACK driver routine --
350* -- LAPACK is a software package provided by Univ. of Tennessee, --
351* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
352*
353* .. Scalar Arguments ..
354 CHARACTER EQUED, FACT, TRANS
355 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
356 REAL RCOND
357* ..
358* .. Array Arguments ..
359 INTEGER IPIV( * ), IWORK( * )
360 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
361 $ BERR( * ), C( * ), FERR( * ), R( * ),
362 $ WORK( * ), X( LDX, * )
363* ..
364*
365* =====================================================================
366*
367* .. Parameters ..
368 REAL ZERO, ONE
369 parameter( zero = 0.0e+0, one = 1.0e+0 )
370* ..
371* .. Local Scalars ..
372 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
373 CHARACTER NORM
374 INTEGER I, INFEQU, J
375 REAL AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
376 $ ROWCND, RPVGRW, SMLNUM
377* ..
378* .. External Functions ..
379 LOGICAL LSAME
380 REAL SLAMCH, SLANGE, SLANTR
381 EXTERNAL lsame, slamch, slange, slantr
382* ..
383* .. External Subroutines ..
384 EXTERNAL sgecon, sgeequ, sgerfs, sgetrf, sgetrs,
385 $ 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.
411 $ .NOT.equil .AND.
412 $ .NOT.lsame( fact, 'F' ) )
413 $ THEN
414 info = -1
415 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
416 $ lsame( trans, 'C' ) ) THEN
417 info = -2
418 ELSE IF( n.LT.0 ) THEN
419 info = -3
420 ELSE IF( nrhs.LT.0 ) THEN
421 info = -4
422 ELSE IF( lda.LT.max( 1, n ) ) THEN
423 info = -6
424 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
425 info = -8
426 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
427 $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
428 info = -10
429 ELSE
430 IF( rowequ ) THEN
431 rcmin = bignum
432 rcmax = zero
433 DO 10 j = 1, n
434 rcmin = min( rcmin, r( j ) )
435 rcmax = max( rcmax, r( j ) )
436 10 CONTINUE
437 IF( rcmin.LE.zero ) THEN
438 info = -11
439 ELSE IF( n.GT.0 ) THEN
440 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
441 ELSE
442 rowcnd = one
443 END IF
444 END IF
445 IF( colequ .AND. info.EQ.0 ) THEN
446 rcmin = bignum
447 rcmax = zero
448 DO 20 j = 1, n
449 rcmin = min( rcmin, c( j ) )
450 rcmax = max( rcmax, c( j ) )
451 20 CONTINUE
452 IF( rcmin.LE.zero ) THEN
453 info = -12
454 ELSE IF( n.GT.0 ) THEN
455 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
456 ELSE
457 colcnd = one
458 END IF
459 END IF
460 IF( info.EQ.0 ) THEN
461 IF( ldb.LT.max( 1, n ) ) THEN
462 info = -14
463 ELSE IF( ldx.LT.max( 1, n ) ) THEN
464 info = -16
465 END IF
466 END IF
467 END IF
468*
469 IF( info.NE.0 ) THEN
470 CALL xerbla( 'SGESVX', -info )
471 RETURN
472 END IF
473*
474 IF( equil ) THEN
475*
476* Compute row and column scalings to equilibrate the matrix A.
477*
478 CALL sgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax,
479 $ infequ )
480 IF( infequ.EQ.0 ) THEN
481*
482* Equilibrate the matrix.
483*
484 CALL slaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
485 $ equed )
486 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
487 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
488 END IF
489 END IF
490*
491* Scale the right hand side.
492*
493 IF( notran ) THEN
494 IF( rowequ ) THEN
495 DO 40 j = 1, nrhs
496 DO 30 i = 1, n
497 b( i, j ) = r( i )*b( i, j )
498 30 CONTINUE
499 40 CONTINUE
500 END IF
501 ELSE IF( colequ ) THEN
502 DO 60 j = 1, nrhs
503 DO 50 i = 1, n
504 b( i, j ) = c( i )*b( i, j )
505 50 CONTINUE
506 60 CONTINUE
507 END IF
508*
509 IF( nofact .OR. equil ) THEN
510*
511* Compute the LU factorization of A.
512*
513 CALL slacpy( 'Full', n, n, a, lda, af, ldaf )
514 CALL sgetrf( n, n, af, ldaf, ipiv, info )
515*
516* Return if INFO is non-zero.
517*
518 IF( info.GT.0 ) THEN
519*
520* Compute the reciprocal pivot growth factor of the
521* leading rank-deficient INFO columns of A.
522*
523 rpvgrw = slantr( 'M', 'U', 'N', info, info, af, ldaf,
524 $ work )
525 IF( rpvgrw.EQ.zero ) THEN
526 rpvgrw = one
527 ELSE
528 rpvgrw = slange( 'M', n, info, a, lda, work ) / rpvgrw
529 END IF
530 work( 1 ) = rpvgrw
531 rcond = zero
532 RETURN
533 END IF
534 END IF
535*
536* Compute the norm of the matrix A and the
537* reciprocal pivot growth factor RPVGRW.
538*
539 IF( notran ) THEN
540 norm = '1'
541 ELSE
542 norm = 'I'
543 END IF
544 anorm = slange( norm, n, n, a, lda, work )
545 rpvgrw = slantr( 'M', 'U', 'N', n, n, af, ldaf, work )
546 IF( rpvgrw.EQ.zero ) THEN
547 rpvgrw = one
548 ELSE
549 rpvgrw = slange( 'M', n, n, a, lda, work ) / rpvgrw
550 END IF
551*
552* Compute the reciprocal of the condition number of A.
553*
554 CALL sgecon( norm, n, af, ldaf, anorm, rcond, work, iwork,
555 $ info )
556*
557* Compute the solution matrix X.
558*
559 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
560 CALL sgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
561*
562* Use iterative refinement to improve the computed solution and
563* compute error bounds and backward error estimates for it.
564*
565 CALL sgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
566 $ ldx, ferr, berr, work, iwork, info )
567*
568* Transform the solution matrix X to a solution of the original
569* system.
570*
571 IF( notran ) THEN
572 IF( colequ ) THEN
573 DO 80 j = 1, nrhs
574 DO 70 i = 1, n
575 x( i, j ) = c( i )*x( i, j )
576 70 CONTINUE
577 80 CONTINUE
578 DO 90 j = 1, nrhs
579 ferr( j ) = ferr( j ) / colcnd
580 90 CONTINUE
581 END IF
582 ELSE IF( rowequ ) THEN
583 DO 110 j = 1, nrhs
584 DO 100 i = 1, n
585 x( i, j ) = r( i )*x( i, j )
586 100 CONTINUE
587 110 CONTINUE
588 DO 120 j = 1, nrhs
589 ferr( j ) = ferr( j ) / rowcnd
590 120 CONTINUE
591 END IF
592*
593* Set INFO = N+1 if the matrix is singular to working precision.
594*
595 IF( rcond.LT.slamch( 'Epsilon' ) )
596 $ info = n + 1
597*
598 work( 1 ) = rpvgrw
599 RETURN
600*
601* End of SGESVX
602*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgecon
subroutine sgecon(norm, n, a, lda, anorm, rcond, work, iwork, info)
SGECON
Definition sgecon.f:130
sgeequ
subroutine sgeequ(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
SGEEQU
Definition sgeequ.f:137
sgerfs
subroutine sgerfs(trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SGERFS
Definition sgerfs.f:184
sgetrf
subroutine sgetrf(m, n, a, lda, ipiv, info)
SGETRF
Definition sgetrf.f:106
sgetrs
subroutine sgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
SGETRS
Definition sgetrs.f:119
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:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:112
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:140
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:140
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: