fml  0.1-0
Fused Matrix Library
linalg_factorizations.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_FACTORIZATIONS_H
6 #define FML_CPU_LINALG_LINALG_FACTORIZATIONS_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 "../internals/cpu_utils.hh"
17 
18 #include "../copy.hh"
19 #include "../cpumat.hh"
20 #include "../cpuvec.hh"
21 
22 #include "lapack.hh"
23 #include "linalg_blas.hh"
24 
25 
26 namespace fml
27 {
29 namespace linalg
30 {
52  template <typename REAL>
53  void lu(cpumat<REAL> &x, cpuvec<int> &p, int &info)
54  {
55  info = 0;
56  const len_t m = x.nrows();
57  const len_t lipiv = std::min(m, x.ncols());
58 
59  p.resize(lipiv);
60 
61  fml::lapack::getrf(m, x.ncols(), x.data_ptr(), m, p.data_ptr(), &info);
62  }
63 
65  template <typename REAL>
66  void lu(cpumat<REAL> &x)
67  {
68  cpuvec<int> p;
69  int info;
70 
71  lu(x, p, info);
72 
73  fml::linalgutils::check_info(info, "getrf");
74  }
75 
76 
77 
78  namespace
79  {
80  template <typename REAL>
81  int svd_internals(const int nu, const int nv, cpumat<REAL> &x, cpuvec<REAL> &s,
82  cpumat<REAL> &u, cpumat<REAL> &vt)
83  {
84  int info = 0;
85  char jobz;
86  int ldvt;
87 
88  const len_t m = x.nrows();
89  const len_t n = x.ncols();
90  const len_t minmn = std::min(m, n);
91 
92  s.resize(minmn);
93 
94  if (nu == 0 && nv == 0)
95  {
96  jobz = 'N';
97  ldvt = 1; // value is irrelevant, but must exist!
98  }
99  else if (nu <= minmn && nv <= minmn)
100  {
101  jobz = 'S';
102  ldvt = minmn;
103 
104  u.resize(m, minmn);
105  vt.resize(minmn, n);
106  }
107  else
108  {
109  jobz = 'A';
110  ldvt = n;
111  }
112 
113  cpuvec<int> iwork(8*minmn);
114 
115  REAL tmp;
116  fml::lapack::gesdd(jobz, m, n, x.data_ptr(), m, s.data_ptr(), u.data_ptr(), m, vt.data_ptr(), ldvt, &tmp, -1, iwork.data_ptr(), &info);
117  int lwork = (int) tmp;
118  cpuvec<REAL> work(lwork);
119 
120  fml::lapack::gesdd(jobz, m, n, x.data_ptr(), m, s.data_ptr(), u.data_ptr(), m, vt.data_ptr(), ldvt, work.data_ptr(), lwork, iwork.data_ptr(), &info);
121 
122  return info;
123  }
124  }
125 
145  template <typename REAL>
147  {
148  cpumat<REAL> ignored;
149  int info = svd_internals(0, 0, x, s, ignored, ignored);
150  fml::linalgutils::check_info(info, "gesdd");
151  }
152 
154  template <typename REAL>
156  {
157  int info = svd_internals(1, 1, x, s, u, vt);
158  fml::linalgutils::check_info(info, "gesdd");
159  }
160 
161 
162 
163  namespace
164  {
165  template <typename REAL>
166  int eig_sym_internals(const bool only_values, cpumat<REAL> &x,
167  cpuvec<REAL> &values, cpumat<REAL> &vectors)
168  {
169  if (!x.is_square())
170  throw std::runtime_error("'x' must be a square matrix");
171 
172  int info = 0;
173  int nfound;
174  char jobz;
175 
176  len_t n = x.nrows();
177  values.resize(n);
178  cpuvec<int> support;
179 
180  if (only_values)
181  jobz = 'N';
182  else
183  {
184  jobz = 'V';
185  vectors.resize(n, n);
186  support.resize(2*n);
187  }
188 
189  REAL worksize;
190  int lwork, liwork;
191  fml::lapack::syevr(jobz, 'A', 'L', n, x.data_ptr(), n, (REAL) 0.f, (REAL) 0.f,
192  0, 0, (REAL) 0.f, &nfound, values.data_ptr(), vectors.data_ptr(), n,
193  support.data_ptr(), &worksize, -1, &liwork, -1,
194  &info);
195 
196  lwork = (int) worksize;
197  cpuvec<REAL> work(lwork);
198  cpuvec<int> iwork(liwork);
199 
200  fml::lapack::syevr(jobz, 'A', 'L', n, x.data_ptr(), n, (REAL) 0.f, (REAL) 0.f,
201  0, 0, (REAL) 0.f, &nfound, values.data_ptr(), vectors.data_ptr(), n,
202  support.data_ptr(), work.data_ptr(), lwork, iwork.data_ptr(), liwork,
203  &info);
204 
205  return info;
206  }
207  }
208 
229  template <typename REAL>
231  {
232  cpumat<REAL> ignored;
233 
234  int info = eig_sym_internals(true, x, values, ignored);
235  fml::linalgutils::check_info(info, "syevr");
236  }
237 
239  template <typename REAL>
240  void eigen_sym(cpumat<REAL> &x, cpuvec<REAL> &values, cpumat<REAL> &vectors)
241  {
242  int info = eig_sym_internals(false, x, values, vectors);
243  fml::linalgutils::check_info(info, "syevr");
244  }
245 
246 
247 
264  template <typename REAL>
266  {
267  const len_t n = x.nrows();
268  if (!x.is_square())
269  throw std::runtime_error("'x' must be a square matrix");
270 
271  // Factor x = LU
272  cpuvec<int> p;
273  int info;
274  lu(x, p, info);
275  fml::linalgutils::check_info(info, "getrf");
276 
277  // Invert
278  REAL tmp;
279  fml::lapack::getri(n, x.data_ptr(), n, p.data_ptr(), &tmp, -1, &info);
280  int lwork = (int) tmp;
281  cpuvec<REAL> work(lwork);
282 
283  fml::lapack::getri(n, x.data_ptr(), n, p.data_ptr(), work.data_ptr(), lwork, &info);
284  fml::linalgutils::check_info(info, "getri");
285  }
286 
287 
288 
289  namespace
290  {
291  template <typename REAL>
292  void solver(cpumat<REAL> &x, len_t ylen, len_t nrhs, REAL *y_d)
293  {
294  const len_t n = x.nrows();
295  if (!x.is_square())
296  throw std::runtime_error("'x' must be a square matrix");
297  if (n != ylen)
298  throw std::runtime_error("rhs 'y' must be compatible with data matrix 'x'");
299 
300  int info;
301  cpuvec<int> p(n);
302  fml::lapack::gesv(n, nrhs, x.data_ptr(), n, p.data_ptr(), y_d, n, &info);
303  fml::linalgutils::check_info(info, "gesv");
304  }
305  }
306 
325  template <typename REAL>
327  {
328  solver(x, y.size(), 1, y.data_ptr());
329  }
330 
332  template <typename REAL>
334  {
335  solver(x, y.nrows(), y.ncols(), y.data_ptr());
336  }
337 
338 
339 
340  namespace
341  {
342  template <typename REAL>
343  void qr_internals(const bool pivot, cpumat<REAL> &x, cpuvec<REAL> &qraux, cpuvec<REAL> &work)
344  {
345  const len_t m = x.nrows();
346  const len_t n = x.ncols();
347  const len_t minmn = std::min(m, n);
348 
349  int info = 0;
350  qraux.resize(minmn);
351 
352  REAL tmp;
353  if (pivot)
354  fml::lapack::geqp3(m, n, NULL, m, NULL, NULL, &tmp, -1, &info);
355  else
356  fml::lapack::geqrf(m, n, NULL, m, NULL, &tmp, -1, &info);
357 
358  int lwork = std::max((int) tmp, 1);
359  if (lwork > work.size())
360  work.resize(lwork);
361 
362  if (pivot)
363  {
364  cpuvec<int> p(n);
365  p.fill_zero();
366  fml::lapack::geqp3(m, n, x.data_ptr(), m, p.data_ptr(), qraux.data_ptr(), work.data_ptr(), lwork, &info);
367  }
368  else
369  fml::lapack::geqrf(m, n, x.data_ptr(), m, qraux.data_ptr(), work.data_ptr(), lwork, &info);
370 
371  if (info != 0)
372  {
373  if (pivot)
374  fml::linalgutils::check_info(info, "geqp3");
375  else
376  fml::linalgutils::check_info(info, "geqrf");
377  }
378  }
379  }
380 
404  template <typename REAL>
405  void qr(const bool pivot, cpumat<REAL> &x, cpuvec<REAL> &qraux)
406  {
407  cpuvec<REAL> work;
408  qr_internals(pivot, x, qraux, work);
409  }
410 
430  template <typename REAL>
431  void qr_Q(const cpumat<REAL> &QR, const cpuvec<REAL> &qraux, cpumat<REAL> &Q,
432  cpuvec<REAL> &work)
433  {
434  const len_t m = QR.nrows();
435  const len_t n = QR.ncols();
436  const len_t minmn = std::min(m, n);
437 
438  int info = 0;
439  REAL tmp;
440  fml::lapack::ormqr('L', 'N', m, minmn, m, QR.data_ptr(), m, NULL,
441  NULL, m, &tmp, -1, &info);
442 
443  int lwork = (int) tmp;
444  if (lwork > work.size())
445  work.resize(lwork);
446 
447  Q.resize(m, minmn);
448  Q.fill_eye();
449 
450  fml::lapack::ormqr('L', 'N', m, minmn, m, QR.data_ptr(), m, qraux.data_ptr(),
451  Q.data_ptr(), m, work.data_ptr(), lwork, &info);
452  fml::linalgutils::check_info(info, "ormqr");
453  }
454 
472  template <typename REAL>
473  void qr_R(const cpumat<REAL> &QR, cpumat<REAL> &R)
474  {
475  const len_t m = QR.nrows();
476  const len_t n = QR.ncols();
477  const len_t minmn = std::min(m, n);
478 
479  R.resize(minmn, n);
480  R.fill_zero();
481  fml::lapack::lacpy('U', m, n, QR.data_ptr(), m, R.data_ptr(), minmn);
482  }
483 
484 
485 
486  namespace
487  {
488  template <typename REAL>
489  void lq_internals(cpumat<REAL> &x, cpuvec<REAL> &lqaux, cpuvec<REAL> &work)
490  {
491  const len_t m = x.nrows();
492  const len_t n = x.ncols();
493  const len_t minmn = std::min(m, n);
494 
495  int info = 0;
496  lqaux.resize(minmn);
497 
498  REAL tmp;
499  fml::lapack::gelqf(m, n, NULL, m, NULL, &tmp, -1, &info);
500  int lwork = std::max((int) tmp, 1);
501  if (lwork > work.size())
502  work.resize(lwork);
503 
504  fml::lapack::gelqf(m, n, x.data_ptr(), m, lqaux.data_ptr(), work.data_ptr(),
505  lwork, &info);
506 
507  if (info != 0)
508  fml::linalgutils::check_info(info, "gelqf");
509  }
510  }
511 
533  template <typename REAL>
534  void lq(cpumat<REAL> &x, cpuvec<REAL> &lqaux)
535  {
536  cpuvec<REAL> work;
537  lq_internals(x, lqaux, work);
538  }
539 
557  template <typename REAL>
558  void lq_L(const cpumat<REAL> &LQ, cpumat<REAL> &L)
559  {
560  const len_t m = LQ.nrows();
561  const len_t n = LQ.ncols();
562  const len_t minmn = std::min(m, n);
563 
564  L.resize(m, minmn);
565  L.fill_zero();
566 
567  fml::lapack::lacpy('L', m, n, LQ.data_ptr(), m, L.data_ptr(), m);
568  }
569 
589  template <typename REAL>
590  void lq_Q(const cpumat<REAL> &LQ, const cpuvec<REAL> &lqaux, cpumat<REAL> &Q,
591  cpuvec<REAL> &work)
592  {
593  const len_t m = LQ.nrows();
594  const len_t n = LQ.ncols();
595  const len_t minmn = std::min(m, n);
596 
597  int info = 0;
598  REAL tmp;
599  fml::lapack::ormlq('R', 'N', minmn, n, minmn, LQ.data_ptr(), m, NULL,
600  NULL, minmn, &tmp, -1, &info);
601 
602  int lwork = (int) tmp;
603  if (lwork > work.size())
604  work.resize(lwork);
605 
606  Q.resize(minmn, n);
607  Q.fill_eye();
608 
609  fml::lapack::ormlq('R', 'N', minmn, n, minmn, LQ.data_ptr(), m, lqaux.data_ptr(),
610  Q.data_ptr(), minmn, work.data_ptr(), lwork, &info);
611  fml::linalgutils::check_info(info, "ormlq");
612  }
613 
614 
615 
635  template <typename REAL>
637  {
638  const len_t n = x.nrows();
639  if (n != x.ncols())
640  throw std::runtime_error("'x' must be a square matrix");
641 
642  int info = 0;
643  fml::lapack::potrf('L', n, x.data_ptr(), n, &info);
644 
645  if (info < 0)
646  fml::linalgutils::check_info(info, "potrf");
647  else if (info > 0)
648  throw std::runtime_error("chol: leading minor of order " + std::to_string(info) + " is not positive definite");
649 
650  fml::cpu_utils::tri2zero('U', false, n, n, x.data_ptr(), n);
651  }
652 }
653 }
654 
655 
656 #endif
fml::cpumat
Matrix class for data held on a single CPU.
Definition: cpumat.hh:36
fml::cpumat::fill_eye
void fill_eye()
Set diagonal entries to 1 and non-diagonal entries to 0.
Definition: cpumat.hh:401
fml::linalg::solve
void solve(cpumat< REAL > &x, cpuvec< REAL > &y)
Solve a system of equations.
Definition: linalg_factorizations.hh:326
fml::unimat::is_square
bool is_square() const
Is the matrix square?
Definition: unimat.hh:34
fml::cpumat::fill_zero
void fill_zero()
Set all values to zero.
Definition: cpumat.hh:321
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::cpuvec::resize
void resize(len_t size)
Resize the internal object storage.
Definition: cpuvec.hh:185
fml::linalg::qr
void qr(const bool pivot, cpumat< REAL > &x, cpuvec< REAL > &qraux)
Computes the QR decomposition.
Definition: linalg_factorizations.hh:405
fml::linalg::lu
void lu(cpumat< REAL > &x, cpuvec< int > &p, int &info)
Computes the PLU factorization with partial pivoting.
Definition: linalg_factorizations.hh:53
fml::cpumat::resize
void resize(len_t nrows, len_t ncols)
Resize the internal object storage.
Definition: cpumat.hh:207
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::linalg::svd
void svd(cpumat< REAL > &x, cpuvec< REAL > &s)
Computes the singular value decomposition.
Definition: linalg_factorizations.hh:146
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_factorizations.hh:558
fml
Core namespace.
Definition: linalgutils.hh:15
fml::linalg::qr_Q
void qr_Q(const cpumat< REAL > &QR, const cpuvec< REAL > &qraux, cpumat< REAL > &Q, cpuvec< REAL > &work)
Recover the Q matrix from a QR decomposition.
Definition: linalg_factorizations.hh:431
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_factorizations.hh:473
fml::linalg::chol
void chol(cpumat< REAL > &x)
Compute the Choleski factorization.
Definition: linalg_factorizations.hh:636
fml::linalg::lq
void lq(cpumat< REAL > &x, cpuvec< REAL > &lqaux)
Computes the LQ decomposition.
Definition: linalg_factorizations.hh:534
fml::linalg::eigen_sym
void eigen_sym(cpumat< REAL > &x, cpuvec< REAL > &values)
Compute the eigenvalues and optionally the eigenvectors for a symmetric matrix.
Definition: linalg_factorizations.hh:230
fml::linalg::invert
void invert(cpumat< REAL > &x)
Compute the matrix inverse.
Definition: linalg_factorizations.hh:265
fml::linalg::lq_Q
void lq_Q(const cpumat< REAL > &LQ, const cpuvec< REAL > &lqaux, cpumat< REAL > &Q, cpuvec< REAL > &work)
Recover the Q matrix from an LQ decomposition.
Definition: linalg_factorizations.hh:590