fml  0.1-0
Fused Matrix Library
linalg_qr.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_QR_H
6 #define FML_MPI_LINALG_LINALG_QR_H
7 #pragma once
8 
9 
10 #include <stdexcept>
11 
12 #include "../../_internals/linalgutils.hh"
13 #include "../../cpu/cpuvec.hh"
14 
15 #include "../internals/bcutils.hh"
16 #include "../internals/mpi_utils.hh"
17 
18 #include "../copy.hh"
19 #include "../mpimat.hh"
20 
21 #include "linalg_blas.hh"
22 #include "linalg_err.hh"
23 #include "scalapack.hh"
24 
25 
26 namespace fml
27 {
28 namespace linalg
29 {
30  namespace
31  {
32  template <typename REAL>
33  void qr_internals(const bool pivot, mpimat<REAL> &x, cpuvec<REAL> &qraux, cpuvec<REAL> &work)
34  {
35  const len_t m = x.nrows();
36  const len_t n = x.ncols();
37  const len_t minmn = std::min(m, n);
38 
39  const int *descx = x.desc_ptr();
40 
41  int info = 0;
42  qraux.resize(minmn);
43 
44  REAL tmp;
45  if (pivot)
46  fml::scalapack::geqpf(m, n, NULL, descx, NULL, NULL, &tmp, -1, &info);
47  else
48  fml::scalapack::geqrf(m, n, NULL, descx, NULL, &tmp, -1, &info);
49 
50  int lwork = std::max((int) tmp, 1);
51  if (lwork > work.size())
52  work.resize(lwork);
53 
54  if (pivot)
55  {
56  cpuvec<int> p(n);
57  p.fill_zero();
58  fml::scalapack::geqpf(m, n, x.data_ptr(), descx, p.data_ptr(),
59  qraux.data_ptr(), work.data_ptr(), lwork, &info);
60  }
61  else
62  fml::scalapack::geqrf(m, n, x.data_ptr(), descx, qraux.data_ptr(),
63  work.data_ptr(), lwork, &info);
64 
65  if (info != 0)
66  {
67  if (pivot)
68  fml::linalgutils::check_info(info, "geqpf");
69  else
70  fml::linalgutils::check_info(info, "geqrf");
71  }
72  }
73  }
74 
100  template <typename REAL>
101  void qr(const bool pivot, mpimat<REAL> &x, cpuvec<REAL> &qraux)
102  {
103  cpuvec<REAL> work;
104  qr_internals(pivot, x, qraux, work);
105  }
106 
128  template <typename REAL>
129  void qr_Q(const mpimat<REAL> &QR, const cpuvec<REAL> &qraux, mpimat<REAL> &Q, cpuvec<REAL> &work)
130  {
131  err::check_grid(QR, Q);
132 
133  const len_t m = QR.nrows();
134  const len_t n = QR.ncols();
135  const len_t minmn = std::min(m, n);
136 
137  const int *descQR = QR.desc_ptr();
138 
139  Q.resize(m, minmn);
140  const int *descQ = Q.desc_ptr();
141 
142  int info = 0;
143  REAL tmp;
144  fml::scalapack::orgqr(m, minmn, minmn, NULL, descQR, NULL,
145  &tmp, -1, &info);
146 
147  int lwork = (int) tmp;
148  if (lwork > work.size())
149  work.resize(lwork);
150 
151  fml::scalapack::lacpy('A', m, minmn, QR.data_ptr(), descQR, Q.data_ptr(),
152  descQ);
153 
154  fml::scalapack::orgqr(m, minmn, minmn, Q.data_ptr(), descQR,
155  qraux.data_ptr(), work.data_ptr(), lwork, &info);
156  fml::linalgutils::check_info(info, "orgqr");
157  }
158 
178  template <typename REAL>
179  void qr_R(const mpimat<REAL> &QR, mpimat<REAL> &R)
180  {
181  err::check_grid(QR, R);
182 
183  const len_t m = QR.nrows();
184  const len_t n = QR.ncols();
185  const len_t minmn = std::min(m, n);
186 
187  R.resize(minmn, n);
188  R.fill_zero();
189  fml::scalapack::lacpy('U', m, n, QR.data_ptr(), QR.desc_ptr(), R.data_ptr(),
190  R.desc_ptr());
191  }
192 
193 
194 
195  namespace
196  {
197  template <typename REAL>
198  void lq_internals(mpimat<REAL> &x, cpuvec<REAL> &lqaux, cpuvec<REAL> &work)
199  {
200  const len_t m = x.nrows();
201  const len_t n = x.ncols();
202  const len_t minmn = std::min(m, n);
203 
204  const int *descx = x.desc_ptr();
205 
206  int info = 0;
207  lqaux.resize(minmn);
208 
209  REAL tmp;
210  fml::scalapack::gelqf(m, n, NULL, descx, NULL, &tmp, -1, &info);
211  int lwork = std::max((int) tmp, 1);
212  if (lwork > work.size())
213  work.resize(lwork);
214 
215  fml::scalapack::gelqf(m, n, x.data_ptr(), descx, lqaux.data_ptr(),
216  work.data_ptr(), lwork, &info);
217 
218  if (info != 0)
219  fml::linalgutils::check_info(info, "gelqf");
220  }
221  }
222 
246  template <typename REAL>
247  void lq(mpimat<REAL> &x, cpuvec<REAL> &lqaux)
248  {
249  cpuvec<REAL> work;
250  lq_internals(x, lqaux, work);
251  }
252 
272  template <typename REAL>
273  void lq_L(const mpimat<REAL> &LQ, mpimat<REAL> &L)
274  {
275  err::check_grid(LQ, L);
276 
277  const len_t m = LQ.nrows();
278  const len_t n = LQ.ncols();
279  const len_t minmn = std::min(m, n);
280 
281  L.resize(m, minmn);
282  L.fill_zero();
283 
284  fml::scalapack::lacpy('L', m, n, LQ.data_ptr(), LQ.desc_ptr(), L.data_ptr(),
285  L.desc_ptr());
286  }
287 
309  template <typename REAL>
310  void lq_Q(const mpimat<REAL> &LQ, const cpuvec<REAL> &lqaux, mpimat<REAL> &Q, cpuvec<REAL> &work)
311  {
312  err::check_grid(LQ, Q);
313 
314  const len_t m = LQ.nrows();
315  const len_t n = LQ.ncols();
316  const len_t minmn = std::min(m, n);
317 
318  const int *descLQ = LQ.desc_ptr();
319 
320  Q.resize(minmn, n);
321  const int *descQ = Q.desc_ptr();
322 
323  int info = 0;
324  REAL tmp;
325  fml::scalapack::orglq(minmn, n, minmn, NULL, descLQ, NULL,
326  &tmp, -1, &info);
327 
328  int lwork = (int) tmp;
329  if (lwork > work.size())
330  work.resize(lwork);
331 
332  fml::scalapack::lacpy('A', minmn, n, LQ.data_ptr(), descLQ, Q.data_ptr(),
333  descQ);
334 
335  fml::scalapack::orglq(minmn, n, minmn, Q.data_ptr(), descQ,
336  lqaux.data_ptr(), work.data_ptr(), lwork, &info);
337  fml::linalgutils::check_info(info, "orglq");
338  }
339 }
340 }
341 
342 
343 #endif
fml::mpimat::resize
void resize(len_t nrows, len_t ncols)
Resize the internal object storage.
Definition: mpimat.hh:323
fml::mpimat
Matrix class for data distributed over MPI in the 2-d block cyclic format.
Definition: mpimat.hh:40
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:209
fml::linalg::qr
void qr(const bool pivot, cpumat< REAL > &x, cpuvec< REAL > &qraux)
Computes the QR decomposition.
Definition: linalg_qr.hh:94
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::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::mpimat::fill_zero
void fill_zero()
Set all values to zero.
Definition: mpimat.hh:562
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_qr.hh:120
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::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_qr.hh:279