SUBROUTINE DSGESV(N,NRHS,A,LDA,IPIV,X,LDX,WORK,SWORK,ITER,INFO)
*
* -- LAPACK driver routine (version @(version)) --
* @(who)
* @(date)
*
* .. Scalar Arguments ..
INTEGER INFO,ITER,LDA,LDX,N,NRHS
DOUBLE PRECISION BWD
* ..
* .. Array Arguments ..
INTEGER IPIV(*)
REAL SWORK(*)
DOUBLE PRECISION A(LDA,*),WORK(*),X(LDX,*)
* ..
*
* Purpose
* =======
*
* DSGESV computes the solution to a real system of linear equations
* A * X = B,
* where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
*
* DSGESV first tries to factorize the matrix in SINGLE PRECISION and
* use this factorization within an iterative refinement procedure to
* produce a solution with DOUBLE PRECISION normwise backward error
* quality (see below). If the approach fails the method switch to a
* DOUBLE PRECISION factorization and solve.
*
* The iterative refinement process is stopped if
* ITER < ITERMAX = 30
* or
* BWD = RNRM /
* ( XNRM * ANRM * EPS * MIN(4.0E00,SQRT(FLOAT(N))/6E00) )
* < 1
* where
* o ITER is the number of iteration in the iterative refinement
* process
* o RNRM is the 2-norm of the residual
* o XNRM is the 2-norm of the solution
* o ANRM is the Frobenius-norm of the matrix A
* o EPS is the relative machine precision returned by DLAMCH
*
* Arguments
* =========
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. M >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrix B. NRHS >= 0.
*
* A (input or input/ouptut) DOUBLE PRECISION array,
* dimension (LDA,N)
* On entry, the M-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.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
*
* IPIV (output) INTEGER array, dimension (N)
* The pivot indices that define the permutation matrix P;
* row i of the matrix was interchanged with row IPIV(i).
*
* X (input/output) DOUBLE PRECISION array, dimension (LDX,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.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* WORK (workspace) DOUBLE PRECISION array, dimension (2*N*NRHS)
*
* SWORK (workspace) SINGLE PRECISION array, dimension (N*(N+NRHS))
*
* ITER (output) INTEGER
* < 0: iterative refinement has failed
* -1 : SPDRAT is too small, it is not worth working in
* SINGLE PRECISION
* -2 : overflow of an entry when moving from double to
* SINGLE PRECISION
* -3 : failure of SGETRF
* -4 : failure of SGETRS
* -31: stop the iterative refinement after the 30th
* iterations
* > 0: iterative refinement has been used successfully.
* Returns the number of iteration
*
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is
* exactly zero. The factorization has been completed,
* but the factor U is exactly singular, so the solution
* could not be computed.
*
* =========
*
* .. Parameters ..
DOUBLE PRECISION MONE,ONE
PARAMETER (MONE=-1.0D0,ONE=1.0D0)
*
* .. Local Scalars ..
INTEGER I,IITER,ITERMAX,DPTB,DPTR,OK,SPTSA,SPTSX
DOUBLE PRECISION ANRM,BWDCTE,EPS,RNRM,SPDRAT,XNRM
*
* .. External Subroutines ..
EXTERNAL DAXPY,DGEMM,DLACPY,DSLACONVERTGED2S,DSLACONVERTGES2D,
+ SGETRF,SGETRS,XERBLA
* ..
* .. External Functions ..
INTEGER ILAENV
DOUBLE PRECISION DLAMCH,DLANGE,DNRM2
EXTERNAL DLAMCH,DLANGE,DNRM2,ILAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX,SQRT
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
ITERMAX = 30
SPDRAT = 100.00E+00
*
BWD = 10.00E00
ITER = 0
INFO = 0
*
IF (N.LT.0) THEN
INFO = -1
ELSE IF (NRHS.LT.0) THEN
INFO = -2
ELSE IF (LDA.LT.MAX(1,N)) THEN
INFO = -4
ELSE IF (LDX.LT.MAX(1,N)) THEN
INFO = -7
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSGESV ',-INFO)
RETURN
END IF
*
* Quick return if N .EQ. 0
*
IF (N.EQ.0) RETURN
*
* The iterative refinement is not going to be a winning strategy if
* SPDRAT, the speed ratio low_prec_perf/high_prec_perf is too small.
* (will need to take the number of right-hand sides into account as
* well)
*
IF (SPDRAT.LT.1.5) THEN
ITER = -1
GO TO 160
END IF
*
ANRM = DLANGE('Frobenius',N,N,A,LDA,WORK)
EPS = DLAMCH('P')
BWDCTE = ANRM*EPS*MIN(4.0E00,SQRT(FLOAT(N))/6E00)
*
SPTSA = 1
SPTSX = SPTSA + N*N
*
DPTB = 1
DPTR = DPTB + N*NRHS
*
CALL DLACPY('All',N,NRHS,X,LDX,WORK(DPTB),N)
*
CALL DSLACONVERTGED2S(N,N,A,LDA,SWORK(SPTSA),N,INFO)
*
IF (INFO.NE.0) THEN
ITER = -2
GO TO 160
END IF
*
* Compute the LU factorization of SA.
*
CALL SGETRF(N,N,SWORK(SPTSA),N,IPIV,INFO)
*
IF (INFO.NE.0) THEN
ITER = -3
GO TO 160
END IF
*
* Solve the system A*X = B, overwriting B with X.
*
CALL DSLACONVERTGED2S(N,NRHS,WORK(DPTB),N,SWORK(SPTSX),N,INFO)
*
IF (INFO.NE.0) THEN
ITER = -2
GO TO 160
END IF
*
CALL SGETRS('No transpose',N,NRHS,SWORK(SPTSA),N,IPIV,
+ SWORK(SPTSX),N,INFO)
*
IF (INFO.NE.0) THEN
ITER = -4
GO TO 160
END IF
*
CALL DSLACONVERTGES2D(N,NRHS,SWORK(SPTSX),N,X,LDX,INFO)
*
CALL DLACPY('All',N,NRHS,WORK(DPTB),N,WORK(DPTR),N)
*
CALL DGEMM('No Transpose','No Transpose',N,NRHS,N,MONE,A,LDA,X,
+ LDX,ONE,WORK(DPTR),N)
*
OK = 1
DO I = 1,NRHS
XNRM = DNRM2(N,X(1,I),1)
RNRM = DNRM2(N,WORK(DPTR+ (I-1)*N),1)
BWD = RNRM/XNRM/BWDCTE
IF (BWD.GT.1) OK = 0
END DO
*
ITER = 1
*
DO 90 IITER = 1,ITERMAX
*
IF (OK.EQ.1) GO TO 160
*
CALL DSLACONVERTGED2S(N,NRHS,WORK(DPTR),N,SWORK(SPTSX),N,INFO)
*
IF (INFO.NE.0) THEN
ITER = -2
GO TO 160
END IF
*
CALL SGETRS('No transpose',N,NRHS,SWORK(SPTSA),N,IPIV,
+ SWORK(SPTSX),N,INFO)
*
IF (INFO.NE.0) THEN
ITER = -4
GO TO 160
END IF
*
CALL DSLACONVERTGES2D(N,NRHS,SWORK(SPTSX),N,WORK(DPTR),N,INFO)
*
CALL DAXPY(N*NRHS,ONE,WORK(DPTR),1,X,1)
*
CALL DLACPY('All',N,NRHS,WORK(DPTB),N,WORK(DPTR),N)
*
CALL DGEMM('No Transpose','No Transpose',N,NRHS,N,MONE,A,LDA,
+ X,LDX,ONE,WORK(DPTR),N)
*
*
OK = 1
DO I = 1,NRHS
XNRM = DNRM2(N,X(1,I),1)
RNRM = DNRM2(N,WORK(DPTR+ (I-1)*N),1)
BWD = RNRM/XNRM/BWDCTE
IF (BWD.GT.1) OK = 0
END DO
*
ITER = ITER + 1
*
90 CONTINUE
160 CONTINUE
*
IF (BWD.LT.1) RETURN
*
BWD = 0.0E0
IF (ITER.GE.0) THEN
ITER = -ITERMAX - 1
END IF
*
CALL DGETRF(N,N,A,LDA,IPIV,INFO)
*
CALL DLACPY('All',N,NRHS,WORK(DPTB),N,X,LDX)
*
IF (INFO.EQ.0) THEN
*
* Solve the system A*X = B, overwriting B with X.
*
CALL DGETRS('No transpose',N,NRHS,A,LDA,IPIV,X,LDX,INFO)
END IF
RETURN
*
* End of DGESV
*
END