1DSGESV ‐ the solution to a real system of linear equations A * X
2= B, SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX,
3WORK,
4 + SWORK, ITER, INFO)
5 INTEGER INFO,ITER,LDA,LDB,LDX,N,NRHS
6 INTEGER IPIV(*)
7 REAL SWORK(*)
8 DOUBLE PRECISION A(LDA,*),B(LDB,*),WORK(N,*),X(LDX,*) DSGESV
9computes the solution to a real system of linear equations
10 A * X = B, where A is an N‐by‐N matrix and X and B are N‐by‐
11NRHS matrices.
12DSGESV first attempts to factorize the matrix in SINGLE PRECISION
14and use this factorization within an iterative refinement proce‐
15dure to produce a solution with DOUBLE PRECISION normwise back‐
16ward error quality (see below). If the approach fails the method
17switches to a DOUBLE PRECISION factorization and solve.
18The iterative refinement is not going to be a winning strategy if
20the ratio SINGLE PRECISION performance over DOUBLE PRECISION per‐
21formance is too small. A reasonable strategy should take the num‐
22ber of right‐hand sides and the size of the matrix into account.
23This might be done with a call to ILAENV in the future. Up to
24now, we always try iterative refinement.
25The iterative refinement process is stopped if
27 ITER > ITERMAX
28or for all the RHS we have:
29 RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
30where
31 o ITER is the number of the current iteration in the itera‐
32tive
33 refinement process
34 o RNRM is the infinity‐norm of the residual
35 o XNRM is the infinity‐norm of the solution
36 o ANRM is the infinity‐operator‐norm of the matrix A
37 o EPS is the machine epsilon returned by DLAMCH('Epsilon')
38The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respec‐
39tively.
40N (input) INTEGER The number of linear equations, i.e., the
42order of the matrix A. N >= 0. NRHS (input) INTEGER The num‐
43ber of right hand sides, i.e., the number of columns of the ma‐
44trix B. NRHS >= 0. A (input or input/ouptut) DOUBLE PRE‐
45CISION array, dimension (LDA,N) On entry, the N‐by‐N coefficient
46matrix A. On exit, if iterative refinement has been successfully
47used (INFO.EQ.0 and ITER.GE.0, see description below), then A is
48unchanged, if double precision factorization has been used (IN‐
49FO.EQ.0 and ITER.LT.0, see description below), then the array A
50contains the factors L and U from the factorization A = P*L*U;
51the unit diagonal elements of L are not stored. LDA (input)
52INTEGER The leading dimension of the array A. LDA >= max(1,N).
53IPIV (output) INTEGER array, dimension (N) The pivot indices
54that define the permutation matrix P; row i of the matrix was in‐
55terchanged with row IPIV(i). Corresponds either to the single
56precision factorization (if INFO.EQ.0 and ITER.GE.0) or the dou‐
57ble precision factorization (if INFO.EQ.0 and ITER.LT.0). B
58(input) DOUBLE PRECISION array, dimension (LDB,NRHS) The N‐by‐
59NRHS matrix of right hand side matrix B. LDB (input) INTEGER
60The leading dimension of the array B. LDB >= max(1,N). X
61(output) DOUBLE PRECISION array, dimension (LDX,NRHS) If INFO =
620, the N‐by‐NRHS solution matrix X. LDX (input) INTEGER The
63leading dimension of the array X. LDX >= max(1,N). WORK
64(workspace) DOUBLE PRECISION array, dimension (N*NRHS) This array
65is used to hold the residual vectors. SWORK (workspace) REAL
66array, dimension (N*(N+NRHS)) This array is used to use the sin‐
67gle precision matrix and the right‐hand sides or solutions in
68single precision. ITER (output) INTEGER < 0: iterative re‐
69finement has failed, double precision factorization has been per‐
70formed ‐1 : taking into account machine parameters, N, NRHS, it
71is a priori not worth working in SINGLE PRECISION ‐2 : overflow
72of an entry when moving from double to SINGLE PRECISION ‐3 :
73failure of SGETRF
74‐31: stop the iterative refinement after the 30th iterations > 0:
75iterative refinement has been sucessfully used. Returns the num‐
76ber of iterations INFO (output) INTEGER = 0: successful exit
77< 0: if INFO = ‐i, the i‐th argument had an illegal value
78> 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is exactly
79zero. The factorization has been completed, but the factor U is
80exactly singular, so the solution could not be computed.
81=========
83