ViennaCL - The Vienna Computing Library  1.7.1
Free open-source GPU-accelerated linear algebra and solver library.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gerschgorin.hpp
Go to the documentation of this file.
1 #ifndef _VIENNACL_LINALG_DETAIL_BISECT_GERSCHORIN_HPP_
2 #define _VIENNACL_LINALG_DETAIL_BISECT_GERSCHORIN_HPP_
3 
4 
5 /* =========================================================================
6  Copyright (c) 2010-2016, Institute for Microelectronics,
7  Institute for Analysis and Scientific Computing,
8  TU Wien.
9  Portions of this software are copyright by UChicago Argonne, LLC.
10 
11  -----------------
12  ViennaCL - The Vienna Computing Library
13  -----------------
14 
15  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
16 
17  (A list of authors and contributors can be found in the manual)
18 
19  License: MIT (X11), see file LICENSE in the base directory
20 ============================================================================= */
21 
22 
31 #include <cstdio>
32 #include <cstdlib>
33 #include <cmath>
34 #include <cfloat>
35 
37 #include "viennacl/vector.hpp"
38 
39 namespace viennacl
40 {
41 namespace linalg
42 {
43 namespace detail
44 {
53  template<typename NumericT>
54  void
55  computeGerschgorin(std::vector<NumericT> & d, std::vector<NumericT> & s, unsigned int n, NumericT &lg, NumericT &ug)
56  {
57  // compute bounds
58  for (unsigned int i = 1; i < (n - 1); ++i)
59  {
60 
61  // sum over the absolute values of all elements of row i
62  NumericT sum_abs_ni = fabsf(s[i]) + fabsf(s[i + 1]);
63 
64  lg = min(lg, d[i] - sum_abs_ni);
65  ug = max(ug, d[i] + sum_abs_ni);
66  }
67 
68  // first and last row, only one superdiagonal element
69 
70  // first row
71  lg = min(lg, d[0] - fabsf(s[1]));
72  ug = max(ug, d[0] + fabsf(s[1]));
73 
74  // last row
75  lg = min(lg, d[n-1] - fabsf(s[n-1]));
76  ug = max(ug, d[n-1] + fabsf(s[n-1]));
77 
78  // increase interval to avoid side effects of fp arithmetic
79  NumericT bnorm = max(fabsf(ug), fabsf(lg));
80 
81  // these values depend on the implmentation of floating count that is
82  // employed in the following
83  NumericT psi_0 = 11 * FLT_EPSILON * bnorm;
84  NumericT psi_n = 11 * FLT_EPSILON * bnorm;
85 
86  lg = lg - bnorm * 2 * static_cast<NumericT>(n) * FLT_EPSILON - psi_0;
87  ug = ug + bnorm * 2 * static_cast<NumericT>(n) * FLT_EPSILON + psi_n;
88 
89  ug = max(lg, ug);
90  }
91 } // namespace detail
92 } // namespace linalg
93 } // namespace viennacl
94 #endif // _VIENNACL_LINALG_DETAIL_GERSCHORIN_H_
Utility functions.
void computeGerschgorin(std::vector< NumericT > &d, std::vector< NumericT > &s, unsigned int n, NumericT &lg, NumericT &ug)
Definition: gerschgorin.hpp:55
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
float NumericT
Definition: bisect.cpp:40
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
T min(const T &lhs, const T &rhs)
Minimum.
Definition: util.hpp:45