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
|
<!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>Linear Equations</TITLE>
<META NAME="description" CONTENT="Linear Equations">
<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="tex2html2794" HREF="node52.html"><IMG WIDTH=37 HEIGHT=24 ALIGN=BOTTOM ALT="next" SRC="http://www.netlib.org/utk/icons/next_motif.gif"></A> <A NAME="tex2html2792" HREF="node50.html"><IMG WIDTH=26 HEIGHT=24 ALIGN=BOTTOM ALT="up" SRC="http://www.netlib.org/utk/icons/up_motif.gif"></A> <A NAME="tex2html2786" HREF="node50.html"><IMG WIDTH=63 HEIGHT=24 ALIGN=BOTTOM ALT="previous" SRC="http://www.netlib.org/utk/icons/previous_motif.gif"></A> <A NAME="tex2html2796" 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="tex2html2797" 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="tex2html2795" HREF="node52.html">Orthogonal Factorizations and Linear </A>
<B>Up:</B> <A NAME="tex2html2793" HREF="node50.html">Computational Routines</A>
<B> Previous:</B> <A NAME="tex2html2787" HREF="node50.html">Computational Routines</A>
<BR> <P>
<H2><A NAME="SECTION04331000000000000000">Linear Equations</A></H2>
<A NAME="subseccomplineq"> </A>
<P>
We use the standard notation for a system of simultaneous
linear<A NAME="1160"> </A>
equations<A NAME="1161"> </A>:
<BR><A NAME="Axb1"> </A><IMG WIDTH=500 HEIGHT=17 ALIGN=BOTTOM ALT="equation1162" SRC="img95.gif"><BR>
where <I>A</I> is the <B>coefficient matrix</B><A NAME="1166"> </A>,
<I>b</I> is the <B>right-hand side</B><A NAME="1168"> </A><A NAME="1169"> </A>,
and <I>x</I> is the <B>solution</B><A NAME="1171"> </A><A NAME="1172"> </A>.
In (<A HREF="node51.html#Axb1">3.2</A>) <I>A</I> is assumed to be a square matrix of order <I>n</I>,
but some of the individual routines allow <I>A</I> to be rectangular.
If there are several right-hand sides
we write
<BR><A NAME="AXB"> </A><IMG WIDTH=500 HEIGHT=17 ALIGN=BOTTOM ALT="equation1174" SRC="img96.gif"><BR>
where the columns of <I>B</I> are the individual right-hand sides,
and the columns of
<I>X</I> are the corresponding solutions.
The basic task is to compute <I>X</I>, given <I>A</I> and <I>B</I>.
<P>
If <I>A</I> is upper or lower triangular, (<A HREF="node51.html#Axb1">3.2</A>) can be solved by a
straightforward
process of backward or forward substitution.
Otherwise, the solution is obtained after first factorizing <I>A</I> as a
product of
triangular matrices (and possibly also a diagonal matrix or permutation
matrix).
<P>
The form of the factorization depends on the properties of the matrix
<I>A</I>.
ScaLAPACK provides routines for the following types of matrices, based on
the stated<A NAME="1178"> </A> factorizations:
<P>
<UL>
<LI> <B>general</B> matrices<A NAME="1181"> </A><A NAME="1182"> </A> (<I>LU</I> factorization with partial pivoting)<A NAME="1183"> </A><A NAME="1184"> </A>:
<BR><IMG WIDTH=289 HEIGHT=16 ALIGN=BOTTOM ALT="displaymath12933" SRC="img97.gif"><BR>
where <I>P</I> is a permutation matrix, <I>L</I> is lower triangular with unit
diagonal elements (lower trapezoidal if <I>m</I> > <I>n</I>), and U is upper
triangular (upper trapezoidal if <I>m</I> < <I>n</I>).
<LI> <B>symmetric and Hermitian positive definite</B> matrices
(Cholesky factorization)<A NAME="1186"> </A>:
<BR><IMG WIDTH=428 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12934" SRC="img98.gif"><BR>
<BR><IMG WIDTH=433 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12935" SRC="img99.gif"><BR>
where <I>U</I> is an upper triangular matrix and <I>L</I> is lower triangular.
<LI> <B>general band</B> matrices (<I>LU</I> factorization with partial pivoting):<BR>
<P>
If <I>A</I> is <I>m</I>-by-<I>n</I> with <I>bwl</I> subdiagonals and <I>bwu</I> superdiagonals,
the factorization is
<BR><IMG WIDTH=297 HEIGHT=16 ALIGN=BOTTOM ALT="displaymath12936" SRC="img100.gif"><BR>
where <I>P</I> and <I>Q</I> are permutation matrices and <I>L</I> and <I>U</I> are banded
lower and upper triangular matrices, respectively.
<LI> <B>general diagonally dominant-like band</B> matrices<A NAME="1197"> </A> including <B>general tridiagonal</B> matrices (<I>LU</I> factorization without pivoting):<BR>
<P>
A <B>diagonally dominant-like</B><A NAME="1200"> </A> matrix
is one for which it is known
<EM>a priori</EM> that pivoting for stability is NOT required in the <I>LU</I>
factorization of the matrix. Diagonally dominant matrices themselves are
examples of diagonally dominant-like matrices.<BR>
<P>
If <I>A</I> is <I>m</I>-by-<I>n</I> with <I>bwl</I> subdiagonals and <I>bwu</I> superdiagonals,
the factorization is
<BR><IMG WIDTH=302 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12937" SRC="img101.gif"><BR>
where <I>P</I> is a permutation matrix and <I>L</I> and <I>U</I> are banded lower
and upper triangular matrices respectively.
<LI> <B>symmetric and Hermitian positive definite band</B> matrices
(Cholesky factorization)<A NAME="1204"> </A>:
<BR><IMG WIDTH=466 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12938" SRC="img102.gif"><BR>
<BR><IMG WIDTH=471 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12939" SRC="img103.gif"><BR>
where <I>P</I> is a permutation matrix and <I>U</I> and <I>L</I> are banded upper and lower
triangular matrices, respectively.
<LI> <B>symmetric and Hermitian positive definite tridiagonal</B> matrices (<IMG WIDTH=48 HEIGHT=15 ALIGN=BOTTOM ALT="tex2html_wrap_inline13039" SRC="img104.gif"> factorization):<A NAME="1219"> </A><A NAME="1220"> </A>
<BR><IMG WIDTH=481 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12940" SRC="img106.gif"><BR>
<BR><IMG WIDTH=486 HEIGHT=21 ALIGN=BOTTOM ALT="displaymath12941" SRC="img107.gif"><BR>
where <I>P</I> is a permutation matrix and <I>U</I> and <I>L</I> are bidiagonal upper and
lower triangular matrices respectively.
<P>
</UL>
<P>
<B>Note</B>: In the banded and tridiagonal factorizations
(PxDBTRF, PxDTTRF, PxGBTRF, PxPBTRF, and PxPTTRF), the resulting
factorization is <EM>not</EM> the same factorization as returned from
LAPACK. Additional permutations are performed on the matrix for the
sake of parallelism.
Further details of the algorithmic implementations can be found
in [<A HREF="node189.html#lawn125">32</A>].
<P>
The factorization for a general diagonally dominant-like<A NAME="1237"> </A> tridiagonal
matrix is like that for
a general diagonally dominant-like band matrix with <I>bwl</I> = 1 and <I>bwu</I> = 1.
Band matrices use the band
storage scheme described in section <A HREF="node84.html#seclocalband">4.4.3</A>.
<P>
While the primary use of a matrix factorization is to solve a system
of equations, other related tasks are provided as well.
Wherever possible, ScaLAPACK provides routines to perform each of these
tasks
for each type of matrix and storage scheme (see
table <A HREF="node51.html#tabcomplineq1">3.6</A>).
The following list relates the tasks
to the last three characters of the name of the corresponding
computational routine:
<P>
<DL ><DT><STRONG>PxyyTRF:</STRONG>
<DD> factorize (obviously not needed for triangular matrices);
<A NAME="1241"> </A><A NAME="1242"> </A><A NAME="1243"> </A><A NAME="1244"> </A>
<A NAME="1245"> </A><A NAME="1246"> </A><A NAME="1247"> </A><A NAME="1248"> </A>
<A NAME="1249"> </A><A NAME="1250"> </A><A NAME="1251"> </A><A NAME="1252"> </A>
<A NAME="1253"> </A><A NAME="1254"> </A><A NAME="1255"> </A><A NAME="1256"> </A>
<A NAME="1257"> </A><A NAME="1258"> </A>
<A NAME="1259"> </A><A NAME="1260"> </A><A NAME="1261"> </A><A NAME="1262"> </A>
<A NAME="1263"> </A><A NAME="1264"> </A><A NAME="1265"> </A><A NAME="1266"> </A>
<P>
<DT><STRONG>PxyyTRS:</STRONG>
<DD> use the factorization (or the matrix <I>A</I> itself if it is
triangular) to
solve (<A HREF="node51.html#AXB">3.3</A>) by forward or backward substitution;
<A NAME="1268"> </A><A NAME="1269"> </A><A NAME="1270"> </A><A NAME="1271"> </A>
<A NAME="1272"> </A><A NAME="1273"> </A><A NAME="1274"> </A><A NAME="1275"> </A>
<A NAME="1276"> </A><A NAME="1277"> </A><A NAME="1278"> </A><A NAME="1279"> </A>
<A NAME="1280"> </A><A NAME="1281"> </A><A NAME="1282"> </A><A NAME="1283"> </A>
<A NAME="1284"> </A><A NAME="1285"> </A>
<A NAME="1286"> </A><A NAME="1287"> </A><A NAME="1288"> </A><A NAME="1289"> </A>
<A NAME="1290"> </A><A NAME="1291"> </A><A NAME="1292"> </A><A NAME="1293"> </A>
<A NAME="1294"> </A><A NAME="1295"> </A><A NAME="1296"> </A><A NAME="1297"> </A>
<P>
<DT><STRONG>PxyyCON:</STRONG>
<DD> estimate the reciprocal of the condition number
<IMG WIDTH=141 HEIGHT=30 ALIGN=MIDDLE ALT="tex2html_wrap_inline13057" SRC="img108.gif">;
Higham's modification [<A HREF="node189.html#nick2">81</A>] of Hager's method [<A HREF="node189.html#hager">72</A>]
is used to estimate <IMG WIDTH=45 HEIGHT=30 ALIGN=MIDDLE ALT="tex2html_wrap_inline13059" SRC="img109.gif"> (not provided for band or
tridiagonal matrices);
<A NAME="1302"> </A><A NAME="1303"> </A><A NAME="1304"> </A><A NAME="1305"> </A>
<A NAME="1306"> </A><A NAME="1307"> </A><A NAME="1308"> </A><A NAME="1309"> </A>
<A NAME="1310"> </A><A NAME="1311"> </A><A NAME="1312"> </A><A NAME="1313"> </A>
<P>
<DT><STRONG>PxyyRFS:</STRONG>
<DD> compute bounds on the error in the computed solution (returned
by the PxyyTRS routine), and
refine the solution to reduce the backward error (see below) (not
provided for band or tridiagonal matrices);
<A NAME="1314"> </A><A NAME="1315"> </A><A NAME="1316"> </A><A NAME="1317"> </A>
<A NAME="1318"> </A><A NAME="1319"> </A><A NAME="1320"> </A><A NAME="1321"> </A>
<A NAME="1322"> </A><A NAME="1323"> </A><A NAME="1324"> </A><A NAME="1325"> </A>
<P>
<DT><STRONG>PxyyTRI:</STRONG>
<DD> use the factorization (or the matrix <I>A</I> itself if it is
triangular)
to compute <IMG WIDTH=29 HEIGHT=15 ALIGN=BOTTOM ALT="tex2html_wrap_inline13063" SRC="img110.gif"> (not provided for band matrices, because the inverse
does not in general preserve bandedness);
<A NAME="1327"> </A><A NAME="1328"> </A><A NAME="1329"> </A><A NAME="1330"> </A>
<A NAME="1331"> </A><A NAME="1332"> </A><A NAME="1333"> </A><A NAME="1334"> </A>
<A NAME="1335"> </A><A NAME="1336"> </A><A NAME="1337"> </A><A NAME="1338"> </A>
<A NAME="1339"> </A>
<A NAME="1340"> </A>
<P>
<DT><STRONG>PxyyEQU:</STRONG>
<DD> compute scaling factors to equilibrate<A NAME="1341"> </A> <I>A</I>
(not provided for band, tridiagonal, or triangular matrices). These
routines do not actually scale
the matrices: auxiliary routines PxLAQyy may be used for that purpose --
see the code of the driver routines PxyySVX for sample usage<A NAME="1342"> </A>.
<A NAME="1343"> </A><A NAME="1344"> </A><A NAME="1345"> </A><A NAME="1346"> </A>
<A NAME="1347"> </A><A NAME="1348"> </A><A NAME="1349"> </A><A NAME="1350"> </A>
<P>
</DL>
<P>
Note that some of the above routines depend on the output of others:
<P>
<DL ><DT><STRONG>PxyyTRF:</STRONG>
<DD> may work on an equilibrated matrix produced by
PxyyEQU and PxLAQyy, if yy is one of {GE, PO};
<P>
<DT><STRONG>PxyyTRS:</STRONG>
<DD> requires the factorization returned by PxyyTRF;
<P>
<DT><STRONG>PxyyCON:</STRONG>
<DD> requires the norm of the original matrix <I>A</I> and the
factorization returned by PxyyTRF;
<P>
<DT><STRONG>PxyyRFS:</STRONG>
<DD> requires the original matrices <I>A</I> and <I>B</I>, the factorization
returned by PxyyTRF, and the solution <I>X</I> returned by PxyyTRS;
<P>
<DT><STRONG>PxyyTRI:</STRONG>
<DD> requires the factorization returned by PxyyTRF.
<P>
</DL>
<P>
The RFS (``refine solution'') routines perform iterative
refinement<A NAME="1354"> </A>
and compute backward and forward error<A NAME="1355"> </A> bounds for the solution.
Iterative refinement is done in the same precision as the input data.
In particular, the residual is <EM>not</EM> computed with extra precision,
as has been traditionally done.
The benefit of this procedure is discussed in section <A HREF="node139.html#secAxb">6.5</A>.
<P>
<P><A NAME="1359"> </A><A NAME="tabcomplineq1"> </A><IMG WIDTH=739 HEIGHT=720 ALIGN=BOTTOM ALT="table1358" SRC="img111.gif"><BR>
<STRONG>Table 3.6:</STRONG> Computational routines for linear equations<BR>
<P>
<P>
<HR><A NAME="tex2html2794" HREF="node52.html"><IMG WIDTH=37 HEIGHT=24 ALIGN=BOTTOM ALT="next" SRC="http://www.netlib.org/utk/icons/next_motif.gif"></A> <A NAME="tex2html2792" HREF="node50.html"><IMG WIDTH=26 HEIGHT=24 ALIGN=BOTTOM ALT="up" SRC="http://www.netlib.org/utk/icons/up_motif.gif"></A> <A NAME="tex2html2786" HREF="node50.html"><IMG WIDTH=63 HEIGHT=24 ALIGN=BOTTOM ALT="previous" SRC="http://www.netlib.org/utk/icons/previous_motif.gif"></A> <A NAME="tex2html2796" 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="tex2html2797" 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="tex2html2795" HREF="node52.html">Orthogonal Factorizations and Linear </A>
<B>Up:</B> <A NAME="tex2html2793" HREF="node50.html">Computational Routines</A>
<B> Previous:</B> <A NAME="tex2html2787" HREF="node50.html">Computational Routines</A>
<P><ADDRESS>
<I>Susan Blackford <BR>
Tue May 13 09:21:01 EDT 1997</I>
</ADDRESS>
</BODY>
</HTML>
|
