fml  0.1-0
Fused Matrix Library
norm.hh
1 // This file is part of fml which is released under the Boost Software
2 // License, Version 1.0. See accompanying file LICENSE or copy at
3 // https://www.boost.org/LICENSE_1_0.txt
4 
5 #ifndef FML_CPU_LINALG_LINALG_NORM_H
6 #define FML_CPU_LINALG_LINALG_NORM_H
7 #pragma once
8 
9 
10 #include <cmath>
11 #include <stdexcept>
12 
13 #include "../../_internals/linalgutils.hh"
14 #include "../../_internals/omp.hh"
15 
16 #include "../cpumat.hh"
17 #include "../cpuvec.hh"
18 
19 #include "internals/lapack.hh"
20 #include "qr.hh"
21 #include "svd.hh"
22 
23 
24 namespace fml
25 {
26 namespace linalg
27 {
37  template <typename REAL>
38  REAL norm_1(const cpumat<REAL> &x)
39  {
40  const len_t m = x.nrows();
41  const len_t n = x.ncols();
42  const REAL *x_d = x.data_ptr();
43 
44  REAL norm = 0;
45 
46  #pragma omp parallel for if(m*n > fml::omp::OMP_MIN_SIZE) reduction(max:norm)
47  for (len_t j=0; j<n; j++)
48  {
49  REAL tmp = 0;
50 
51  #pragma omp simd reduction(+:tmp)
52  for (len_t i=0; i<m; i++)
53  tmp += fabs(x_d[i + m*j]);
54 
55  if (tmp > norm)
56  norm = tmp;
57  }
58 
59  return norm;
60  }
61 
62 
63 
78  template <typename REAL>
79  REAL norm_I(const cpumat<REAL> &x)
80  {
81  const len_t m = x.nrows();
82  const len_t n = x.ncols();
83  const REAL *x_d = x.data_ptr();
84 
85  REAL norm = 0;
86 
87  cpuvec<REAL> mars(m);
88  mars.fill_zero();
89  REAL *mars_d = mars.data_ptr();
90 
91  for (len_t j=0; j<n; j++)
92  {
93  for (len_t i=0; i<m; i++)
94  mars_d[i] += fabs(x_d[i + m*j]);
95  }
96 
97  for (len_t i=0; i<m; i++)
98  {
99  if (mars_d[i] > norm)
100  norm = mars_d[i];
101  }
102 
103  return norm;
104  }
105 
106 
107 
117  template <typename REAL>
118  REAL norm_F(const cpumat<REAL> &x)
119  {
120  const len_t m = x.nrows();
121  const len_t n = x.ncols();
122  const REAL *x_d = x.data_ptr();
123 
124  REAL scale = 0;
125  REAL sumsq = 1;
126 
127  for (len_t j=0; j<n; j++)
128  lapack::lassq(m, x_d + m*j, 1, &scale, &sumsq);
129 
130  return scale * sqrtf(sumsq);
131  }
132 
133 
134 
144  template <typename REAL>
145  REAL norm_M(const cpumat<REAL> &x)
146  {
147  const len_t m = x.nrows();
148  const len_t n = x.ncols();
149  const REAL *x_d = x.data_ptr();
150 
151  REAL norm = 0;
152 
153  #pragma omp parallel for if(m*n > fml::omp::OMP_MIN_SIZE) reduction(max:norm)
154  for (len_t j=0; j<n; j++)
155  {
156  for (len_t i=0; i<m; i++)
157  {
158  REAL tmp = fabs(x_d[i + m*j]);
159  if (tmp > norm)
160  norm = tmp;
161  }
162  }
163 
164  return norm;
165  }
166 
167 
168 
187  template <typename REAL>
189  {
190  REAL ret;
191  cpuvec<REAL> s;
192  svd(x, s);
193  ret = s.get(0);
194 
195  return ret;
196  }
197 
198 
199 
200  namespace
201  {
202  template <typename REAL>
203  REAL cond_square_internals(const char norm, cpumat<REAL> &x)
204  {
205  const len_t n = x.nrows();
206 
207  REAL ret;
208  int info;
209 
210  lu(x);
211 
212  cpuvec<REAL> work(4*n);
213  cpuvec<REAL> work2(n);
214  REAL xnorm = norm_1(x);
215  lapack::gecon(norm, n, x.data_ptr(), n, xnorm, &ret, work.data_ptr(),
216  work2.data_ptr(), &info);
217 
218  fml::linalgutils::check_info(info, "gecon");
219 
220  return ((REAL)1)/ret;
221  }
222 
223  template <typename REAL>
224  REAL cond_nonsquare_internals(const char norm, cpumat<REAL> &x)
225  {
226  const len_t m = x.nrows();
227  const len_t n = x.ncols();
228 
229  REAL ret;
230  int info;
231 
232  cpuvec<REAL> aux;
233 
234  if (m > n)
235  {
236  cpumat<REAL> R;
237  qr(false, x, aux);
238  qr_R(x, R);
239 
240  aux.resize(R.nrows());
241  lapack::trcon(norm, 'U', 'N', n, R.data_ptr(), n, &ret,
242  x.data_ptr(), aux.data_ptr(), &info);
243  }
244  else
245  {
246  cpumat<REAL> L;
247  lq(x, aux);
248  lq_L(x, L);
249 
250  aux.resize(L.nrows());
251  lapack::trcon(norm, 'L', 'N', m, L.data_ptr(), m, &ret,
252  x.data_ptr(), aux.data_ptr(), &info);
253  }
254 
255  fml::linalgutils::check_info(info, "trcon");
256 
257  return ((REAL)1)/ret;
258  }
259  }
260 
280  template <typename REAL>
282  {
283  if (x.is_square())
284  return cond_square_internals('1', x);
285  else
286  return cond_nonsquare_internals('1', x);
287  }
288 
289 
290 
310  template <typename REAL>
312  {
313  if (x.is_square())
314  return cond_square_internals('I', x);
315  else
316  return cond_nonsquare_internals('I', x);
317  }
318 
319 
320 
338  template <typename REAL>
340  {
341  cpuvec<REAL> s;
342  svd(x, s);
343 
344  REAL *s_d = s.data_ptr();
345 
346  REAL max = s_d[0];
347  REAL min = s_d[0];
348 
349  #pragma omp parallel for if(s.size() > fml::omp::OMP_MIN_SIZE) reduction(max:max) reduction(min:min)
350  for (len_t i=1; i<s.size(); i++)
351  {
352  if (s_d[i] > max)
353  max = s_d[i];
354  if (s_d[i] > 0 && s_d[i] < min)
355  min = s_d[i];
356  }
357 
358  if (max == 0)
359  return 0;
360  else
361  return max/min;
362  }
363 }
364 }
365 
366 
367 #endif
fml::cpumat
Matrix class for data held on a single CPU.
Definition: cpumat.hh:36
fml::linalg::norm_I
REAL norm_I(const cpumat< REAL > &x)
Computes the infinity matrix norm, the maximum absolute row sum.
Definition: norm.hh:79
fml::linalg::cond_2
REAL cond_2(cpumat< REAL > &x)
Estimates the condition number under the 2 norm.
Definition: norm.hh:339
fml::unimat::is_square
bool is_square() const
Is the matrix square?
Definition: unimat.hh:34
fml::univec::data_ptr
T * data_ptr()
Pointer to the internal array.
Definition: univec.hh:28
fml::unimat::nrows
len_t nrows() const
Number of rows.
Definition: unimat.hh:36
fml::linalg::cond_1
REAL cond_1(cpumat< REAL > &x)
Estimates the condition number under the 1-norm.
Definition: norm.hh:281
fml::linalg::qr
void qr(const bool pivot, cpumat< REAL > &x, cpuvec< REAL > &qraux)
Computes the QR decomposition.
Definition: qr.hh:94
fml::linalg::norm_1
REAL norm_1(const cpumat< REAL > &x)
Computes the 1 matrix norm, the maximum absolute column sum.
Definition: norm.hh:38
fml::linalg::lu
void lu(cpumat< REAL > &x, cpuvec< int > &p, int &info)
Computes the PLU factorization with partial pivoting.
Definition: lu.hh:48
fml::cpuvec
Vector class for data held on a single CPU.
Definition: cpuvec.hh:31
fml::unimat::ncols
len_t ncols() const
Number of columns.
Definition: unimat.hh:38
fml::cpuvec::get
T get(const len_t i) const
Get the specified value.
Definition: cpuvec.hh:524
fml::linalg::cond_I
REAL cond_I(cpumat< REAL > &x)
Estimates the condition number under the infinity norm.
Definition: norm.hh:311
fml::linalg::svd
void svd(cpumat< REAL > &x, cpuvec< REAL > &s)
Computes the singular value decomposition.
Definition: svd.hh:101
fml::linalg::norm_2
REAL norm_2(cpumat< REAL > &x)
Computes the 2/spectral matrix norm.
Definition: norm.hh:188
fml::unimat::data_ptr
REAL * data_ptr()
Pointer to the internal array.
Definition: unimat.hh:40
fml::linalg::lq_L
void lq_L(const cpumat< REAL > &LQ, cpumat< REAL > &L)
Recover the L matrix from an LQ decomposition.
Definition: qr.hh:247
fml
Core namespace.
Definition: dimops.hh:10
fml::univec::size
len_t size() const
Number of elements in the vector.
Definition: univec.hh:26
fml::linalg::qr_R
void qr_R(const cpumat< REAL > &QR, cpumat< REAL > &R)
Recover the R matrix from a QR decomposition.
Definition: qr.hh:162
fml::linalg::lq
void lq(cpumat< REAL > &x, cpuvec< REAL > &lqaux)
Computes the LQ decomposition.
Definition: qr.hh:223
fml::cpuvec::fill_zero
void fill_zero()
Set all values to zero.
Definition: cpuvec.hh:317
fml::linalg::norm_M
REAL norm_M(const cpumat< REAL > &x)
Computes the maximum modulus matrix norm.
Definition: norm.hh:145
fml::linalg::norm_F
REAL norm_F(const cpumat< REAL > &x)
Computes the Frobenius/Euclidean matrix norm.
Definition: norm.hh:118