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
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
|
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.6"/>
<title>ViennaCL - The Vienna Computing Library: Algorithms</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<script type="text/javascript">
$(document).ready(initResizable);
$(window).load(resizeHeight);
</script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
$(document).ready(function() { searchBox.OnSelectItem(0); });
</script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<td style="padding-left: 0.5em;">
<div id="projectname">ViennaCL - The Vienna Computing Library
 <span id="projectnumber">1.7.1</span>
</div>
<div id="projectbrief">Free open-source GPU-accelerated linear algebra and solver library.</div>
</td>
<td> <div id="MSearchBox" class="MSearchBoxInactive">
<span class="left">
<img id="MSearchSelect" src="search/mag_sel.png"
onmouseover="return searchBox.OnSearchSelectShow()"
onmouseout="return searchBox.OnSearchSelectHide()"
alt=""/>
<input type="text" id="MSearchField" value="Search" accesskey="S"
onfocus="searchBox.OnSearchFieldFocus(true)"
onblur="searchBox.OnSearchFieldFocus(false)"
onkeyup="searchBox.OnSearchFieldChange(event)"/>
</span><span class="right">
<a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.6 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
</div><!-- top -->
<div id="side-nav" class="ui-resizable side-nav-resizable">
<div id="nav-tree">
<div id="nav-tree-contents">
<div id="nav-sync" class="sync"></div>
</div>
</div>
<div id="splitbar" style="-moz-user-select:none;"
class="ui-resizable-handle">
</div>
</div>
<script type="text/javascript">
$(document).ready(function(){initNavTree('manual-algorithms.html','');});
</script>
<div id="doc-content">
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
onmouseover="return searchBox.OnSearchSelectShow()"
onmouseout="return searchBox.OnSearchSelectHide()"
onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark"> </span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark"> </span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark"> </span>Namespaces</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark"> </span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark"> </span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark"> </span>Variables</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(6)"><span class="SelectionMark"> </span>Typedefs</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(7)"><span class="SelectionMark"> </span>Enumerations</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(8)"><span class="SelectionMark"> </span>Enumerator</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(9)"><span class="SelectionMark"> </span>Friends</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(10)"><span class="SelectionMark"> </span>Macros</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(11)"><span class="SelectionMark"> </span>Pages</a></div>
<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0"
name="MSearchResults" id="MSearchResults">
</iframe>
</div>
<div class="header">
<div class="headertitle">
<div class="title">Algorithms </div> </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><p>This chapter gives an overview over the available algorithms in ViennaCL. The focus of ViennaCL is on iterative solvers, for which generic implementations that allows the use of the same code on the CPU (either using Boost.uBLAS, Eigen, MTL4, or ViennaCL types) and on the GPU (using ViennaCL types) are provided.</p>
<h1><a class="anchor" id="manual-algorithms-direct-solvers"></a>
Direct Solvers</h1>
<p>ViennaCL provides triangular solvers and LU factorization without pivoting for the solution of dense linear systems. The interface is similar to that of Boost.uBLAS. </p>
<div class="fragment"><div class="line"><span class="keyword">using namespace </span>viennacl::linalg; <span class="comment">//to keep solver calls short</span></div>
<div class="line"><a class="code" href="classviennacl_1_1matrix.html">viennacl::matrix<float></a> vcl_matrix;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<float></a> vcl_rhs;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<float></a> vcl_result;</div>
<div class="line"></div>
<div class="line"><span class="comment">// Set up matrix and vectors here</span></div>
<div class="line"></div>
<div class="line"><span class="comment">// solution of an upper triangular system:</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">solve</a>(vcl_matrix, vcl_rhs, <a class="code" href="structviennacl_1_1linalg_1_1upper__tag.html">upper_tag</a>());</div>
<div class="line"><span class="comment">//solution of a lower triangular system:</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">solve</a>(vcl_matrix, vcl_rhs, <a class="code" href="structviennacl_1_1linalg_1_1lower__tag.html">lower_tag</a>());</div>
<div class="line"></div>
<div class="line"><span class="comment">// solution of a full system right into the load vector vcl_rhs:</span></div>
<div class="line"><a class="code" href="namespaceviennacl_1_1linalg.html#a1050039fb153ae8ad818390a30f96bc4">lu_factorize</a>(vcl_matrix);</div>
<div class="line"><a class="code" href="namespaceviennacl_1_1linalg.html#a118a9164bc965f8ec728eda73a5e2450">lu_substitute</a>(vcl_matrix, vcl_rhs);</div>
</div><!-- fragment --><p> In ViennaCL there is no pivoting included in the LU factorization process, hence the computation may break down or yield results with poor accuracy. However, for certain classes of matrices (like diagonal dominant matrices) good results can be obtained without pivoting.</p>
<p>It is also possible to solve for multiple right hand sides: </p>
<div class="fragment"><div class="line"><span class="keyword">using namespace </span>viennacl::linalg; <span class="comment">//to keep solver calls short</span></div>
<div class="line"><a class="code" href="classviennacl_1_1matrix.html">viennacl::matrix<float></a> vcl_matrix;</div>
<div class="line"><a class="code" href="classviennacl_1_1matrix.html">viennacl::matrix<float></a> vcl_rhs_matrix;</div>
<div class="line"><a class="code" href="classviennacl_1_1matrix.html">viennacl::matrix<float></a> vcl_result;</div>
<div class="line"></div>
<div class="line"><span class="comment">// Set up matrices here</span></div>
<div class="line"></div>
<div class="line"><span class="comment">// solution of an upper triangular system:</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">solve</a>(vcl_matrix, vcl_rhs_matrix, <a class="code" href="structviennacl_1_1linalg_1_1upper__tag.html">upper_tag</a>());</div>
<div class="line"></div>
<div class="line"><span class="comment">// solution of a lower triangular system:</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">solve</a>(vcl_matrix, vcl_rhs_matrix, <a class="code" href="structviennacl_1_1linalg_1_1lower__tag.html">lower_tag</a>());</div>
</div><!-- fragment --><h1><a class="anchor" id="manual-algorithms-iterative-solvers"></a>
Iterative Solvers</h1>
<p>Iterative solvers approximately solve a (usually sparse) system <img class="formulaInl" alt="$ Ax = b $" src="form_127.png"/> through iterated application of the matrix <img class="formulaInl" alt="$ A $" src="form_1.png"/> to vectors. ViennaCL provides the following iterative solvers for use with both assembled matrices as well as in a so-called matrix-free setting, where to user provides the operator <img class="formulaInl" alt="$ A $" src="form_1.png"/>:</p>
<center> <table class="doxtable">
<tr>
<th>Method </th><th>Matrix class </th><th>ViennaCL </th></tr>
<tr>
<td>Conjugate Gradient (CG) </td><td>symmetric positive definite </td><td><code>y = solve(A, x, cg_tag());</code> </td></tr>
<tr>
<td>Mixed-Precision Conjugate Gradient (Mixed-CG) </td><td>symmetric positive definite </td><td><code>y = solve(A, x, mixed_precision_cg_tag());</code> </td></tr>
<tr>
<td>Stabilized Bi-CG (BiCGStab) </td><td>non-symmetric </td><td><code>y = solve(A, x, bicgstab_tag());</code> </td></tr>
<tr>
<td>Generalized Minimum Residual (GMRES) </td><td>general </td><td><code>y = solve(A, x, gmres_tag());</code> </td></tr>
</table>
<p><b>Linear solver routines in ViennaCL for the computation of <img class="formulaInl" alt="$ x $" src="form_22.png"/> in the expression <img class="formulaInl" alt="$ Ax = b $" src="form_127.png"/> with given <img class="formulaInl" alt="$ A $" src="form_1.png"/>, <img class="formulaInl" alt="$ b $" src="form_128.png"/>.</b> </center><p> Pipelined versions of CG, BiCGStab as well as GMRES are implemented for the case that no preconditioner is provided. This provides performance benefits for medium-sized systems of about 10k to 100k unknowns, because kernel launch and data transfer latencies are reduced by a factor of two to three.</p>
<p>Unlike direct solvers, the convergence of iterative solvers relies on certain properties of the system matrix. Keep in mind that an iterative solver may fail to converge, especially if the matrix is ill conditioned or a wrong solver is chosen.</p>
<dl class="section note"><dt>Note</dt><dd>The iterative solvers can also be used for Armadillo, Boost.uBLAS, Eigen and MTL4 objects! Have a look at <a class="el" href="manual-interfacing.html">Interfacing Other Libraries</a> and the respective tutorials in the <code>examples/tutorials/</code> folder.</dd>
<dd>
An example of a matrix-free use of the iterative solvers can be found in <a class="el" href="matrix-free_8cpp.html">examples/tutorial/matrix-free.cpp</a>, where additional informations on the interface requirements are given.</dd></dl>
<h2><a class="anchor" id="manual-algorithms-iterative-solvers-cg"></a>
Conjugate Gradients (CG)</h2>
<p>The conjugate gradient method is the method of choice for many symmetric positive definite systems. A minimal snippet for running a CG solver is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1compressed__matrix.html">viennacl::compressed_matrix<T></a> A;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> b;</div>
<div class="line"></div>
<div class="line"><span class="comment">// Set up A and b here</span></div>
<div class="line"></div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> x = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b, <a class="code" href="classviennacl_1_1linalg_1_1cg__tag.html">viennacl::linalg::cg_tag</a>());</div>
</div><!-- fragment --><p> A preconditioner object can be passed as optional fourth argument to the <a class="el" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">solve()</a> routine. Customized iteration counts as well as relative or absolute tolerances can be provided through the respective member functions of <code>cg_tag</code>: The convention is that solver tags take the relative error tolerance as first argument and the maximum number of iteration steps as second argument. Furthermore, after the solver run the number of iterations and the estimated error can be obtained from the solver tags as follows: </p>
<div class="fragment"><div class="line"><span class="comment">// conjugate gradient solver with tolerance 1e10</span></div>
<div class="line"><span class="comment">// and at most 100 iterations:</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1cg__tag.html">viennacl::linalg::cg_tag</a> custom_cg(1e-10, 100);</div>
<div class="line">x = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b, custom_cg);</div>
<div class="line"></div>
<div class="line"><span class="comment">// print number of iterations taken and estimated error:</span></div>
<div class="line">std::cout << <span class="stringliteral">"No. of iters: "</span> << custom_cg.iters() << std::endl;</div>
<div class="line">std::cout << <span class="stringliteral">"Est. error: "</span> << custom_cg.error() << std::endl;</div>
</div><!-- fragment --><p> Absolute tolerances can be set via the member function <code>abs_tolerance()</code>.</p>
<p>An initial guess <img class="formulaInl" alt="$ x_0 $" src="form_129.png"/> for the iterative solution of the system <img class="formulaInl" alt="$ Ax = b $" src="form_127.png"/> can be incorporated manually by considering </p>
<p class="formulaDsp">
<img class="formulaDsp" alt="\[ x = x_0 + y \\ \Rightarrow A(x_0 + y) = b \Leftrightarrow Ay = b - Ax_0 := \tilde{b} \]" src="form_130.png"/>
</p>
<p> Thus, one may provide a nonzero initial guess by supplying the modified right hand side <img class="formulaInl" alt="$ b - Ax_0 $" src="form_131.png"/> instead of <img class="formulaInl" alt="$ b $" src="form_128.png"/>. The obtained solution <img class="formulaInl" alt="$ y $" src="form_20.png"/> then needs to be added to <img class="formulaInl" alt="$ x_0 $" src="form_129.png"/> to obtain <img class="formulaInl" alt="$ x $" src="form_22.png"/>.</p>
<p>If a custom monitor callback function should be provided, or ViennaCL should deal with the initial guess internally, then the extended CG interface is required. This extended interface relies on the solver class <code>cg_solver<T></code>, which takes the vector type as template argument and the <code>cg_tag</code> as constructor argument. A custom monitor callback routine returns a bool and takes the following three parameters:</p>
<ul>
<li>The current approximation to the solution <img class="formulaInl" alt="$ x $" src="form_22.png"/></li>
<li>The current relative residual norm estimate</li>
<li>An optional user-provided void pointer to additional user data If the monitor returns <code>true</code>, the iterative solver terminates. If the monitor returns <code>false</code>, the iterative solver may still terminate because the maximum iteration count, or the absolute or relative tolerance is reached.</li>
</ul>
<p>An exemplary code snippet for a CG solver with custom monitor callback and initial guess is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1compressed__matrix.html">viennacl::compressed_matrix<T></a> A;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> b;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> init_guess;</div>
<div class="line"></div>
<div class="line"><span class="comment">// Set up A, b, and init_guess here</span></div>
<div class="line"></div>
<div class="line"><span class="comment">// Set up CG solver object</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1cg__tag.html">viennacl::linalg::cg_tag</a> my_cg_tag(1e-5, 20);</div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1cg__solver.html">viennacl::linalg::cg_solver<viennacl::vector<T></a> > my_cg_solver(my_cg_tag);</div>
<div class="line"></div>
<div class="line"><span class="comment">// Register monitor and initial guess:</span></div>
<div class="line">my_cg_solver.set_monitor(my_custom_monitor, &my_monitor_data);</div>
<div class="line">my_cg_solver.set_initial_guess(init_guess);</div>
<div class="line"></div>
<div class="line"><span class="comment">// Run solver:</span></div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> x = my_cg_solver(A, b);</div>
</div><!-- fragment --><dl class="section note"><dt>Note</dt><dd>The monitor routine is called for the CG method applied to the system <img class="formulaInl" alt="$ Ay = b - Ax_0 $" src="form_132.png"/> if an initial guess is provided.</dd></dl>
<h2><a class="anchor" id="manual-algorithms-iterative-solvers-mixed-cg"></a>
Mixed-Precision Conjugate Gradients</h2>
<p>A mixed precision CG solver is available for symmetric positive systems with low to moderate condition number in double precision. The solver operates in a two-stage manner: The residual is computed in double precision, for which a correction is then computed in low precision.</p>
<p>While single precision floating point numbers only require half the number of bytes as double precision floating point numbers, in practice one also needs to transfer indices. This limits the possible performance gains over a conventional CG method to about 20 percent in practice. Also, a mixed precision solver is fairly sensitive with respect to high condition numbers, hence we advise users to only use a mixed precision CG if the system matrix is known to be well-behaved.</p>
<p>Sample code for running a mixed-precision CG solver is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1linalg_1_1mixed__precision__cg__tag.html">viennacl::linalg::mixed_precision_cg_tag</a> mixed_prec_cg_config;</div>
<div class="line">x = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b, mixed_prec_cg_config);</div>
</div><!-- fragment --><p> As usual, the first parameter to the constructor of <code>mixed_precision_cg_tag</code> is the relative tolerance for the residual, while the second parameter denotes the maximum number of solver iterations. The third parameter denotes the relative tolerance for the inner low-precision CG iterations and defaults to <code>0.01</code>.</p>
<dl class="section note"><dt>Note</dt><dd>Have a look at <code>examples/banchmarks/solver.cpp</code> for an example.</dd>
<dd>
A mixed-precision solver makes sense only if the matrix and right-hand-side vector are supplied in <code>double</code> precision.</dd></dl>
<p>Currently no extended interface for passing monitors or initial guesses is available for the mixed precision CG solver.</p>
<h2><a class="anchor" id="manual-algorithms-iterative-solvers-bicgstab"></a>
Stabilized Bi-CG (BiCGStab)</h2>
<p>The BiCGStab method is an attractive option for non-symmetric systems. It can be used for general systems, but it is not guaranteed to converge. In comparison to GMRES, BiCGStab has constant costs per iteration and lower memory footprint.</p>
<p>In its shortest form, a BiCGStab solver is called as follows: </p>
<div class="fragment"><div class="line">x = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b, <a class="code" href="classviennacl_1_1linalg_1_1bicgstab__tag.html">viennacl::linalg::bicgstab_tag</a>());</div>
</div><!-- fragment --><p> Similar to the CG method, relative tolerances and maximum iteration counts can be passed to the constructor of the <code>bicgstab_tag</code> object.</p>
<p>An extended interface for custom monitors and initial guesses is available. The templated class <code><a class="el" href="classviennacl_1_1linalg_1_1bicgstab__solver.html">viennacl::linalg::bicgstab_solver</a><VectorT></code> is used in exactly the same way as <code>cg_solver</code> above.</p>
<h2><a class="anchor" id="manual-algorithms-iterative-solvers-gmres"></a>
Generalized Minimum Residual (GMRES)</h2>
<p>ViennaCL provides an implementation of the GMRES method with (optional) restart. GMRES is guaranteed to converge after <img class="formulaInl" alt="$ n $" src="form_133.png"/> steps for a system with <img class="formulaInl" alt="$ n $" src="form_133.png"/> unknowns, even though in practice one usually needs significantly less iterations. The computational cost of GMRES grows with the number of iterations, hence GMRES is either used with a restart after every <img class="formulaInl" alt="$ k $" src="form_134.png"/> steps, or with a good preconditioner.</p>
<p>A sample code snippet for running the GMRES solver in ViennaCL is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1compressed__matrix.html">viennacl::compressed_matrix<T></a> A;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> b;</div>
<div class="line"></div>
<div class="line"><span class="comment">// Set up A and b here</span></div>
<div class="line"></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1gmres__tag.html">viennacl::linalg::gmres_tag</a> my_gmres_tag(1e-5, 100, 20); <span class="comment">// up to 100 iterations, restart after 20 iterations</span></div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<T></a> x = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b, my_gmres_tag);</div>
<div class="line"></div>
<div class="line"><span class="comment">// print number of iterations taken and estimated error:</span></div>
<div class="line">std::cout << <span class="stringliteral">"No. of iters: "</span> << my_gmres_tag.iters() << std::endl;</div>
<div class="line">std::cout << <span class="stringliteral">"Est. error: "</span> << my_gmres_tag.error() << std::endl;</div>
</div><!-- fragment --><h1><a class="anchor" id="manual-algorithms-preconditioners"></a>
Preconditioners</h1>
<p>ViennaCL provides (partially) generic implementations of several preconditioners. Due to the need to dynamically allocate memory, preconditioner setup is usually carried out on the CPU host. Thus, one may not obtain an overall performance benefit if too much time is spent on the preconditioner setup. Exceptions are simple diagonal preconditioners and the Chow-Patel ILU and ICC preconditioners, which have static allocation patterns and can be computed on the backend device.</p>
<dl class="section note"><dt>Note</dt><dd>Some preconditioners also work for Boost.uBLAS types!</dd></dl>
<p>An overview of preconditioners available for the various sparse matrix types is as follows: </p>
<center> <table class="doxtable">
<tr>
<th>Matrix Type </th><th>ICHOL </th><th>(Block-)ILU[0/T] </th><th>Jacobi </th><th>Row-scaling </th><th>AMG </th><th>SPAI </th></tr>
<tr>
<td><code>compressed_matrix</code> </td><td>yes </td><td>yes </td><td>yes </td><td>yes </td><td>yes </td><td>yes </td></tr>
<tr>
<td><code>coordinate_matrix</code> </td><td>no </td><td>no </td><td>yes </td><td>yes </td><td>no </td><td>no </td></tr>
<tr>
<td><code>ell_matrix</code> </td><td>no </td><td>no </td><td>no </td><td>no </td><td>no </td><td>no </td></tr>
<tr>
<td><code>hyb_matrix</code> </td><td>no </td><td>no </td><td>no </td><td>no </td><td>no </td><td>no </td></tr>
<tr>
<td><code>sliced_ell_matrix</code> </td><td>no </td><td>no </td><td>no </td><td>no </td><td>no </td><td>no </td></tr>
</table>
</center><p> We aim to provide broader support for preconditioners using other sparse matrix formats in future releases. Sparse approximate inverse (SPAI) preconditioners are described in <a class="el" href="manual-additional-algorithms.html">Additional Algorithms</a> section.</p>
<h2><a class="anchor" id="manual-algorithms-preconditioners-parallel-ilu0"></a>
Parallel Incomplete LU Factorization with Static Pattern (Chow-Patel-ILU0)</h2>
<p>Incomplete LU (ILU) factorizations are popular black-box preconditioners and may work in cases where more advanced and problem-specific techniques such as multigrid approaches are not possible. A drawback of classical ILU algorithms is the sequential nature, which prohibits good performance on massively parallel hardware such as GPUs.</p>
<p>ViennaCL provides an implementation of the recently proposed parallel ILU factorization algorithm proposed by Chow and Patel <a class="el" href="citelist.html#CITEREF_chow:fine-grained-ilu">[7]</a>. While the authors propose an asynchronous algorithm, we use a synchronous variant for reasons of better maintainability and user support in the case of failures. Internal evaluations suggest that the performance difference of the synchronous and asynchronous variants are negligible in practice. Also, the static pattern of <img class="formulaInl" alt="$ L $" src="form_24.png"/> and <img class="formulaInl" alt="$ U $" src="form_25.png"/> are taken from the system matrix <img class="formulaInl" alt="$ A $" src="form_1.png"/> for efficiency reasons.</p>
<p>The preconditioner exposes two parameters: First, the incomplete factors <img class="formulaInl" alt="$ L $" src="form_24.png"/> and <img class="formulaInl" alt="$ U $" src="form_25.png"/> with <img class="formulaInl" alt="$ A \approx LU $" src="form_26.png"/> are computed in parallel using a nonlinear iteration scheme. The number of these sweeps can be adjusted in the configuration class <code>chow_patel_tag</code>. Second, the forward- and backward solves in the preconditioner application is approximated by a Jacobi method for the truncated Neumann series. That is, instead of the forward- or backward substition when solving <img class="formulaInl" alt="$ R x = r $" src="form_135.png"/> for some vector <img class="formulaInl" alt="$ r $" src="form_136.png"/> and lower or upper triangular matrix <img class="formulaInl" alt="$ R $" src="form_37.png"/>, one or several Jacobi iterations of the form </p>
<p class="formulaDsp">
<img class="formulaDsp" alt="\[ x_{k+1} = (I - R)x_k + D^{-1} b \]" src="form_137.png"/>
</p>
<p> are used. The number of Jacobi iterations for each triangular factor is a second parameter. For full details of the underlying algorithm we refer to the original paper <a class="el" href="citelist.html#CITEREF_chow:fine-grained-ilu">[7]</a>.</p>
<p>The preconditioner class is <code>chow_patel_ilu_precond<MatrixT></code> with sparse matrix type <code>MatrixT</code>. <code>MatrixT</code> is currently limited to <code>compressed_matrix</code>, which denotes the system matrix type passed as first parameter. A sample code snippet for using the Chow-Patel-ILU0 preconditioner with a BiCGStab solver is as follows: </p>
<div class="fragment"><div class="line"><span class="comment">// configuration of preconditioner:</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1chow__patel__tag.html">viennacl::linalg::chow_patel_tag</a> chow_patel_ilu_config;</div>
<div class="line">chow_patel_ilu_config.<a class="code" href="classviennacl_1_1linalg_1_1chow__patel__tag.html#a646bc4b6df30814e9f64d96d9b872738">sweeps</a>(3); <span class="comment">// three nonlinear sweeps</span></div>
<div class="line">chow_patel_ilu_config.<a class="code" href="classviennacl_1_1linalg_1_1chow__patel__tag.html#a57b522b1e889fe09b1e5b9ed47b7e036">jacobi_iters</a>(2); <span class="comment">// two Jacobi iterations per triangular 'solve' Rx=r</span></div>
<div class="line"></div>
<div class="line"><span class="comment">// create and compute preconditioner:</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1chow__patel__ilu__precond.html">viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix<ScalarType></a> > chow_patel_ilu(A, chow_patel_ilu_config);</div>
<div class="line"></div>
<div class="line"><span class="comment">// use in solver (e.g. using BiCGStab solver)</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b,</div>
<div class="line"> <a class="code" href="classviennacl_1_1linalg_1_1bicgstab__tag.html">viennacl::linalg::bicgstab_tag</a>(), <span class="comment">// solver here</span></div>
<div class="line"> chow_patel_ilu); <span class="comment">// preconditioner here</span></div>
</div><!-- fragment --><p>Even though the setup and solve phase are executed on the GPU when using the CUDA or OpenCL backend, any overall performance gains are highly problem-specific. The number of nonlinear sweeps and Jacobi iterations need to be set problem-specific for best performance. Values between one and four are likely to give best results.</p>
<h2><a class="anchor" id="manual-algorithms-preconditioners-parallel-icc0"></a>
Parallel Incomplete Cholesky Factorization with Static Pattern (Chow-Patel-IChol0)</h2>
<p>The parallel incomplete Cholesky factorization is similar to the parallel Chow-Patel-ILU0 preconditioner in the previous section, but with a factorization $LL^T$ for symmetric positive definite matrices. The associated preconditioner class is <code>chow_patel_ichol_precond<MatrixT></code> and used in the same way as the class chow_patel_ilu_precond` described above: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1linalg_1_1chow__patel__icc__precond.html">viennacl::linalg::chow_patel_icc_precond< viennacl::compressed_matrix<ScalarType></a> > vcl_chow_patel_icc(A, <a class="code" href="classviennacl_1_1linalg_1_1chow__patel__tag.html">viennacl::linalg::chow_patel_tag</a>());</div>
</div><!-- fragment --><p> The configuration tag of type <code>chow_patel_tag</code> is of the same type for both the incomplete LU and the incomplete Cholesky factorization case. Thus, the same parameters are exposed for tweaking.</p>
<h2><a class="anchor" id="manual-algorithms-preconditioners-ilut"></a>
Incomplete LU Factorization with Threshold (ILUT)</h2>
<p>The incomplete LU factorization preconditioner aims at computing sparse matrices lower and upper triangular matrices <img class="formulaInl" alt="$ L $" src="form_24.png"/> and <img class="formulaInl" alt="$ U $" src="form_25.png"/> such that the sparse system matrix is approximately given by <img class="formulaInl" alt="$ A \approx LU $" src="form_26.png"/>. In order to control the sparsity pattern of <img class="formulaInl" alt="$ L $" src="form_24.png"/> and <img class="formulaInl" alt="$ U $" src="form_25.png"/>, a threshold strategy is used (ILUT) <a class="el" href="citelist.html#CITEREF_saad-iterative-solution">[26]</a> . Due to the serial nature of the preconditioner, the setup of ILUT is always computed on the CPU using the respective ViennaCL backend.</p>
<div class="fragment"><div class="line"><span class="comment">// compute ILUT preconditioner:</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1ilut__tag.html">viennacl::linalg::ilut_tag</a> ilut_config;</div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1ilut__precond.html">viennacl::linalg::ilut_precond< SparseMatrix ></a> vcl_ilut(vcl_matrix, ilut_config);</div>
<div class="line"></div>
<div class="line"><span class="comment">// solve (e.g. using conjugate gradient solver)</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(vcl_matrix, vcl_rhs,</div>
<div class="line"> <a class="code" href="classviennacl_1_1linalg_1_1bicgstab__tag.html">viennacl::linalg::bicgstab_tag</a>(), <span class="comment">// solver here</span></div>
<div class="line"> vcl_ilut); <span class="comment">// preconditioner here</span></div>
</div><!-- fragment --><p> The triangular substitutions may be applied in parallel on GPUs by enabling <em>level-scheduling</em> <a class="el" href="citelist.html#CITEREF_saad-iterative-solution">[26]</a> via the member function call <code>use_level_scheduling(true)</code> in the <code>ilut_config</code> object.</p>
<p>Three parameters can be passed to the constructor of <code>ilut_tag</code>: The first specifies the maximum number of entries per row in <img class="formulaInl" alt="$ L $" src="form_24.png"/> and <img class="formulaInl" alt="$ U $" src="form_25.png"/>, while the second parameter specifies the drop tolerance. The third parameter is the boolean specifying whether level scheduling should be used.</p>
<dl class="section note"><dt>Note</dt><dd>The performance of level scheduling depends strongly on the matrix pattern and is thus disabled by default.</dd></dl>
<h2><a class="anchor" id="manual-algorithms-preconditioners-ilu0"></a>
Incomplete LU Factorization with Static Pattern (ILU0)</h2>
<p>Similar to ILUT, ILU0 computes an approximate LU factorization with sparse factors L and U. While ILUT determines the location of nonzero entries on the fly, ILU0 uses the sparsity pattern of A for the sparsity pattern of L and U <a class="el" href="citelist.html#CITEREF_saad-iterative-solution">[26]</a> Due to the serial nature of the preconditioner, the setup of ILU0 is computed on the CPU. </p>
<div class="fragment"><div class="line"><span class="comment">// compute ILU0 preconditioner:</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1ilu0__tag.html">viennacl::linalg::ilu0_tag</a> ilu0_config;</div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1ilu0__precond.html">viennacl::linalg::ilu0_precond< SparseMatrix ></a> vcl_ilut(vcl_matrix, ilu0_config);</div>
<div class="line"></div>
<div class="line"><span class="comment">// solve (e.g. using conjugate gradient solver)</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(vcl_matrix, vcl_rhs,</div>
<div class="line"> <a class="code" href="classviennacl_1_1linalg_1_1bicgstab__tag.html">viennacl::linalg::bicgstab_tag</a>(), <span class="comment">// solver here</span></div>
<div class="line"> vcl_ilut); <span class="comment">// preconditioner here</span></div>
</div><!-- fragment --><p> The triangular substitutions may be applied in parallel on GPUs by enabling <em>level-scheduling</em> <a class="el" href="citelist.html#CITEREF_saad-iterative-solution">[26]</a> via the member function call <code>use_level_scheduling(true)</code> in the <code>ilu0_config</code> object.</p>
<p>One parameter can be passed to the constructor of <code>ilu0_tag</code>, being the boolean specifying whether level scheduling should be used.</p>
<dl class="section note"><dt>Note</dt><dd>The performance of level scheduling depends strongly on the matrix pattern and is thus disabled by default.</dd></dl>
<h2><a class="anchor" id="manual-algorithms-preconditioners-icc0"></a>
Incomplete Cholesky Factorization with Static Pattern (IChol0)</h2>
<p>Similar to the Chow-Patel-family of incomplete LU and incomplete Cholesky factorization preconditioners, ViennaCL also provides a 'classical' incomplete Cholesky factorization preconditioner with static pattern identical to the system matrix. A sample code snippet is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1linalg_1_1ichol0__tag.html">viennacl::linalg::ichol0_tag</a> ichol0_config;</div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1ichol0__precond.html">viennacl::linalg::ichol0_precond< SparseMatrix ></a> vcl_ilut(A, ichol0_config);</div>
</div><!-- fragment --><p> No level scheduling is currently available for this preconditioner.</p>
<h2><a class="anchor" id="manual-algorithms-preconditioners-block-ilu"></a>
Block-ILU</h2>
<p>To overcome the serial nature of ILUT and ILU0 applied to the full system matrix, a parallel variant is to apply ILU to diagonal blocks of the system matrix. This is accomplished by the <code>block_ilu</code> preconditioner, which takes the system matrix type as first template argument and the respective ILU-tag type as second template argument (either <code>ilut_tag</code> or <code>ilu0_tag</code>). Support for accelerators using CUDA or OpenCL is provided.</p>
<div class="fragment"><div class="line"><span class="comment">// compute block-ILU preconditioner using ILU0 for each block:</span></div>
<div class="line">block_ilu_precond<SparseMatrix, ilu0_tag> vcl_block_ilu0(vcl_matrix, ilu0_tag());</div>
<div class="line"></div>
<div class="line"><span class="comment">// solve</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(vcl_matrix, vcl_rhs,</div>
<div class="line"> <a class="code" href="classviennacl_1_1linalg_1_1bicgstab__tag.html">viennacl::linalg::bicgstab_tag</a>(),</div>
<div class="line"> vcl_block_ilu0);</div>
</div><!-- fragment --><p> A third argument can be passed to the constructor of <code>block_ilu_precond</code>: Either the number of blocks to be used (defaults to <code>8</code>), or an index vector with fine-grained control over the blocks.</p>
<dl class="section note"><dt>Note</dt><dd>The number of blocks is a design parameter for your sparse linear system at hand. Higher number of blocks leads to better memory bandwidth utilization on GPUs, but may increase the number of solver iterations.</dd></dl>
<h2><a class="anchor" id="manual-algorithms-preconditioners-jacobi"></a>
Jacobi Preconditioner</h2>
<p>A Jacobi preconditioner is a simple diagonal preconditioner given by the reciprocals of the diagonal entries of the system matrix. Use the preconditioner as follows: </p>
<div class="fragment"><div class="line"><span class="comment">//compute Jacobi preconditioner:</span></div>
<div class="line">jacobi_precond< SparseMatrix > vcl_jacobi(vcl_matrix, <a class="code" href="classviennacl_1_1linalg_1_1jacobi__tag.html">viennacl::linalg::jacobi_tag</a>());</div>
<div class="line"></div>
<div class="line"><span class="comment">//solve (e.g. using conjugate gradient solver)</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(vcl_matrix, vcl_rhs,</div>
<div class="line"> <a class="code" href="classviennacl_1_1linalg_1_1cg__tag.html">viennacl::linalg::cg_tag</a>(),</div>
<div class="line"> vcl_jacobi);</div>
</div><!-- fragment --><h2><a class="anchor" id="manual-algorithms-preconditioners-row-scaling"></a>
Row-Scaling Preconditioner</h2>
<p>A row scaling preconditioner is a simple diagonal preconditioner given by the reciprocals of the norms of the rows of the system matrix. Use the preconditioner as follows: </p>
<div class="fragment"><div class="line"><span class="comment">//compute row scaling preconditioner:</span></div>
<div class="line">row_scaling< SparseMatrix > vcl_row_scaling(vcl_matrix, <a class="code" href="classviennacl_1_1linalg_1_1row__scaling__tag.html">viennacl::linalg::row_scaling_tag</a>());</div>
<div class="line"></div>
<div class="line"><span class="comment">//solve (e.g. using conjugate gradient solver)</span></div>
<div class="line">vcl_result = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(vcl_matrix, vcl_rhs,</div>
<div class="line"> <a class="code" href="classviennacl_1_1linalg_1_1cg__tag.html">viennacl::linalg::cg_tag</a>(),</div>
<div class="line"> vcl_row_scaling);</div>
</div><!-- fragment --><p> The tag <code><a class="el" href="classviennacl_1_1linalg_1_1row__scaling__tag.html" title="A tag for a row scaling preconditioner which merely normalizes the equation system such that each row...">viennacl::linalg::row_scaling_tag()</a></code> can be supplied with a parameter denoting the norm to be used. A value of <code>1</code> specifies the <img class="formulaInl" alt="$ l^1 $" src="form_27.png"/>-norm, while a value of <img class="formulaInl" alt="$ 2 $" src="form_28.png"/> selects the <img class="formulaInl" alt="$ l^2 $" src="form_29.png"/>-norm (default).</p>
<h2><a class="anchor" id="manual-algorithms-preconditioners-amg"></a>
Algebraic Multigrid Preconditioners</h2>
<p>Algebraic multigrid (AMG) mimics the behavior of geometric multigrid on the algebraic level and is thus suited for black-box purposes, where only the system matrix and the right hand side vector are available <a class="el" href="citelist.html#CITEREF_trottenberg:multigrid">[28]</a> . Many different flavors of the individual multigrid ingredients exist <a class="el" href="citelist.html#CITEREF_yang:parallel-amg">[30]</a>, of which some (with an emphasis on fine-grained parallelism) are implemented in ViennaCL.</p>
<dl class="section note"><dt>Note</dt><dd>Despite the asymptotic optimality of AMG for e.g. systems arising from finite difference or finite element discretizations of elliptic PDEs, AMG methods are not a silver bullet and may just as well fail for general matrices.</dd></dl>
<p>The two main ingredients of algebraic multigrid are a coarsening algorithm and an interpolation algorithm. The available coarsening methods are:</p>
<center> <table class="doxtable">
<tr>
<th>Description </th><th>ViennaCL option in <code><a class="el" href="namespaceviennacl_1_1linalg.html" title="Provides all linear algebra operations which are not covered by operator overloads. ">viennacl::linalg</a>::</code> </th></tr>
<tr>
<td>One-Pass Classical Coarsening </td><td><code>AMG_COARSENING_METHOD_ONEPASS</code> </td></tr>
<tr>
<td>Sequential Aggregation </td><td><code>AMG_COARSENING_METHOD_AGGREGATION</code> </td></tr>
<tr>
<td>Parallel MIS-2 aggregation (default) </td><td><code>AMG_COARSENING_METHOD_MIS2_AGGREGATION</code> </td></tr>
</table>
<p><b>AMG coarsening methods available in ViennaCL. </b> </center><p>The parallel maximum independent set (MIS) algorithm for distance 2 was initially described in <a class="el" href="citelist.html#CITEREF_Bell:AMG">[5]</a> for CUDA. ViennaCL's implementation also supports OpenCL and OpenMP, thus making the method available for all major computing platforms.</p>
<p>The available interpolation methods are: </p>
<center> <table class="doxtable">
<tr>
<th>Description </th><th>ViennaCL option constant </th></tr>
<tr>
<td>Direct interpolation </td><td><code>AMG_INTERPOLATION_METHOD_DIRECT</code> </td></tr>
<tr>
<td>Aggregation </td><td><code>AMG_INTERPOLATION_METHOD_AGGREGATION</code> </td></tr>
<tr>
<td>Smoothed aggregation (default) </td><td><code>AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION</code> </td></tr>
</table>
<p><b>AMG interpolation methods available in ViennaCL.</b> </center><p>The preconditioner setup follows a two-stage process: Parameters for customizing the preconditioner are configured through the <code><a class="el" href="classviennacl_1_1linalg_1_1amg__tag.html" title="A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...">viennacl::linalg::amg_tag</a></code>. Since AMG preconditioners are not a silver bullet, some user customization are typically required for best results. These customizations require a certain familiarity with the concept of multigrid methods. A list of parameters available for tweaks is as follows:</p>
<ul>
<li><b>Strong connection threshold</b>: A relative threshold value above which two nodes in the algebraic graph are considered to be strongly connected.</li>
<li><b>Jacobi smoother weight</b>: Damping parameter for the damped Jacobi method. Parameter values of 0.67 or 1.0 are good starting points for experimentation.</li>
<li><b>Number of pre-smoothing steps</b>: Number of smoother applications on the fine level before restricting the residual to the coarse level.</li>
<li><b>Number of post-smoothing steps</b>: Number of smoother applications after the coarse grid correction has been interpolated back to the fine level.</li>
<li><b>Maximum number of coarse levels</b>: Maximum number of coarse levels to use when setting up the hierarchy. A direct solver is employed on the coarsest level.</li>
<li><b>Coarse level cut-off</b>: Number of unknowns below which the coarsening stops and a direct solver is employed.</li>
<li><b>Context for the preconditioner setup</b>: Explicitly specify the backend to be used for setting up the preconditioner. This way one can e.g. run the setup on the CPU and the preconditioner applications on the GPU.</li>
<li><b>Context for the preconditioner application</b>: Explicitly specify the backend to be used for applying the preconditioner to a vector. This way one can e.g. run the setup on the CPU and the preconditioner applications on the GPU.</li>
</ul>
<p>A typical customization code snippet for running the setup on the CPU and the preconditioner application on a CUDA-enabled GPU is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1context.html">viennacl::context</a> host_ctx(<a class="code" href="namespaceviennacl.html#a00b40450b6b2fd87aad3527d9b2084b8a427356f0fb1b8d32b28f37e36b272df4">viennacl::MAIN_MEMORY</a>);</div>
<div class="line"><a class="code" href="classviennacl_1_1context.html">viennacl::context</a> target_ctx(<a class="code" href="namespaceviennacl.html#a00b40450b6b2fd87aad3527d9b2084b8ab58facda25e2c7e20d9fe1b5e62f46d2">viennacl::CUDA_MEMORY</a>);</div>
<div class="line"></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html">viennacl::linalg::amg_tag</a> amg_tag_direct;</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#a773dab9bf7799da8332b0fe916cc4771">set_coarsening_method</a>(<a class="code" href="namespaceviennacl_1_1linalg.html#a3ba810acdca541a5eada4560982a645ca5167faa7286ca297508f581b26a048dc">viennacl::linalg::AMG_COARSENING_METHOD_ONEPASS</a>);</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#ab33a64319b2dc25a9b2bc38771c7c5a4">set_interpolation_method</a>(<a class="code" href="namespaceviennacl_1_1linalg.html#a9933216144a64dbd433ec02c95bbfdd7a48d5bb166d266affee4d526e8671bec2">viennacl::linalg::AMG_INTERPOLATION_METHOD_DIRECT</a>);</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#ad5b5468b2991f795dabf0d2cc618a72b">set_strong_connection_threshold</a>(0.25);</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#a83d5125ff6aacce5400f7efc3e3fc4d6">set_jacobi_weight</a>(0.67);</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#a39e4042f7fb17c99c33ede3522c9361c">set_presmooth_steps</a>(1);</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#ab00bf522866559edae47e527cb6d6c26">set_postsmooth_steps</a>(1);</div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#a150f1b177671b9daef9ea0d4391c27d8">set_setup_context</a>(host_ctx); <span class="comment">// run setup on host</span></div>
<div class="line">amg_tag_direct.<a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html#a7513d4c67b8878c034dfcd7f4a8f180c">set_target_context</a>(target_ctx); <span class="comment">// run solver cycles on device</span></div>
</div><!-- fragment --><p>The actual preconditioner class is <code><a class="el" href="classviennacl_1_1linalg_1_1amg__precond.html" title="AMG preconditioner class, can be supplied to solve()-routines. ">viennacl::linalg::amg_precond</a><MatrixT></code> for the sparse matrix type <code>MatrixT</code>. Unlike other preconditioner objects, the AMG preconditioner requires an explicit call to the member function <code>.setup()</code> for constructing interpolations, coarse grid hierarchies, etc. A minimal code for using an AMG preconditioner within an iterative solver is as follows: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1compressed__matrix.html">viennacl::compressed_matrix<NumericT></a> A;</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<NumericT></a> b;</div>
<div class="line"></div>
<div class="line"><span class="comment">// Setup A and b here</span></div>
<div class="line"></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1amg__tag.html">viennacl::linalg::amg_tag</a> my_amg_tag;</div>
<div class="line"><span class="comment">// customize my_amg_tag here if needed</span></div>
<div class="line"><a class="code" href="classviennacl_1_1linalg_1_1amg__precond.html">viennacl::linalg::amg_precond<viennacl::compressed_matrix<NumericT></a> > my_amg(A, my_amg_tag);</div>
<div class="line">my_amg.setup();</div>
<div class="line"><a class="code" href="classviennacl_1_1vector.html">viennacl::vector<NumericT></a> x = <a class="code" href="namespaceviennacl_1_1linalg.html#a6e9b329b64ac782e6a5687ad2fc47a2a">viennacl::linalg::solve</a>(A, b, <a class="code" href="classviennacl_1_1linalg_1_1cg__tag.html">viennacl::linalg::cg_tag</a>(), my_amg);</div>
</div><!-- fragment --><dl class="section note"><dt>Note</dt><dd>Note that the efficiency of the various AMG flavors are typically highly problem-specific. Therefore, failure of one method for a particular problem does NOT imply that other coarsening or interpolation strategies will fail as well.</dd></dl>
<h1><a class="anchor" id="manual-algorithms-eigenvalues"></a>
Eigenvalue Computations</h1>
<p>Two algorithms for the computations of the eigenvalues of a sparse matrix are implemented in ViennaCL:</p>
<ul>
<li>The Power Iteration <a class="el" href="citelist.html#CITEREF_golub:matrix-computations">[13]</a></li>
<li>The Lanczos Algorithm <a class="el" href="citelist.html#CITEREF_simon:lanczos-pro">[27]</a></li>
</ul>
<p>The algorithms are called for a matrix object <code>A</code> by </p>
<div class="fragment"><div class="line">std::vector<double> largest_eigenvalues = <a class="code" href="namespaceviennacl_1_1linalg.html#af5002a88fa4cc3fbe65a904797a0ecba">viennacl::linalg::eig</a>(A, ltag);</div>
<div class="line"><span class="keywordtype">double</span> largest_eigenvalue = <a class="code" href="namespaceviennacl_1_1linalg.html#af5002a88fa4cc3fbe65a904797a0ecba">viennacl::linalg::eig</a>(A, ptag);</div>
</div><!-- fragment --><p> Depending on the second parameter <code>tag</code> either of the two methods is called. Both algorithms can be used for either Boost.uBLAS or ViennaCL compressed matrices. In order to get the eigenvalue with the greatest absolut value, the power iteration should be called. The Lanczos algorithm returns a vector of the largest eigenvalues with the same type as the entries of the matrix.</p>
<h2><a class="anchor" id="manual-algorithms-eigenvalues-power"></a>
Power Iteration</h2>
<p>The Power iteration aims at computing the eigenvalues of a matrix by calculating the product of the matrix and a vector for several times, where the resulting vector is used for the next product of the matrix and so on. The computation stops as soon as the norm of the vector changes by less than the prescribed tolerance. The final vector is the eigenvector to the eigenvalue with the greatest absolut value. To call this algorithm, <code>piter_tag</code> must be used. This tag has only one parameter: <code>terminationfactor</code> defines the accuracy of the computation, i.e. if the new norm of the eigenvector changes less than this parameter the computation stops and returns the corresponding eigenvalue (default: <img class="formulaInl" alt="$ 10^{-10} $" src="form_23.png"/> ) The call of the constructor may look as follows: </p>
<div class="fragment"><div class="line">viennacl::linalg::piter_tag ptag(1e-8);</div>
</div><!-- fragment --><p> A minimal example is as follows: </p>
<div class="fragment"><div class="line"><span class="preprocessor">#include "<a class="code" href="compressed__matrix_8hpp.html">viennacl/compressed_matrix.hpp</a>"</span></div>
<div class="line"><span class="preprocessor">#include "<a class="code" href="vector_8hpp.html">viennacl/vector.hpp</a>"</span></div>
<div class="line"><span class="preprocessor">#include "<a class="code" href="power__iter_8hpp.html">viennacl/linalg/power_iter.hpp</a>"</span></div>
<div class="line"></div>
<div class="line"><span class="keywordtype">int</span> <a class="code" href="tests_2src_2bisect_8cpp.html#ae66f6b31b5ad750f1fe042a706a4e3d4">main</a>()</div>
<div class="line">{</div>
<div class="line"> <a class="code" href="classviennacl_1_1compressed__matrix.html">viennacl::compressed_matrix<double></a> A(N, N);</div>
<div class="line"> <a class="code" href="classviennacl_1_1vector.html">viennacl::vector<double></a> eigenvector(N);</div>
<div class="line"></div>
<div class="line"> <span class="comment">// fill A with data here, e.g. by reading a MatrixMarket file</span></div>
<div class="line"></div>
<div class="line"> viennacl::linalg::piter_tag ptag(1e-8);</div>
<div class="line"> std::cout << <span class="stringliteral">"Eigenvalue: "</span> << <a class="code" href="namespaceviennacl_1_1linalg.html#af5002a88fa4cc3fbe65a904797a0ecba">viennacl::linalg::eig</a>(A, ptag, eigenvector) << std::endl;</div>
<div class="line"> std::cout << <span class="stringliteral">"Eigenvector: "</span> << eigenvector << std::endl;</div>
<div class="line">}</div>
</div><!-- fragment --><p> Note that the eigenvalue and eigenvector returned are only approximate. Their accuracy depends strongly on the prescribed tolerance and separation of the largest eigenvalue from the second-largest eigenvalue (both in modulus).</p>
<dl class="section note"><dt>Note</dt><dd>Example code can be found in <code><a class="el" href="power-iter_8cpp.html">examples/tutorial/power-iter.cpp</a></code></dd></dl>
<h2><a class="anchor" id="manual-algorithms-eigenvalues-lanczos"></a>
The Lanczos Algorithm</h2>
<p>In order to compute the eigenvalues of a sparse high-dimensional matrix the Lanczos algorithm can be used to find these. This algorithm reformulates the given high-dimensional matrix in a way such that the matrix can be rewritten in a tridiagonal matrix at much lower dimension. The eigenvalues of this tridiagonal matrix are equal to the largest eigenvalues of the original matrix and calculated by using the bisection method <a class="el" href="citelist.html#CITEREF_golub:matrix-computations">[13]</a> . To call this Lanczos algorithm, <code>lanczos_tag</code> must be used. This tag has several parameters that can be passed to the constructor:</p>
<ul>
<li>The exponent of epsilon for the tolerance of the reorthogonalization, defined by the parameter <code>factor</code> (default: <code>0.75</code>)</li>
<li>The method of the Lanczos algorithm: <code>0</code> uses partial reorthogonalization, <code>1</code> full reothogonalization, and <code>2</code> does not use reorthogonalization (default: <code>0</code>)</li>
<li>The number of eigenvalues that are returned is specified by <code>num_eigenvalues</code> (default: <code>10</code>)</li>
<li>The size of the Krylov space used for the computations can be set by the parameter <code>krylov_size</code> (default: <code>100</code>). The maximum number of iterations can be equal or less this parameter.</li>
</ul>
<p>The call of the constructor may look like the following: </p>
<div class="fragment"><div class="line"><a class="code" href="classviennacl_1_1linalg_1_1lanczos__tag.html">viennacl::linalg::lanczos_tag</a> ltag(0.85, 15, 0, 200);</div>
</div><!-- fragment --><dl class="section note"><dt>Note</dt><dd>Example code can be found in <code><a class="el" href="lanczos_8cpp.html">examples/tutorial/lanczos.cpp</a></code></dd></dl>
<h1><a class="anchor" id="manual-algorithms-qr-factorization"></a>
QR Factorization</h1>
<dl class="section note"><dt>Note</dt><dd>The current QR factorization implementation depends on Boost.uBLAS.</dd></dl>
<p>A matrix <img class="formulaInl" alt="$ A \in \mathbb{R}^{n\times m} $" src="form_30.png"/> can be factored into <img class="formulaInl" alt="$ A = Q R $" src="form_31.png"/>, where <img class="formulaInl" alt="$ Q \in \mathbb{R}^{n\times n}$" src="form_32.png"/> is an orthogonal matrix and <img class="formulaInl" alt="$ R \in \mathbb{R}^{n \times m}$" src="form_33.png"/> is upper triangular. This so-called QR-factorization is important for eigenvalue computations as well as for the solution of least-squares problems <a class="el" href="citelist.html#CITEREF_golub:matrix-computations">[13]</a> . ViennaCL provides a generic implementation of the QR-factorization using Householder reflections in file <code><a class="el" href="qr_8hpp.html" title="Provides a QR factorization using a block-based approach. ">viennacl/linalg/qr.hpp</a></code>. An example application can be found in <code><a class="el" href="tutorial_2qr_8cpp.html">examples/tutorial/qr.cpp</a></code>.</p>
<p>The Householder reflectors <img class="formulaInl" alt="$ v_i $" src="form_34.png"/> defining the Householder reflection <img class="formulaInl" alt="$ I - \beta_i v_i v_i^{\mathrm{T}} $" src="form_35.png"/> are stored in the columns below the diagonal of the input matrix <img class="formulaInl" alt="$ A $" src="form_1.png"/> <a class="el" href="citelist.html#CITEREF_golub:matrix-computations">[13]</a> . The normalization coefficients <img class="formulaInl" alt="$ \beta_i $" src="form_36.png"/> are returned by the worker function <code>inplace_qr</code>. The upper triangular matrix <img class="formulaInl" alt="$ R $" src="form_37.png"/> is directly written to the upper triangular part of <img class="formulaInl" alt="$ A $" src="form_1.png"/>. </p>
<div class="fragment"><div class="line">std::vector<ScalarType> betas = <a class="code" href="namespaceviennacl_1_1linalg.html#a83f573386903bae3a5379ab4cbfc5e2d">viennacl::linalg::inplace_qr</a>(A, 12);</div>
</div><!-- fragment --><p> If <img class="formulaInl" alt="$ A $" src="form_1.png"/> is a dense matrix from Boost.uBLAS, the calculation is carried out on the CPU using a single thread. If <img class="formulaInl" alt="$ A $" src="form_1.png"/> is a <code><a class="el" href="classviennacl_1_1matrix.html" title="A dense matrix class. ">viennacl::matrix</a></code>, a hybrid implementation is used: The panel factorization is carried out using Boost.uBLAS, while expensive BLAS level 3 operations are computed on the OpenCL device using multiple threads.</p>
<p>Typically, the orthogonal matrix <img class="formulaInl" alt="$ Q $" src="form_38.png"/> is kept in inplicit form because of computational efficiency. However, if <img class="formulaInl" alt="$ Q $" src="form_38.png"/> and <img class="formulaInl" alt="$ R $" src="form_37.png"/> have to be computed explicitly, the function <code>recoverQ</code> can be used: </p>
<div class="fragment"><div class="line"><a class="code" href="namespaceviennacl_1_1linalg.html#a76fc8bf8f5a20c576b7ddbad32d3e56d">viennacl::linalg::recoverQ</a>(A, betas, Q, R);</div>
</div><!-- fragment --><p> Here, <code>A</code> is the inplace QR-factored matrix, <code>betas</code> are the coefficients of the Householder reflectors as returned by <code>inplace_qr</code>, while <code>Q</code> and <code>R</code> are the destination matrices. However, the explicit formation of <code>Q</code> is expensive and is usually avoided. For a number of applications of the QR factorization it is only required to apply <code>Q^T</code> to a vector <code>b</code>. This is accomplished by </p>
<div class="fragment"><div class="line"><a class="code" href="namespaceviennacl_1_1linalg.html#adf77bb0fe92ef09522d0458ee9df03db">viennacl::linalg::inplace_qr_apply_trans_Q</a>(A, betas, b);</div>
</div><!-- fragment --><p> without setting up <code>Q</code> (or <code>Q^T</code>) explicitly.</p>
<dl class="section note"><dt>Note</dt><dd>Have a look at <code><a class="el" href="least-squares_8cpp.html">examples/tutorial/least-squares.cpp</a></code> for a least-squares computation using QR factorizations. </dd></dl>
</div></div><!-- contents -->
</div><!-- doc-content -->
<!-- start footer part -->
<div id="nav-path" class="navpath"><!-- id is needed for treeview function! -->
<ul>
<li class="footer">Generated on Wed Jan 20 2016 22:32:44 for ViennaCL - The Vienna Computing Library by
<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.8.6 </li>
</ul>
</div>
</body>
</html>
|
