1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
|
<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN">
<!--Converted with LaTeX2HTML 96.1-h (September 30, 1996) by Nikos Drakos (nikos@cbl.leeds.ac.uk), CBLU, University of Leeds -->
<HTML>
<HEAD>
<TITLE>Error Bounds for Linear Equation Solving</TITLE>
<META NAME="description" CONTENT="Error Bounds for Linear Equation Solving">
<META NAME="keywords" CONTENT="slug">
<META NAME="resource-type" CONTENT="document">
<META NAME="distribution" CONTENT="global">
<LINK REL=STYLESHEET HREF="slug.css">
</HEAD>
<BODY LANG="EN" >
<A NAME="tex2html3950" HREF="node140.html"><IMG WIDTH=37 HEIGHT=24 ALIGN=BOTTOM ALT="next" SRC="http://www.netlib.org/utk/icons/next_motif.gif"></A> <A NAME="tex2html3948" HREF="node132.html"><IMG WIDTH=26 HEIGHT=24 ALIGN=BOTTOM ALT="up" SRC="http://www.netlib.org/utk/icons/up_motif.gif"></A> <A NAME="tex2html3942" HREF="node138.html"><IMG WIDTH=63 HEIGHT=24 ALIGN=BOTTOM ALT="previous" SRC="http://www.netlib.org/utk/icons/previous_motif.gif"></A> <A NAME="tex2html3952" HREF="node1.html"><IMG WIDTH=65 HEIGHT=24 ALIGN=BOTTOM ALT="contents" SRC="http://www.netlib.org/utk/icons/contents_motif.gif"></A> <A NAME="tex2html3953" HREF="node190.html"><IMG WIDTH=43 HEIGHT=24 ALIGN=BOTTOM ALT="index" SRC="http://www.netlib.org/utk/icons/index_motif.gif"></A> <BR>
<B> Next:</B> <A NAME="tex2html3951" HREF="node140.html">Error Bounds for Linear </A>
<B>Up:</B> <A NAME="tex2html3949" HREF="node132.html">Accuracy and Stability</A>
<B> Previous:</B> <A NAME="tex2html3943" HREF="node138.html">Improved Error Bounds</A>
<BR> <P>
<H1><A NAME="SECTION04650000000000000000">Error Bounds for Linear Equation Solving</A></H1>
<A NAME="secAxb"> </A>
<P>
Let <I>Ax</I>=<I>b</I> be the system to be solved, and <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> the computed
solution. Let <I>n</I> be the dimension of <I>A</I>.
An approximate error bound<A NAME="5088"> </A>
for <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> may be obtained in one of the following two ways,
depending on whether the solution is computed by a simple driver or
an expert driver:
<P>
<OL>
<LI> Suppose that <I>Ax</I>=<I>b</I> is solved using the simple driver PSGESV
<A NAME="5091"> </A><A NAME="5092"> </A>
(section <A HREF="node44.html#subsecdrivelineq">3.2.1</A>).
Then the approximate error bound<A NAME="tex2html1256" HREF="footnode.html#8057"><IMG ALIGN=BOTTOM ALT="gif" SRC="http://www.netlib.org/utk/icons/foot_motif.gif"></A>
<BR><BR>
<BR><IMG WIDTH=322 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18202" SRC="img570.gif"><BR>
can be computed by the following code fragment.
<P>
<PRE> EPSMCH = PSLAMCH( ICTXT, 'E' )
* Get infinity-norm of A
ANORM = PSLANGE( 'I', N, N, A, IA, JA, DESCA, WORK )
* Solve system; The solution X overwrites B
CALL PSGESV( N, 1, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO )
IF( INFO.GT.0 ) THEN
PRINT *,'Singular Matrix'
ELSE IF( N.GT.0 ) THEN
* Get reciprocal condition number RCOND of A
CALL PSGECON( 'I', N, A, IA, JA, DESCA, ANORM, RCOND, WORK,
$ LWORK, IWORK, LIWORK, INFO )
RCOND = MAX( RCOND, EPSMCH )
ERRBD = EPSMCH / RCOND
END IF</PRE>
<A NAME="5101"> </A>
<P>
For example, suppose
<BR><BR>
<IMG WIDTH=315 HEIGHT=30 ALIGN=MIDDLE ALT="tex2html_wrap_inline18242" SRC="img571.gif">,
<BR><IMG WIDTH=442 HEIGHT=67 ALIGN=BOTTOM ALT="displaymath18203" SRC="img572.gif"><BR>
Then (to 4 decimal places)
<BR><IMG WIDTH=399 HEIGHT=67 ALIGN=BOTTOM ALT="displaymath18204" SRC="img573.gif"><BR>
<IMG WIDTH=144 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18244" SRC="img574.gif">,
<IMG WIDTH=155 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18246" SRC="img575.gif">,
the true reciprocal condition number <IMG WIDTH=103 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18248" SRC="img576.gif">,
<IMG WIDTH=137 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18250" SRC="img577.gif">, and the true error
<IMG WIDTH=86 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18252" SRC="img578.gif">.
<A NAME="5120"> </A>
<LI> Suppose that <I>Ax</I>=<I>b</I> is solved using the expert driver PSGESVX
(section <A HREF="node44.html#subsecdrivelineq">3.2.1</A>).
<A NAME="5122"> </A><A NAME="5123"> </A>
This routine provides an explicit error bound <TT>FERR</TT>, measured
with the infinity-norm:
<A NAME="5125"> </A>
<BR><IMG WIDTH=317 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18205" SRC="img579.gif"><BR>
For example, the following code fragment solves
<I>Ax</I>=<I>b</I> and computes an approximate error bound <TT>FERR</TT>:
<P>
<PRE> CALL PSGESVX( 'E', 'N', N, 1, A, IA, JA, DESCA, AF, IAF, JAF,
$ DESCAF, IPIV, EQUED, R, C, B, IB, JB, DESCB, X, IX,
$ JX, DESCX, RCOND, FERR, BERR, WORK, LWORK, IWORK,
$ LIWORK, INFO )
IF( INFO.GT.0 ) PRINT *,'(Nearly) Singular Matrix'</PRE>
<P>
For the same <TT>A</TT> and <TT>b</TT> as above,
<IMG WIDTH=125 HEIGHT=76 ALIGN=MIDDLE ALT="tex2html_wrap_inline18258" SRC="img580.gif">,
<IMG WIDTH=128 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18260" SRC="img581.gif">,
and the actual error is <IMG WIDTH=68 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18262" SRC="img582.gif">.
<P>
</OL>
<P>
This example illustrates
that the expert driver provides an error bound with less programming
effort than the simple driver, and also that it may produce a significantly
more accurate answer.
<P>
Similar code fragments, with obvious adaptations,
may be used with all the driver routines PxPOSV and PxPOSVX
<A NAME="5139"> </A><A NAME="5140"> </A><A NAME="5141"> </A><A NAME="5142"> </A>
in Table <A HREF="node44.html#tabdrivelineq">3.2</A>.
For example, if a symmetric positive definite or Hermitian positive
definite system is solved by using the simple driver PxPOSV,
<A NAME="5144"> </A><A NAME="5145"> </A><A NAME="5146"> </A><A NAME="5147"> </A>
then PxLANSY
or PxLANHE, respectively, must
be used to compute <TT>ANORM</TT>, and
PxPOCON<A NAME="5149"> </A><A NAME="5150"> </A><A NAME="5151"> </A><A NAME="5152"> </A> must
be used to compute <TT>RCOND</TT>.
<P>
The drivers
PxGBSV
<A NAME="5154"> </A><A NAME="5155"> </A><A NAME="5156"> </A><A NAME="5157"> </A>
(for solving general band matrices with partial pivoting),
PxPBSV <A NAME="5158"> </A><A NAME="5159"> </A><A NAME="5160"> </A><A NAME="5161"> </A>
(for solving positive definite band matrices) and
PxPTSV <A NAME="5162"> </A><A NAME="5163"> </A><A NAME="5164"> </A><A NAME="5165"> </A>
(for solving positive definite tridiagonal matrices),
do not yet have the corresponding routines needed to compute
error bounds, namely,
PxLAnHE to compute <TT>ANORM</TT> and PxyyCON to compute <TT>RCOND</TT>.
<P>
The drivers
PxDBSV
<A NAME="5168"> </A><A NAME="5169"> </A><A NAME="5170"> </A><A NAME="5171"> </A>
(for solving general band matrices) and
PxDTSV <A NAME="5172"> </A><A NAME="5173"> </A><A NAME="5174"> </A><A NAME="5175"> </A>
(for solving general tridiagonal matrices) do not pivot for numerical
stability, and so may be faster but less accurate than their
pivoting counterparts above. These routines may be used safely when
any diagonal pivot sequence leads to a stable factorization;
diagonally dominant matrices and symmetric positive definite matrices
[<A HREF="node189.html#GVL2">71</A>] have this property, for example.
<P>
<B>Further Details: Error Bounds for Linear Equation Solving</B><A NAME="secbackgroundAxb"> </A>
<P>
The conventional error analysis of linear
equation<A NAME="5179"> </A> solving goes as follows.
Let <I>Ax</I>=<I>b</I> be the system to be solved. Let <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> be the solution
computed by ScaLAPACK (or LAPACK) using any of their linear equation solvers.
Let <I>r</I> be
the residual <IMG WIDTH=82 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline18270" SRC="img583.gif">. In the absence of rounding error, <I>r</I>
would be zero and <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> would equal <I>x</I>; with rounding error, one can
only say the following:
<P>
<BLOCKQUOTE> The normwise backward error of the computed solution <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif">,
<A NAME="5185"> </A>
<A NAME="5186"> </A>
with respect to the infinity norm,
is the pair <I>E</I>,<I>f</I>, which minimizes
<BR><IMG WIDTH=329 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18206" SRC="img584.gif"><BR>
subject to the constraint <IMG WIDTH=132 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18075" SRC="img546.gif">.
The minimal value of
<IMG WIDTH=137 HEIGHT=38 ALIGN=MIDDLE ALT="tex2html_wrap_inline18284" SRC="img585.gif">
is given by
<BR><IMG WIDTH=360 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18207" SRC="img586.gif"><BR>
One can show that the computed solution <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif">
satisfies <IMG WIDTH=100 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18288" SRC="img587.gif">,
where <I>p</I>(<I>n</I>) is a modestly growing function of <I>n</I>.
The corresponding condition number is
<IMG WIDTH=192 HEIGHT=30 ALIGN=MIDDLE ALT="tex2html_wrap_inline18294" SRC="img588.gif">.
<A NAME="5200"> </A>
The error <IMG WIDTH=40 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline18296" SRC="img589.gif"> is bounded by
<BR><IMG WIDTH=393 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18208" SRC="img590.gif"><BR>
In the first code fragment in the preceding section, <IMG WIDTH=44 HEIGHT=25 ALIGN=MIDDLE ALT="tex2html_wrap_inline18300" SRC="img591.gif">,
which is <IMG WIDTH=86 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18302" SRC="img592.gif"> in the numerical example,
is approximated by <IMG WIDTH=172 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18304" SRC="img593.gif">.
Approximations<A NAME="5209"> </A>
of <IMG WIDTH=48 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18121" SRC="img555.gif"> -- or, strictly speaking, its reciprocal <TT>RCOND</TT> --
are returned by computational routines
PxyyCON (section <A HREF="node51.html#subseccomplineq">3.3.1</A>) or driver routines
<A NAME="5212"> </A><A NAME="5213"> </A><A NAME="5214"> </A><A NAME="5215"> </A>
<A NAME="5216"> </A><A NAME="5217"> </A><A NAME="5218"> </A><A NAME="5219"> </A>
<A NAME="5220"> </A><A NAME="5221"> </A><A NAME="5222"> </A><A NAME="5223"> </A>
PxyySVX (section <A HREF="node44.html#subsecdrivelineq">3.2.1</A>). The code fragment
makes sure <TT>RCOND</TT> is at least <IMG WIDTH=24 HEIGHT=8 ALIGN=BOTTOM ALT="tex2html_wrap_inline18308" SRC="img594.gif"> <TT>EPSMCH</TT> to
avoid overflow in computing
<TT>ERRBD</TT>.<A NAME="5228"> </A><A NAME="5229"> </A>
This limits
<TT>ERRBD</TT> to a maximum of 1, which is no loss of generality because
a relative error of 1 or more indicates the same thing:
<A NAME="5231"> </A><A NAME="5232"> </A>
a complete loss of accuracy. <A NAME="5233"> </A>
Note that the
value of <TT>RCOND</TT> returned by PxyySVX may apply to a linear
system obtained from <I>Ax</I>=<I>b</I> by <EM>equilibration</EM>, namely,
scaling the rows and columns of <I>A</I> in order to make the
condition number smaller. This is the case in the second
code fragment in the preceding section, where the program
chose to scale the rows by the factors returned in
<IMG WIDTH=191 HEIGHT=30 ALIGN=MIDDLE ALT="tex2html_wrap_inline18314" SRC="img595.gif">,
resulting in <IMG WIDTH=155 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18316" SRC="img596.gif">.
</BLOCKQUOTE>
<P>
As stated in section <A HREF="node138.html#seccomponentwise">6.4.2</A>,
this approach does not respect the presence
of zero or tiny entries in <I>A</I>. In contrast,
the ScaLAPACK computational routines
PxyyRFS (section <A HREF="node51.html#subseccomplineq">3.3.1</A>) or driver routines PxyySVX
(section <A HREF="node44.html#subsecdrivelineq">3.2.1</A>) will (except in rare cases)
compute a solution <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> with the following properties:
<P>
<BLOCKQUOTE> The componentwise backward error
of the computed solution <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> is the pair <I>E</I>,<I>f</I> which minimizes
<A NAME="5247"> </A>
<A NAME="5248"> </A>
<BR><IMG WIDTH=314 HEIGHT=48 ALIGN=BOTTOM ALT="displaymath18209" SRC="img597.gif"><BR>
(where we interpret 0/0 as 0)
subject to the constraint <IMG WIDTH=132 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18075" SRC="img546.gif">.
The minimal value of <IMG WIDTH=139 HEIGHT=40 ALIGN=MIDDLE ALT="tex2html_wrap_inline18330" SRC="img598.gif">
is given by
<BR><IMG WIDTH=351 HEIGHT=42 ALIGN=BOTTOM ALT="displaymath18210" SRC="img599.gif"><BR>
One can show that for most problems the <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> computed by PxyySVX
satisfies <IMG WIDTH=93 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18334" SRC="img600.gif">,
where <I>p</I>(<I>n</I>) is a modestly growing function of <I>n</I>.
In other words, <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> is the exact solution of the
perturbed problem <IMG WIDTH=132 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18075" SRC="img546.gif">,
where <I>E</I> and <I>f</I> are small relative perturbations in each entry of <I>A</I> and
<I>b</I>, respectively.
The corresponding condition number is
<IMG WIDTH=333 HEIGHT=31 ALIGN=MIDDLE ALT="tex2html_wrap_inline18352" SRC="img601.gif">.
<A NAME="5269"> </A>
The error <IMG WIDTH=40 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline18296" SRC="img589.gif"> is bounded by
<BR><IMG WIDTH=354 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18211" SRC="img602.gif"><BR>
<P>
The routines PxyyRFS and PxyySVX return <IMG WIDTH=17 HEIGHT=17 ALIGN=MIDDLE ALT="tex2html_wrap_inline18356" SRC="img603.gif">,
<A NAME="5274"> </A><A NAME="5275"> </A><A NAME="5276"> </A><A NAME="5277"> </A>
<A NAME="5278"> </A><A NAME="5279"> </A><A NAME="5280"> </A><A NAME="5281"> </A>
<A NAME="5282"> </A><A NAME="5283"> </A><A NAME="5284"> </A><A NAME="5285"> </A>
<A NAME="5286"> </A><A NAME="5287"> </A><A NAME="5288"> </A><A NAME="5289"> </A>
<A NAME="5290"> </A><A NAME="5291"> </A><A NAME="5292"> </A><A NAME="5293"> </A>
<A NAME="5294"> </A>
<A NAME="5295"> </A>
which is called <TT>BERR</TT><A NAME="5297"> </A>
(for Backward ERRor),
and a bound on the the actual error
<IMG WIDTH=120 HEIGHT=27 ALIGN=MIDDLE ALT="tex2html_wrap_inline18358" SRC="img604.gif">, called <TT>FERR</TT>
<A NAME="5301"> </A>
(for Forward ERRor), as
in the second code fragment in the last section.
<TT>FERR</TT> is actually calculated by the following formula, which can
be smaller than the bound <IMG WIDTH=104 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18360" SRC="img605.gif"> given above:
<BR><IMG WIDTH=464 HEIGHT=43 ALIGN=BOTTOM ALT="displaymath18212" SRC="img606.gif"><BR>
Here, <IMG WIDTH=8 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline18362" SRC="img607.gif"> is the computed value of the residual <IMG WIDTH=50 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline18364" SRC="img608.gif">, and
the norm in the numerator is estimated by using the same estimation
subroutine used for <TT>RCOND</TT>.
<P>
The value of
<TT>BERR</TT> for the example in the preceding section is <IMG WIDTH=68 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18366" SRC="img609.gif">.
<P>
Even in the rare cases where PxyyRFS fails to make
<TT>BERR</TT> close to its minimum <IMG WIDTH=6 HEIGHT=8 ALIGN=BOTTOM ALT="tex2html_wrap_inline17202" SRC="img429.gif">, the error bound <TT>FERR</TT>
may remain small. See [<A HREF="node189.html#ariolidemmelduff">9</A>]
for details.
<P>
</BLOCKQUOTE><HR><A NAME="tex2html3950" HREF="node140.html"><IMG WIDTH=37 HEIGHT=24 ALIGN=BOTTOM ALT="next" SRC="http://www.netlib.org/utk/icons/next_motif.gif"></A> <A NAME="tex2html3948" HREF="node132.html"><IMG WIDTH=26 HEIGHT=24 ALIGN=BOTTOM ALT="up" SRC="http://www.netlib.org/utk/icons/up_motif.gif"></A> <A NAME="tex2html3942" HREF="node138.html"><IMG WIDTH=63 HEIGHT=24 ALIGN=BOTTOM ALT="previous" SRC="http://www.netlib.org/utk/icons/previous_motif.gif"></A> <A NAME="tex2html3952" HREF="node1.html"><IMG WIDTH=65 HEIGHT=24 ALIGN=BOTTOM ALT="contents" SRC="http://www.netlib.org/utk/icons/contents_motif.gif"></A> <A NAME="tex2html3953" HREF="node190.html"><IMG WIDTH=43 HEIGHT=24 ALIGN=BOTTOM ALT="index" SRC="http://www.netlib.org/utk/icons/index_motif.gif"></A> <BR>
<B> Next:</B> <A NAME="tex2html3951" HREF="node140.html">Error Bounds for Linear </A>
<B>Up:</B> <A NAME="tex2html3949" HREF="node132.html">Accuracy and Stability</A>
<B> Previous:</B> <A NAME="tex2html3943" HREF="node138.html">Improved Error Bounds</A>
<P><ADDRESS>
<I>Susan Blackford <BR>
Tue May 13 09:21:01 EDT 1997</I>
</ADDRESS>
</BODY>
</HTML>
|
