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