dsgesv.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 363 行 · 第 1/2 页
HTML
363 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dsgesv.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.string { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.string a { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DSGESV.1"></a><a href="dsgesv.f.html#DSGESV.1">DSGESV</a>( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
+ SWORK, ITER, INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK PROTOTYPE driver routine (version 3.1.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> February 2007
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. WARNING: PROTOTYPE ..
</span><span class="comment">*</span><span class="comment"> This is an LAPACK PROTOTYPE routine which means that the
</span><span class="comment">*</span><span class="comment"> interface of this routine is likely to be changed in the future
</span><span class="comment">*</span><span class="comment"> based on community feedback.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> INTEGER INFO,ITER,LDA,LDB,LDX,N,NRHS
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> INTEGER IPIV(*)
REAL SWORK(*)
DOUBLE PRECISION A(LDA,*),B(LDB,*),WORK(N,*),X(LDX,*)
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DSGESV.27"></a><a href="dsgesv.f.html#DSGESV.1">DSGESV</a> computes the solution to a real system of linear equations
</span><span class="comment">*</span><span class="comment"> A * X = B,
</span><span class="comment">*</span><span class="comment"> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DSGESV.31"></a><a href="dsgesv.f.html#DSGESV.1">DSGESV</a> first attempts to factorize the matrix in SINGLE PRECISION
</span><span class="comment">*</span><span class="comment"> and use this factorization within an iterative refinement procedure to
</span><span class="comment">*</span><span class="comment"> produce a solution with DOUBLE PRECISION normwise backward error
</span><span class="comment">*</span><span class="comment"> quality (see below). If the approach fails the method switches to a
</span><span class="comment">*</span><span class="comment"> DOUBLE PRECISION factorization and solve.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The iterative refinement is not going to be a winning strategy if
</span><span class="comment">*</span><span class="comment"> the ratio SINGLE PRECISION performance over DOUBLE PRECISION performance
</span><span class="comment">*</span><span class="comment"> is too small. A reasonable strategy should take the number of right-hand
</span><span class="comment">*</span><span class="comment"> sides and the size of the matrix into account. This might be done with a
</span><span class="comment">*</span><span class="comment"> call to <a name="ILAENV.41"></a><a href="hfy-index.html#ILAENV">ILAENV</a> in the future. Up to now, we always try iterative refinement.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The iterative refinement process is stopped if
</span><span class="comment">*</span><span class="comment"> ITER > ITERMAX
</span><span class="comment">*</span><span class="comment"> or for all the RHS we have:
</span><span class="comment">*</span><span class="comment"> RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
</span><span class="comment">*</span><span class="comment"> where
</span><span class="comment">*</span><span class="comment"> o ITER is the number of the current iteration in the iterative
</span><span class="comment">*</span><span class="comment"> refinement process
</span><span class="comment">*</span><span class="comment"> o RNRM is the infinity-norm of the residual
</span><span class="comment">*</span><span class="comment"> o XNRM is the infinity-norm of the solution
</span><span class="comment">*</span><span class="comment"> o ANRM is the infinity-operator-norm of the matrix A
</span><span class="comment">*</span><span class="comment"> o EPS is the machine epsilon returned by <a name="DLAMCH.53"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>('Epsilon')
</span><span class="comment">*</span><span class="comment"> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of linear equations, i.e., the order of the
</span><span class="comment">*</span><span class="comment"> matrix A. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NRHS (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of right hand sides, i.e., the number of columns
</span><span class="comment">*</span><span class="comment"> of the matrix B. NRHS >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> A (input or input/ouptut) DOUBLE PRECISION array,
</span><span class="comment">*</span><span class="comment"> dimension (LDA,N)
</span><span class="comment">*</span><span class="comment"> On entry, the N-by-N coefficient matrix A.
</span><span class="comment">*</span><span class="comment"> On exit, if iterative refinement has been successfully used
</span><span class="comment">*</span><span class="comment"> (INFO.EQ.0 and ITER.GE.0, see description below), then A is
</span><span class="comment">*</span><span class="comment"> unchanged, if double precision factorization has been used
</span><span class="comment">*</span><span class="comment"> (INFO.EQ.0 and ITER.LT.0, see description below), then the
</span><span class="comment">*</span><span class="comment"> array A contains the factors L and U from the factorization
</span><span class="comment">*</span><span class="comment"> A = P*L*U; the unit diagonal elements of L are not stored.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDA (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array A. LDA >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IPIV (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment"> The pivot indices that define the permutation matrix P;
</span><span class="comment">*</span><span class="comment"> row i of the matrix was interchanged with row IPIV(i).
</span><span class="comment">*</span><span class="comment"> Corresponds either to the single precision factorization
</span><span class="comment">*</span><span class="comment"> (if INFO.EQ.0 and ITER.GE.0) or the double precision
</span><span class="comment">*</span><span class="comment"> factorization (if INFO.EQ.0 and ITER.LT.0).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
</span><span class="comment">*</span><span class="comment"> The N-by-NRHS matrix of right hand side matrix B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDB (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array B. LDB >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
</span><span class="comment">*</span><span class="comment"> If INFO = 0, the N-by-NRHS solution matrix X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDX (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array X. LDX >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) DOUBLE PRECISION array, dimension (N*NRHS)
</span><span class="comment">*</span><span class="comment"> This array is used to hold the residual vectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SWORK (workspace) REAL array, dimension (N*(N+NRHS))
</span><span class="comment">*</span><span class="comment"> This array is used to use the single precision matrix and the
</span><span class="comment">*</span><span class="comment"> right-hand sides or solutions in single precision.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ITER (output) INTEGER
</span><span class="comment">*</span><span class="comment"> < 0: iterative refinement has failed, double precision
</span><span class="comment">*</span><span class="comment"> factorization has been performed
</span><span class="comment">*</span><span class="comment"> -1 : taking into account machine parameters, N, NRHS, it
</span><span class="comment">*</span><span class="comment"> is a priori not worth working in SINGLE PRECISION
</span><span class="comment">*</span><span class="comment"> -2 : overflow of an entry when moving from double to
</span><span class="comment">*</span><span class="comment"> SINGLE PRECISION
</span><span class="comment">*</span><span class="comment"> -3 : failure of <a name="SGETRF.113"></a><a href="sgetrf.f.html#SGETRF.1">SGETRF</a>
</span><span class="comment">*</span><span class="comment"> -31: stop the iterative refinement after the 30th
</span><span class="comment">*</span><span class="comment"> iterations
</span><span class="comment">*</span><span class="comment"> > 0: iterative refinement has been sucessfully used.
</span><span class="comment">*</span><span class="comment"> Returns the number of iterations
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: successful exit
</span><span class="comment">*</span><span class="comment"> < 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment"> > 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is
</span><span class="comment">*</span><span class="comment"> exactly zero. The factorization has been completed,
</span><span class="comment">*</span><span class="comment"> but the factor U is exactly singular, so the solution
</span><span class="comment">*</span><span class="comment"> could not be computed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> DOUBLE PRECISION NEGONE,ONE
PARAMETER (NEGONE=-1.0D+0,ONE=1.0D+0)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL DOITREF
INTEGER I,IITER,ITERMAX,OK,PTSA,PTSX
DOUBLE PRECISION ANRM,BWDMAX,CTE,EPS,RNRM,XNRM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL DAXPY,DGEMM,<a name="DLACPY.139"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>,<a name="DLAG2S.139"></a><a href="dlag2s.f.html#DLAG2S.1">DLAG2S</a>,<a name="SLAG2D.139"></a><a href="slag2d.f.html#SLAG2D.1">SLAG2D</a>,
+ <a name="SGETRF.140"></a><a href="sgetrf.f.html#SGETRF.1">SGETRF</a>,<a name="SGETRS.140"></a><a href="sgetrs.f.html#SGETRS.1">SGETRS</a>,<a name="XERBLA.140"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> INTEGER IDAMAX
DOUBLE PRECISION <a name="DLAMCH.144"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>,<a name="DLANGE.144"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
EXTERNAL IDAMAX,<a name="DLAMCH.145"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>,<a name="DLANGE.145"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS,DBLE,MAX,SQRT
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span> ITERMAX = 30
BWDMAX = 1.0E+00
DOITREF = .TRUE.
<span class="comment">*</span><span class="comment">
</span> OK = 0
INFO = 0
ITER = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span> IF (N.LT.0) THEN
INFO = -1
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?