fml  0.1-0
Fused Matrix Library
linalg_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_LINALG_NORM_H
6 #define FML_MPI_LINALG_LINALG_NORM_H
7 #pragma once
8 
9 
10 #include "../../cpu/cpuvec.hh"
11 
12 #include "../mpimat.hh"
13 
14 
15 namespace fml
16 {
17 namespace linalg
18 {
33  template <typename REAL>
34  REAL norm_1(const mpimat<REAL> &x)
35  {
36  const len_t n = x.nrows();
37  const len_t m_local = x.nrows_local();
38  const len_t n_local = x.ncols_local();
39  const REAL *x_d = x.data_ptr();
40  const grid g = x.get_grid();
41 
42  REAL norm = 0;
43 
44  cpuvec<REAL> macs(n);
45  macs.fill_zero();
46  REAL *macs_d = macs.data_ptr();
47 
48  #pragma omp parallel for if(m_local*n_local > fml::omp::OMP_MIN_SIZE)
49  for (len_t j=0; j<n_local; j++)
50  {
51  for (len_t i=0; i<m_local; i++)
52  macs_d[j] += fabs(x_d[i + m_local*j]);
53  }
54 
55  g.allreduce(n, 1, macs_d, 'C');
56 
57  for (len_t j=0; j<n; j++)
58  {
59  if (norm < macs_d[j])
60  norm = macs_d[j];
61  }
62 
63  g.allreduce(1, 1, &norm, 'R', BLACS_MAX);
64 
65  return norm;
66  }
67 
68 
69 
84  template <typename REAL>
85  REAL norm_I(const mpimat<REAL> &x)
86  {
87  const len_t m = x.nrows();
88  const len_t m_local = x.nrows_local();
89  const len_t n_local = x.ncols_local();
90  const REAL *x_d = x.data_ptr();
91  const grid g = x.get_grid();
92 
93  REAL norm = 0;
94 
95  cpuvec<REAL> mars(m);
96  mars.fill_zero();
97  REAL *mars_d = mars.data_ptr();
98 
99  for (len_t j=0; j<n_local; j++)
100  {
101  for (len_t i=0; i<m_local; i++)
102  mars_d[i] += fabs(x_d[i + m_local*j]);
103  }
104 
105  g.allreduce(m, 1, mars_d, 'R');
106 
107  for (len_t i=0; i<m; i++)
108  {
109  if (norm < mars_d[i])
110  norm = mars_d[i];
111  }
112 
113  g.allreduce(1, 1, &norm, 'C', BLACS_MAX);
114 
115  return norm;
116  }
117 
118 
119 
129  template <typename REAL>
130  REAL norm_F(const mpimat<REAL> &x)
131  {
132  const len_t m_local = x.nrows_local();
133  const len_t n_local = x.ncols_local();
134  const REAL *x_d = x.data_ptr();
135 
136  REAL norm = 0;
137 
138  for (len_t j=0; j<n_local; j++)
139  {
140  for (len_t i=0; i<m_local; i++)
141  norm += x_d[i + m_local*j] * x_d[i + m_local*j];
142  }
143 
144  x.get_grid().allreduce(1, 1, &norm, 'A');
145 
146  return sqrt(norm);
147  }
148 
149 
150 
160  template <typename REAL>
161  REAL norm_M(const mpimat<REAL> &x)
162  {
163  const len_t m_local = x.nrows_local();
164  const len_t n_local = x.ncols_local();
165  const REAL *x_d = x.data_ptr();
166 
167  REAL norm = 0;
168 
169  for (len_t j=0; j<n_local; j++)
170  {
171  for (len_t i=0; i<m_local; i++)
172  {
173  REAL tmp = fabs(x_d[i + m_local*j]);
174  if (tmp > norm)
175  norm = tmp;
176  }
177  }
178 
179  x.get_grid().allreduce(1, 1, &norm, 'A', BLACS_MAX);
180 
181  return norm;
182  }
183 
184 
185 
202  template <typename REAL>
204  {
205  REAL ret;
206  cpuvec<REAL> s;
207  svd(x, s);
208  ret = s.get(0);
209 
210  return ret;
211  }
212 
213 
214 
215  namespace
216  {
217  template <typename REAL>
218  REAL cond_square_internals(const char norm, mpimat<REAL> &x)
219  {
220  const len_t n = x.nrows();
221 
222  REAL ret;
223  int info;
224 
225  REAL xnorm = norm_1(x);
226 
227  lu(x);
228 
229  REAL tmp;
230  int liwork;
231  scalapack::gecon(norm, n, x.data_ptr(), x.desc_ptr(), xnorm,
232  &ret, &tmp, -1, &liwork, -1, &info);
233 
234  int lwork = (int) tmp;
235  cpuvec<REAL> work(lwork);
236  cpuvec<int> iwork(liwork);
237 
238  scalapack::gecon(norm, n, x.data_ptr(), x.desc_ptr(), xnorm,
239  &ret, work.data_ptr(), lwork, iwork.data_ptr(), liwork, &info);
240 
241  fml::linalgutils::check_info(info, "gecon");
242 
243  return ((REAL)1)/ret;
244  }
245 
246  template <typename REAL>
247  REAL cond_nonsquare_internals(const char norm, mpimat<REAL> &x)
248  {
249  const len_t m = x.nrows();
250  const len_t n = x.ncols();
251 
252  REAL ret;
253  int info;
254 
255  cpuvec<REAL> aux;
256 
257  if (m > n)
258  {
259  mpimat<REAL> R(x.get_grid(), x.bf_rows(), x.bf_cols());
260  qr(false, x, aux);
261  qr_R(x, R);
262 
263  REAL tmp;
264  int liwork;
265  scalapack::trcon(norm, 'U', 'N', n, R.data_ptr(), R.desc_ptr(), &ret,
266  &tmp, -1, &liwork, -1, &info);
267 
268  int lwork = (int) tmp;
269  cpuvec<REAL> work(lwork);
270  cpuvec<int> iwork(liwork);
271 
272  scalapack::trcon(norm, 'U', 'N', n, R.data_ptr(), R.desc_ptr(), &ret,
273  x.data_ptr(), lwork, iwork.data_ptr(), iwork.size(), &info);
274  }
275  else
276  {
277  mpimat<REAL> L(x.get_grid(), x.bf_rows(), x.bf_cols());
278  lq(x, aux);
279  lq_L(x, L);
280 
281  REAL tmp;
282  int liwork;
283  scalapack::trcon(norm, 'L', 'N', m, L.data_ptr(), L.desc_ptr(), &ret,
284  &tmp, -1, &liwork, -1, &info);
285 
286  int lwork = (int) tmp;
287  cpuvec<REAL> work(lwork);
288  cpuvec<int> iwork(liwork);
289 
290  scalapack::trcon(norm, 'L', 'N', m, L.data_ptr(), L.desc_ptr(), &ret,
291  x.data_ptr(), x.nrows()*x.ncols(), iwork.data_ptr(), iwork.size(), &info);
292  }
293 
294  fml::linalgutils::check_info(info, "trcon");
295 
296  return ((REAL)1)/ret;
297  }
298  }
299 
319  template <typename REAL>
321  {
322  if (x.is_square())
323  return cond_square_internals('1', x);
324  else
325  return cond_nonsquare_internals('1', x);
326  }
327 
328 
329 
349  template <typename REAL>
351  {
352  if (x.is_square())
353  return cond_square_internals('I', x);
354  else
355  return cond_nonsquare_internals('I', x);
356  }
357 
358 
359 
377  template <typename REAL>
379  {
380  cpuvec<REAL> s;
381  svd(x, s);
382 
383  REAL *s_d = s.data_ptr();
384 
385  REAL max = s_d[0];
386  REAL min = s_d[0];
387  for (len_t i=1; i<s.size(); i++)
388  {
389  if (s_d[i] > max)
390  max = s_d[i];
391  if (s_d[i] > 0 && s_d[i] < min)
392  min = s_d[i];
393  }
394 
395  if (max == 0)
396  return 0;
397  else
398  return max/min;
399  }
400 }
401 }
402 
403 
404 #endif
fml::linalg::norm_I
REAL norm_I(const cpumat< REAL > &x)
Computes the infinity matrix norm, the maximum absolute row sum.
Definition: linalg_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: linalg_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: linalg_norm.hh:281
fml::linalg::qr
void qr(const bool pivot, cpumat< REAL > &x, cpuvec< REAL > &qraux)
Computes the QR decomposition.
Definition: linalg_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: linalg_norm.hh:38
fml::linalg::lu
void lu(cpumat< REAL > &x, cpuvec< int > &p, int &info)
Computes the PLU factorization with partial pivoting.
Definition: linalg_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:514
fml::linalg::cond_I
REAL cond_I(cpumat< REAL > &x)
Estimates the condition number under the infinity norm.
Definition: linalg_norm.hh:311
fml::linalg::svd
void svd(cpumat< REAL > &x, cpuvec< REAL > &s)
Computes the singular value decomposition.
Definition: linalg_svd.hh:97
fml::linalg::norm_2
REAL norm_2(cpumat< REAL > &x)
Computes the 2/spectral matrix norm.
Definition: linalg_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: linalg_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: linalg_qr.hh:162
fml::linalg::lq
void lq(cpumat< REAL > &x, cpuvec< REAL > &lqaux)
Computes the LQ decomposition.
Definition: linalg_qr.hh:223
fml::cpuvec::fill_zero
void fill_zero()
Set all values to zero.
Definition: cpuvec.hh:316
fml::linalg::norm_M
REAL norm_M(const cpumat< REAL > &x)
Computes the maximum modulus matrix norm.
Definition: linalg_norm.hh:145
fml::linalg::norm_F
REAL norm_F(const cpumat< REAL > &x)
Computes the Frobenius/Euclidean matrix norm.
Definition: linalg_norm.hh:118