Annotation of rpl/lapack/lapack/dsgesv.f, revision 1.7
1.1 bertrand 1: SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
2: + SWORK, ITER, INFO )
3: *
1.5 bertrand 4: * -- LAPACK PROTOTYPE driver routine (version 3.2.2) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * February 2007
8: *
9: * ..
10: * .. Scalar Arguments ..
11: INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
12: * ..
13: * .. Array Arguments ..
14: INTEGER IPIV( * )
15: REAL SWORK( * )
16: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
17: + X( LDX, * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * DSGESV computes the solution to a real system of linear equations
24: * A * X = B,
25: * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
26: *
27: * DSGESV first attempts to factorize the matrix in SINGLE PRECISION
28: * and use this factorization within an iterative refinement procedure
29: * to produce a solution with DOUBLE PRECISION normwise backward error
30: * quality (see below). If the approach fails the method switches to a
31: * DOUBLE PRECISION factorization and solve.
32: *
33: * The iterative refinement is not going to be a winning strategy if
34: * the ratio SINGLE PRECISION performance over DOUBLE PRECISION
35: * performance is too small. A reasonable strategy should take the
36: * number of right-hand sides and the size of the matrix into account.
37: * This might be done with a call to ILAENV in the future. Up to now, we
38: * always try iterative refinement.
39: *
40: * The iterative refinement process is stopped if
41: * ITER > ITERMAX
42: * or for all the RHS we have:
43: * RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
44: * where
45: * o ITER is the number of the current iteration in the iterative
46: * refinement process
47: * o RNRM is the infinity-norm of the residual
48: * o XNRM is the infinity-norm of the solution
49: * o ANRM is the infinity-operator-norm of the matrix A
50: * o EPS is the machine epsilon returned by DLAMCH('Epsilon')
51: * The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
52: * respectively.
53: *
54: * Arguments
55: * =========
56: *
57: * N (input) INTEGER
58: * The number of linear equations, i.e., the order of the
59: * matrix A. N >= 0.
60: *
61: * NRHS (input) INTEGER
62: * The number of right hand sides, i.e., the number of columns
63: * of the matrix B. NRHS >= 0.
64: *
1.5 bertrand 65: * A (input/output) DOUBLE PRECISION array,
1.1 bertrand 66: * dimension (LDA,N)
67: * On entry, the N-by-N coefficient matrix A.
68: * On exit, if iterative refinement has been successfully used
69: * (INFO.EQ.0 and ITER.GE.0, see description below), then A is
70: * unchanged, if double precision factorization has been used
71: * (INFO.EQ.0 and ITER.LT.0, see description below), then the
72: * array A contains the factors L and U from the factorization
73: * A = P*L*U; the unit diagonal elements of L are not stored.
74: *
75: * LDA (input) INTEGER
76: * The leading dimension of the array A. LDA >= max(1,N).
77: *
78: * IPIV (output) INTEGER array, dimension (N)
79: * The pivot indices that define the permutation matrix P;
80: * row i of the matrix was interchanged with row IPIV(i).
81: * Corresponds either to the single precision factorization
82: * (if INFO.EQ.0 and ITER.GE.0) or the double precision
83: * factorization (if INFO.EQ.0 and ITER.LT.0).
84: *
85: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
86: * The N-by-NRHS right hand side matrix B.
87: *
88: * LDB (input) INTEGER
89: * The leading dimension of the array B. LDB >= max(1,N).
90: *
91: * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
92: * If INFO = 0, the N-by-NRHS solution matrix X.
93: *
94: * LDX (input) INTEGER
95: * The leading dimension of the array X. LDX >= max(1,N).
96: *
1.5 bertrand 97: * WORK (workspace) DOUBLE PRECISION array, dimension (N,NRHS)
1.1 bertrand 98: * This array is used to hold the residual vectors.
99: *
100: * SWORK (workspace) REAL array, dimension (N*(N+NRHS))
101: * This array is used to use the single precision matrix and the
102: * right-hand sides or solutions in single precision.
103: *
104: * ITER (output) INTEGER
105: * < 0: iterative refinement has failed, double precision
106: * factorization has been performed
107: * -1 : the routine fell back to full precision for
108: * implementation- or machine-specific reasons
109: * -2 : narrowing the precision induced an overflow,
110: * the routine fell back to full precision
111: * -3 : failure of SGETRF
112: * -31: stop the iterative refinement after the 30th
113: * iterations
114: * > 0: iterative refinement has been sucessfully used.
115: * Returns the number of iterations
116: *
117: * INFO (output) INTEGER
118: * = 0: successful exit
119: * < 0: if INFO = -i, the i-th argument had an illegal value
120: * > 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is
121: * exactly zero. The factorization has been completed,
122: * but the factor U is exactly singular, so the solution
123: * could not be computed.
124: *
125: * =========
126: *
127: * .. Parameters ..
128: LOGICAL DOITREF
129: PARAMETER ( DOITREF = .TRUE. )
130: *
131: INTEGER ITERMAX
132: PARAMETER ( ITERMAX = 30 )
133: *
134: DOUBLE PRECISION BWDMAX
135: PARAMETER ( BWDMAX = 1.0E+00 )
136: *
137: DOUBLE PRECISION NEGONE, ONE
138: PARAMETER ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
139: *
140: * .. Local Scalars ..
141: INTEGER I, IITER, PTSA, PTSX
142: DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
143: *
144: * .. External Subroutines ..
145: EXTERNAL DAXPY, DGEMM, DLACPY, DLAG2S, SLAG2D, SGETRF,
146: + SGETRS, XERBLA
147: * ..
148: * .. External Functions ..
149: INTEGER IDAMAX
150: DOUBLE PRECISION DLAMCH, DLANGE
151: EXTERNAL IDAMAX, DLAMCH, DLANGE
152: * ..
153: * .. Intrinsic Functions ..
154: INTRINSIC ABS, DBLE, MAX, SQRT
155: * ..
156: * .. Executable Statements ..
157: *
158: INFO = 0
159: ITER = 0
160: *
161: * Test the input parameters.
162: *
163: IF( N.LT.0 ) THEN
164: INFO = -1
165: ELSE IF( NRHS.LT.0 ) THEN
166: INFO = -2
167: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
168: INFO = -4
169: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
170: INFO = -7
171: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
172: INFO = -9
173: END IF
174: IF( INFO.NE.0 ) THEN
175: CALL XERBLA( 'DSGESV', -INFO )
176: RETURN
177: END IF
178: *
179: * Quick return if (N.EQ.0).
180: *
181: IF( N.EQ.0 )
182: + RETURN
183: *
184: * Skip single precision iterative refinement if a priori slower
185: * than double precision factorization.
186: *
187: IF( .NOT.DOITREF ) THEN
188: ITER = -1
189: GO TO 40
190: END IF
191: *
192: * Compute some constants.
193: *
194: ANRM = DLANGE( 'I', N, N, A, LDA, WORK )
195: EPS = DLAMCH( 'Epsilon' )
196: CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
197: *
198: * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
199: *
200: PTSA = 1
201: PTSX = PTSA + N*N
202: *
203: * Convert B from double precision to single precision and store the
204: * result in SX.
205: *
206: CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
207: *
208: IF( INFO.NE.0 ) THEN
209: ITER = -2
210: GO TO 40
211: END IF
212: *
213: * Convert A from double precision to single precision and store the
214: * result in SA.
215: *
216: CALL DLAG2S( N, N, A, LDA, SWORK( PTSA ), N, INFO )
217: *
218: IF( INFO.NE.0 ) THEN
219: ITER = -2
220: GO TO 40
221: END IF
222: *
223: * Compute the LU factorization of SA.
224: *
225: CALL SGETRF( N, N, SWORK( PTSA ), N, IPIV, INFO )
226: *
227: IF( INFO.NE.0 ) THEN
228: ITER = -3
229: GO TO 40
230: END IF
231: *
232: * Solve the system SA*SX = SB.
233: *
234: CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
235: + SWORK( PTSX ), N, INFO )
236: *
237: * Convert SX back to double precision
238: *
239: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
240: *
241: * Compute R = B - AX (R is WORK).
242: *
243: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
244: *
245: CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, A,
246: + LDA, X, LDX, ONE, WORK, N )
247: *
248: * Check whether the NRHS normwise backward errors satisfy the
249: * stopping criterion. If yes, set ITER=0 and return.
250: *
251: DO I = 1, NRHS
252: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
253: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
254: IF( RNRM.GT.XNRM*CTE )
255: + GO TO 10
256: END DO
257: *
258: * If we are here, the NRHS normwise backward errors satisfy the
259: * stopping criterion. We are good to exit.
260: *
261: ITER = 0
262: RETURN
263: *
264: 10 CONTINUE
265: *
266: DO 30 IITER = 1, ITERMAX
267: *
268: * Convert R (in WORK) from double precision to single precision
269: * and store the result in SX.
270: *
271: CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
272: *
273: IF( INFO.NE.0 ) THEN
274: ITER = -2
275: GO TO 40
276: END IF
277: *
278: * Solve the system SA*SX = SR.
279: *
280: CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
281: + SWORK( PTSX ), N, INFO )
282: *
283: * Convert SX back to double precision and update the current
284: * iterate.
285: *
286: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
287: *
288: DO I = 1, NRHS
289: CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
290: END DO
291: *
292: * Compute R = B - AX (R is WORK).
293: *
294: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
295: *
296: CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE,
297: + A, LDA, X, LDX, ONE, WORK, N )
298: *
299: * Check whether the NRHS normwise backward errors satisfy the
300: * stopping criterion. If yes, set ITER=IITER>0 and return.
301: *
302: DO I = 1, NRHS
303: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
304: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
305: IF( RNRM.GT.XNRM*CTE )
306: + GO TO 20
307: END DO
308: *
309: * If we are here, the NRHS normwise backward errors satisfy the
310: * stopping criterion, we are good to exit.
311: *
312: ITER = IITER
313: *
314: RETURN
315: *
316: 20 CONTINUE
317: *
318: 30 CONTINUE
319: *
320: * If we are at this place of the code, this is because we have
321: * performed ITER=ITERMAX iterations and never satisified the
322: * stopping criterion, set up the ITER flag accordingly and follow up
323: * on double precision routine.
324: *
325: ITER = -ITERMAX - 1
326: *
327: 40 CONTINUE
328: *
329: * Single-precision iterative refinement failed to converge to a
330: * satisfactory solution, so we resort to double precision.
331: *
332: CALL DGETRF( N, N, A, LDA, IPIV, INFO )
333: *
334: IF( INFO.NE.0 )
335: + RETURN
336: *
337: CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
338: CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
339: + INFO )
340: *
341: RETURN
342: *
343: * End of DSGESV.
344: *
345: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsgesv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack