 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dsgesv()

 subroutine dsgesv ( integer N, integer NRHS, double precision, dimension( lda, * ) A, integer LDA, integer, dimension( * ) IPIV, double precision, dimension( ldb, * ) B, integer LDB, double precision, dimension( ldx, * ) X, integer LDX, double precision, dimension( n, * ) WORK, real, dimension( * ) SWORK, integer ITER, integer INFO )

DSGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement)

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 attempts 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 switches to a
DOUBLE PRECISION factorization and solve.

The iterative refinement is not going to be a winning strategy if
the ratio SINGLE PRECISION performance over DOUBLE PRECISION
performance is too small. A reasonable strategy should take the
number of right-hand sides and the size of the matrix into account.
This might be done with a call to ILAENV in the future. Up to now, we
always try iterative refinement.

The iterative refinement process is stopped if
ITER > ITERMAX
or for all the RHS we have:
RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
where
o ITER is the number of the current iteration in the iterative
refinement process
o RNRM is the infinity-norm of the residual
o XNRM is the infinity-norm of the solution
o ANRM is the infinity-operator-norm of the matrix A
o EPS is the machine epsilon returned by DLAMCH('Epsilon')
The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
respectively.```
Parameters
 [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 matrix B. NRHS >= 0.``` [in,out] A ``` A is DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N coefficient matrix A. On exit, if iterative refinement has been successfully used (INFO = 0 and ITER >= 0, see description below), then A is unchanged, if double precision factorization has been used (INFO = 0 and ITER < 0, see description below), then the array A contains the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] IPIV ``` IPIV is INTEGER array, dimension (N) The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i). Corresponds either to the single precision factorization (if INFO = 0 and ITER >= 0) or the double precision factorization (if INFO = 0 and ITER < 0).``` [in] B ``` B is DOUBLE PRECISION array, dimension (LDB,NRHS) The N-by-NRHS right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [out] X ``` X is DOUBLE PRECISION array, dimension (LDX,NRHS) If INFO = 0, the N-by-NRHS solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] WORK ``` WORK is DOUBLE PRECISION array, dimension (N,NRHS) This array is used to hold the residual vectors.``` [out] SWORK ``` SWORK is REAL array, dimension (N*(N+NRHS)) This array is used to use the single precision matrix and the right-hand sides or solutions in single precision.``` [out] ITER ``` ITER is INTEGER < 0: iterative refinement has failed, double precision factorization has been performed -1 : the routine fell back to full precision for implementation- or machine-specific reasons -2 : narrowing the precision induced an overflow, the routine fell back to full precision -3 : failure of SGETRF -31: stop the iterative refinement after the 30th iterations > 0: iterative refinement has been successfully used. Returns the number of iterations``` [out] INFO ``` INFO is 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.```

Definition at line 193 of file dsgesv.f.

195*
196* -- LAPACK driver routine --
197* -- LAPACK is a software package provided by Univ. of Tennessee, --
198* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199*
200* .. Scalar Arguments ..
201 INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
202* ..
203* .. Array Arguments ..
204 INTEGER IPIV( * )
205 REAL SWORK( * )
206 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
207 \$ X( LDX, * )
208* ..
209*
210* =====================================================================
211*
212* .. Parameters ..
213 LOGICAL DOITREF
214 parameter( doitref = .true. )
215*
216 INTEGER ITERMAX
217 parameter( itermax = 30 )
218*
219 DOUBLE PRECISION BWDMAX
220 parameter( bwdmax = 1.0e+00 )
221*
222 DOUBLE PRECISION NEGONE, ONE
223 parameter( negone = -1.0d+0, one = 1.0d+0 )
224*
225* .. Local Scalars ..
226 INTEGER I, IITER, PTSA, PTSX
227 DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
228*
229* .. External Subroutines ..
230 EXTERNAL daxpy, dgemm, dlacpy, dlag2s, dgetrf, dgetrs,
232* ..
233* .. External Functions ..
234 INTEGER IDAMAX
235 DOUBLE PRECISION DLAMCH, DLANGE
236 EXTERNAL idamax, dlamch, dlange
237* ..
238* .. Intrinsic Functions ..
239 INTRINSIC abs, dble, max, sqrt
240* ..
241* .. Executable Statements ..
242*
243 info = 0
244 iter = 0
245*
246* Test the input parameters.
247*
248 IF( n.LT.0 ) THEN
249 info = -1
250 ELSE IF( nrhs.LT.0 ) THEN
251 info = -2
252 ELSE IF( lda.LT.max( 1, n ) ) THEN
253 info = -4
254 ELSE IF( ldb.LT.max( 1, n ) ) THEN
255 info = -7
256 ELSE IF( ldx.LT.max( 1, n ) ) THEN
257 info = -9
258 END IF
259 IF( info.NE.0 ) THEN
260 CALL xerbla( 'DSGESV', -info )
261 RETURN
262 END IF
263*
264* Quick return if (N.EQ.0).
265*
266 IF( n.EQ.0 )
267 \$ RETURN
268*
269* Skip single precision iterative refinement if a priori slower
270* than double precision factorization.
271*
272 IF( .NOT.doitref ) THEN
273 iter = -1
274 GO TO 40
275 END IF
276*
277* Compute some constants.
278*
279 anrm = dlange( 'I', n, n, a, lda, work )
280 eps = dlamch( 'Epsilon' )
281 cte = anrm*eps*sqrt( dble( n ) )*bwdmax
282*
283* Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
284*
285 ptsa = 1
286 ptsx = ptsa + n*n
287*
288* Convert B from double precision to single precision and store the
289* result in SX.
290*
291 CALL dlag2s( n, nrhs, b, ldb, swork( ptsx ), n, info )
292*
293 IF( info.NE.0 ) THEN
294 iter = -2
295 GO TO 40
296 END IF
297*
298* Convert A from double precision to single precision and store the
299* result in SA.
300*
301 CALL dlag2s( n, n, a, lda, swork( ptsa ), n, info )
302*
303 IF( info.NE.0 ) THEN
304 iter = -2
305 GO TO 40
306 END IF
307*
308* Compute the LU factorization of SA.
309*
310 CALL sgetrf( n, n, swork( ptsa ), n, ipiv, info )
311*
312 IF( info.NE.0 ) THEN
313 iter = -3
314 GO TO 40
315 END IF
316*
317* Solve the system SA*SX = SB.
318*
319 CALL sgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
320 \$ swork( ptsx ), n, info )
321*
322* Convert SX back to double precision
323*
324 CALL slag2d( n, nrhs, swork( ptsx ), n, x, ldx, info )
325*
326* Compute R = B - AX (R is WORK).
327*
328 CALL dlacpy( 'All', n, nrhs, b, ldb, work, n )
329*
330 CALL dgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone, a,
331 \$ lda, x, ldx, one, work, n )
332*
333* Check whether the NRHS normwise backward errors satisfy the
334* stopping criterion. If yes, set ITER=0 and return.
335*
336 DO i = 1, nrhs
337 xnrm = abs( x( idamax( n, x( 1, i ), 1 ), i ) )
338 rnrm = abs( work( idamax( n, work( 1, i ), 1 ), i ) )
339 IF( rnrm.GT.xnrm*cte )
340 \$ GO TO 10
341 END DO
342*
343* If we are here, the NRHS normwise backward errors satisfy the
344* stopping criterion. We are good to exit.
345*
346 iter = 0
347 RETURN
348*
349 10 CONTINUE
350*
351 DO 30 iiter = 1, itermax
352*
353* Convert R (in WORK) from double precision to single precision
354* and store the result in SX.
355*
356 CALL dlag2s( n, nrhs, work, n, swork( ptsx ), n, info )
357*
358 IF( info.NE.0 ) THEN
359 iter = -2
360 GO TO 40
361 END IF
362*
363* Solve the system SA*SX = SR.
364*
365 CALL sgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
366 \$ swork( ptsx ), n, info )
367*
368* Convert SX back to double precision and update the current
369* iterate.
370*
371 CALL slag2d( n, nrhs, swork( ptsx ), n, work, n, info )
372*
373 DO i = 1, nrhs
374 CALL daxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
375 END DO
376*
377* Compute R = B - AX (R is WORK).
378*
379 CALL dlacpy( 'All', n, nrhs, b, ldb, work, n )
380*
381 CALL dgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone,
382 \$ a, lda, x, ldx, one, work, n )
383*
384* Check whether the NRHS normwise backward errors satisfy the
385* stopping criterion. If yes, set ITER=IITER>0 and return.
386*
387 DO i = 1, nrhs
388 xnrm = abs( x( idamax( n, x( 1, i ), 1 ), i ) )
389 rnrm = abs( work( idamax( n, work( 1, i ), 1 ), i ) )
390 IF( rnrm.GT.xnrm*cte )
391 \$ GO TO 20
392 END DO
393*
394* If we are here, the NRHS normwise backward errors satisfy the
395* stopping criterion, we are good to exit.
396*
397 iter = iiter
398*
399 RETURN
400*
401 20 CONTINUE
402*
403 30 CONTINUE
404*
405* If we are at this place of the code, this is because we have
406* performed ITER=ITERMAX iterations and never satisfied the
407* stopping criterion, set up the ITER flag accordingly and follow up
408* on double precision routine.
409*
410 iter = -itermax - 1
411*
412 40 CONTINUE
413*
414* Single-precision iterative refinement failed to converge to a
415* satisfactory solution, so we resort to double precision.
416*
417 CALL dgetrf( n, n, a, lda, ipiv, info )
418*
419 IF( info.NE.0 )
420 \$ RETURN
421*
422 CALL dlacpy( 'All', n, nrhs, b, ldb, x, ldx )
423 CALL dgetrs( 'No transpose', n, nrhs, a, lda, ipiv, x, ldx,
424 \$ info )
425*
426 RETURN
427*
428* End of DSGESV
429*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine slag2d(M, N, SA, LDSA, A, LDA, INFO)
SLAG2D converts a single precision matrix to a double precision matrix.
Definition: slag2d.f:104
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
subroutine dgetrf(M, N, A, LDA, IPIV, INFO)
DGETRF
Definition: dgetrf.f:108
subroutine dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
DGETRS
Definition: dgetrs.f:121
subroutine dlag2s(M, N, A, LDA, SA, LDSA, INFO)
DLAG2S converts a double precision matrix to a single precision matrix.
Definition: dlag2s.f:108
subroutine sgetrf(M, N, A, LDA, IPIV, INFO)
SGETRF
Definition: sgetrf.f:108
subroutine sgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SGETRS
Definition: sgetrs.f:121
Here is the call graph for this function:
Here is the caller graph for this function: