 |
fml
0.1-0
Fused Matrix Library
|
5 #ifndef FML_PAR_CPU_LINALG_SVD_H
6 #define FML_PAR_CPU_LINALG_SVD_H
10 #include "../parmat.hh"
15 #include "../../../_internals/omp.hh"
17 #include "../../../cpu/linalg/linalg_blas.hh"
18 #include "../../../cpu/linalg/linalg_invert.hh"
19 #include "../../../cpu/linalg/linalg_qr.hh"
20 #include "../../../cpu/linalg/linalg_svd.hh"
22 #include "../../../cpu/copy.hh"
53 template <
typename REAL>
56 const len_t n = x.ncols();
58 auto cp =
crossprod((REAL)1.0, x.data_obj());
59 x.get_comm().allreduce(n*n, cp.data_ptr());
64 template <
typename REAL>
68 const len_t n = x.ncols();
70 auto cp =
crossprod((REAL)1.0, x.data_obj());
71 x.get_comm().allreduce(n*n, cp.data_ptr());
77 for (len_t i=0; i<s.
size(); i++)
78 s_d[i] = sqrt(fabs(s_d[i]));
83 #pragma omp parallel for if(n*n > omp::OMP_MIN_SIZE)
84 for (len_t j=0; j<n; j++)
87 for (len_t i=0; i<n; i++)
88 vt_d[i + n*j] /= s_d[j];
129 template <
typename REAL>
133 qr_R(mpi::REDUCE_TO_ALL, x, R);
138 template <
typename REAL>
142 const len_t n = x.ncols();
143 if (x.nrows() < (len_global_t)n)
144 throw std::runtime_error(
"impossible dimensions");
150 qr(
false, x.data_obj(), qraux);
151 qr_R(x.data_obj(), R_local);
154 tsqr::qr_allreduce(mpi::REDUCE_TO_ALL, n, n, R_local.
data_ptr(),
155 R.
data_ptr(), x.get_comm().get_comm());
162 u.resize(x.nrows(), x.ncols());
165 linalg::matmult(
false,
false, (REAL)1, x.data_obj(), u_R, u.data_obj());
172 template <
typename REAL>
173 void rsvd_A(
const uint32_t seed,
const int k,
const int q,
176 const len_global_t m = x.nrows();
177 const len_t n = x.ncols();
178 if (m < (len_global_t)n)
179 throw std::runtime_error(
"must have m>n");
181 throw std::runtime_error(
"must have k<n");
194 omega.fill_runif(seed);
198 qr_Q(Y, Y_tmp, R, R_local, qraux, QY);
200 for (
int i=0; i<q; i++)
203 linalg::qr_internals(
false, Z, qraux, work);
207 qr_Q(Y, Y_tmp, R, R_local, qraux, QY);
234 template <
typename REAL>
235 void rsvd(
const uint32_t seed,
const int k,
const int q,
242 internals::rsvd_A(seed, k, q, x, QY);
252 template <
typename REAL>
253 void rsvd(
const uint32_t seed,
const int k,
const int q,
261 internals::rsvd_A(seed, k, q, x, QY);
Matrix class for data held on a single CPU.
Definition: cpumat.hh:36
void cpu2cpu(const cpuvec< REAL_IN > &cpu_in, cpuvec< REAL_OUT > &cpu_out)
Copy data from a CPU object to another.
Definition: copy.hh:40
void crossprod(const REAL alpha, const cpumat< REAL > &x, cpumat< REAL > &ret)
Computes lower triangle of alpha*x^T*x.
Definition: crossprod.hh:37
void xpose(const cpumat< REAL > &x, cpumat< REAL > &tx)
Computes the transpose out-of-place (i.e. in a copy).
Definition: xpose.hh:37
T * data_ptr()
Pointer to the internal array.
Definition: univec.hh:28
void resize(len_t size)
Resize the internal object storage.
Definition: cpuvec.hh:210
void qr(const bool pivot, cpumat< REAL > &x, cpuvec< REAL > &qraux)
Computes the QR decomposition.
Definition: qr.hh:94
void rsvd(const uint32_t seed, const int k, const int q, cpumat< REAL > &x, cpuvec< REAL > &s)
Computes the truncated singular value decomposition using the normal projections method of Halko et a...
Definition: svd.hh:462
void rev()
Reverse the vector.
Definition: cpuvec.hh:447
Vector class for data held on a single CPU.
Definition: cpuvec.hh:31
void tssvd(parmat_cpu< REAL > &x, cpuvec< REAL > &s)
Computes the singular value decomposition by first reducing the rectangular matrix to a square matrix...
Definition: svd.hh:130
void svd(cpumat< REAL > &x, cpuvec< REAL > &s)
Computes the singular value decomposition.
Definition: svd.hh:101
REAL * data_ptr()
Pointer to the internal array.
Definition: unimat.hh:40
void rev_cols()
Reverse the columns of the matrix.
Definition: cpumat.hh:638
Core namespace.
Definition: dimops.hh:10
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: qr.hh:120
void trinv(const bool upper, const bool unit_diag, cpumat< REAL > &x)
Compute the matrix inverse of a triangular matrix.
Definition: invert.hh:87
len_t size() const
Number of elements in the vector.
Definition: univec.hh:26
void qr_R(const cpumat< REAL > &QR, cpumat< REAL > &R)
Recover the R matrix from a QR decomposition.
Definition: qr.hh:162
void eigen_sym(cpumat< REAL > &x, cpuvec< REAL > &values)
Compute the eigenvalues and optionally the eigenvectors for a symmetric matrix.
Definition: eigen.hh:95
void matmult(const bool transx, const bool transy, const REAL alpha, const cpumat< REAL > &x, const cpumat< REAL > &y, cpumat< REAL > &ret)
Computes ret = alpha*op(x)*op(y) where op(A) is A or A^T.
Definition: matmult.hh:43
void cpsvd(const cpumat< REAL > &x, cpuvec< REAL > &s, cpumat< REAL > &u, cpumat< REAL > &vt)
Computes the singular value decomposition using the "crossproducts SVD". This method is not numerical...
Definition: svd.hh:323