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_MPI_LINALG_NORM_H
6 #define FML_MPI_LINALG_NORM_H
7 #pragma once
8 
9 
10 #include "../../cpu/cpuvec.hh"
11 
12 #include "../mpimat.hh"
13 
14 #include "lu.hh"
15 #include "qr.hh"
16 #include "svd.hh"
17 
18 
19 namespace fml
20 {
21 namespace linalg
22 {
37  template <typename REAL>
38  REAL norm_1(const mpimat<REAL> &x)
39  {
40  const len_t n = x.nrows();
41  const len_t m_local = x.nrows_local();
42  const len_t n_local = x.ncols_local();
43  const REAL *x_d = x.data_ptr();
44  const grid g = x.get_grid();
45 
46  REAL norm = 0;
47 
48  cpuvec<REAL> macs(n);
49  macs.fill_zero();
50  REAL *macs_d = macs.data_ptr();
51 
52  #pragma omp parallel for if(m_local*n_local > fml::omp::OMP_MIN_SIZE)
53  for (len_t j=0; j<n_local; j++)
54  {
55  for (len_t i=0; i<m_local; i++)
56  macs_d[j] += fabs(x_d[i + m_local*j]);
57  }
58 
59  g.allreduce(n, 1, macs_d, 'C');
60 
61  for (len_t j=0; j<n; j++)
62  {
63  if (norm < macs_d[j])
64  norm = macs_d[j];
65  }
66 
67  g.allreduce(1, 1, &norm, 'R', BLACS_MAX);
68 
69  return norm;
70  }
71 
72 
73 
88  template <typename REAL>
89  REAL norm_I(const mpimat<REAL> &x)
90  {
91  const len_t m = x.nrows();
92  const len_t m_local = x.nrows_local();
93  const len_t n_local = x.ncols_local();
94  const REAL *x_d = x.data_ptr();
95  const grid g = x.get_grid();
96 
97  REAL norm = 0;
98 
99  cpuvec<REAL> mars(m);
100  mars.fill_zero();
101  REAL *mars_d = mars.data_ptr();
102 
103  for (len_t j=0; j<n_local; j++)
104  {
105  for (len_t i=0; i<m_local; i++)
106  mars_d[i] += fabs(x_d[i + m_local*j]);
107  }
108 
109  g.allreduce(m, 1, mars_d, 'R');
110 
111  for (len_t i=0; i<m; i++)
112  {
113  if (norm < mars_d[i])
114  norm = mars_d[i];
115  }
116 
117  g.allreduce(1, 1, &norm, 'C', BLACS_MAX);
118 
119  return norm;
120  }
121 
122 
123 
133  template <typename REAL>
134  REAL norm_F(const mpimat<REAL> &x)
135  {
136  const len_t m_local = x.nrows_local();
137  const len_t n_local = x.ncols_local();
138  const REAL *x_d = x.data_ptr();
139 
140  REAL norm = 0;
141 
142  for (len_t j=0; j<n_local; j++)
143  {
144  for (len_t i=0; i<m_local; i++)
145  norm += x_d[i + m_local*j] * x_d[i + m_local*j];
146  }
147 
148  x.get_grid().allreduce(1, 1, &norm, 'A');
149 
150  return sqrt(norm);
151  }
152 
153 
154 
164  template <typename REAL>
165  REAL norm_M(const mpimat<REAL> &x)
166  {
167  const len_t m_local = x.nrows_local();
168  const len_t n_local = x.ncols_local();
169  const REAL *x_d = x.data_ptr();
170 
171  REAL norm = 0;
172 
173  for (len_t j=0; j<n_local; j++)
174  {
175  for (len_t i=0; i<m_local; i++)
176  {
177  REAL tmp = fabs(x_d[i + m_local*j]);
178  if (tmp > norm)
179  norm = tmp;
180  }
181  }
182 
183  x.get_grid().allreduce(1, 1, &norm, 'A', BLACS_MAX);
184 
185  return norm;
186  }
187 
188 
189 
206  template <typename REAL>
208  {
209  REAL ret;
210  cpuvec<REAL> s;
211  svd(x, s);
212  ret = s.get(0);
213 
214  return ret;
215  }
216 
217 
218 
219  namespace
220  {
221  template <typename REAL>
222  REAL cond_square_internals(const char norm, mpimat<REAL> &x)
223  {
224  const len_t n = x.nrows();
225 
226  REAL ret;
227  int info;
228 
229  REAL xnorm = norm_1(x);
230 
231  lu(x);
232 
233  REAL tmp;
234  int liwork;
235  scalapack::gecon(norm, n, x.data_ptr(), x.desc_ptr(), xnorm,
236  &ret, &tmp, -1, &liwork, -1, &info);
237 
238  int lwork = (int) tmp;
239  cpuvec<REAL> work(lwork);
240  cpuvec<int> iwork(liwork);
241 
242  scalapack::gecon(norm, n, x.data_ptr(), x.desc_ptr(), xnorm,
243  &ret, work.data_ptr(), lwork, iwork.data_ptr(), liwork, &info);
244 
245  fml::linalgutils::check_info(info, "gecon");
246 
247  return ((REAL)1)/ret;
248  }
249 
250  template <typename REAL>
251  REAL cond_nonsquare_internals(const char norm, mpimat<REAL> &x)
252  {
253  const len_t m = x.nrows();
254  const len_t n = x.ncols();
255 
256  REAL ret;
257  int info;
258 
259  cpuvec<REAL> aux;
260 
261  if (m > n)
262  {
263  mpimat<REAL> R(x.get_grid(), x.bf_rows(), x.bf_cols());
264  qr(false, x, aux);
265  qr_R(x, R);
266 
267  REAL tmp;
268  int liwork;
269  scalapack::trcon(norm, 'U', 'N', n, R.data_ptr(), R.desc_ptr(), &ret,
270  &tmp, -1, &liwork, -1, &info);
271 
272  int lwork = (int) tmp;
273  cpuvec<REAL> work(lwork);
274  cpuvec<int> iwork(liwork);
275 
276  scalapack::trcon(norm, 'U', 'N', n, R.data_ptr(), R.desc_ptr(), &ret,
277  x.data_ptr(), lwork, iwork.data_ptr(), iwork.size(), &info);
278  }
279  else
280  {
281  mpimat<REAL> L(x.get_grid(), x.bf_rows(), x.bf_cols());
282  lq(x, aux);
283  lq_L(x, L);
284 
285  REAL tmp;
286  int liwork;
287  scalapack::trcon(norm, 'L', 'N', m, L.data_ptr(), L.desc_ptr(), &ret,
288  &tmp, -1, &liwork, -1, &info);
289 
290  int lwork = (int) tmp;
291  cpuvec<REAL> work(lwork);
292  cpuvec<int> iwork(liwork);
293 
294  scalapack::trcon(norm, 'L', 'N', m, L.data_ptr(), L.desc_ptr(), &ret,
295  x.data_ptr(), x.nrows()*x.ncols(), iwork.data_ptr(), iwork.size(), &info);
296  }
297 
298  fml::linalgutils::check_info(info, "trcon");
299 
300  return ((REAL)1)/ret;
301  }
302  }
303 
323  template <typename REAL>
325  {
326  if (x.is_square())
327  return cond_square_internals('1', x);
328  else
329  return cond_nonsquare_internals('1', x);
330  }
331 
332 
333 
353  template <typename REAL>
355  {
356  if (x.is_square())
357  return cond_square_internals('I', x);
358  else
359  return cond_nonsquare_internals('I', x);
360  }
361 
362 
363 
381  template <typename REAL>
383  {
384  cpuvec<REAL> s;
385  svd(x, s);
386 
387  REAL *s_d = s.data_ptr();
388 
389  REAL max = s_d[0];
390  REAL min = s_d[0];
391  for (len_t i=1; i<s.size(); i++)
392  {
393  if (s_d[i] > max)
394  max = s_d[i];
395  if (s_d[i] > 0 && s_d[i] < min)
396  min = s_d[i];
397  }
398 
399  if (max == 0)
400  return 0;
401  else
402  return max/min;
403  }
404 }
405 }
406 
407 
408 #endif
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::grid::allreduce
void allreduce(const int m, const int n, int *x, const char scope='A', const blacsops op=BLACS_SUM) const
Sum reduce operation across all processes in the grid.
Definition: grid.hh:420
fml::linalg::cond_2
REAL cond_2(cpumat< REAL > &x)
Estimates the condition number under the 2 norm.
Definition: norm.hh:339
fml::grid
2-dimensional MPI process grid.
Definition: grid.hh:70
fml::mpimat
Matrix class for data distributed over MPI in the 2-d block cyclic format.
Definition: mpimat.hh:40
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::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