fml  0.1-0
Fused Matrix Library
linalg_misc.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_MISC_H
6 #define FML_CPU_LINALG_LINALG_MISC_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 "../cpumat.hh"
17 #include "../cpuvec.hh"
18 
19 #include "lapack.hh"
20 #include "linalg_lu.hh"
21 
22 
23 namespace fml
24 {
25 namespace linalg
26 {
45  template <typename REAL>
46  void det(cpumat<REAL> &x, int &sign, REAL &modulus)
47  {
48  const len_t m = x.nrows();
49  if (!x.is_square())
50  throw std::runtime_error("'x' must be a square matrix");
51 
52  cpuvec<int> p;
53  int info;
54  lu(x, p, info);
55 
56  if (info != 0)
57  {
58  if (info > 0)
59  {
60  sign = 1;
61  modulus = -INFINITY;
62  return;
63  }
64  else
65  return;
66  }
67 
68 
69  // get determinant
70  REAL mod = 0.0;
71  int sgn = 1;
72 
73  const int *ipiv = p.data_ptr();
74  for (int i=0; i<m; i++)
75  {
76  if (ipiv[i] != (i + 1))
77  sgn = -sgn;
78  }
79 
80  const REAL *a = x.data_ptr();
81 
82  #pragma omp parallel for reduction(+:mod) reduction(*:sgn)
83  for (int i=0; i<m; i++)
84  {
85  const REAL d = a[i + m*i];
86  if (d < 0)
87  {
88  mod += log(-d);
89  sgn *= -1;
90  }
91  else
92  mod += log(d);
93  }
94 
95  modulus = mod;
96  sign = sgn;
97  }
98 
99 
100 
108  template <typename REAL>
109  REAL trace(const cpumat<REAL> &x)
110  {
111  const REAL *x_d = x.data_ptr();
112  const len_t m = x.nrows();
113  const len_t minmn = std::min(m, x.ncols());
114 
115  REAL tr = 0;
116  #pragma omp parallel for simd if(minmn > fml::omp::OMP_MIN_SIZE) reduction(+:tr)
117  for (len_t i=0; i<minmn; i++)
118  tr += x_d[i + i*m];
119 
120  return tr;
121  }
122 }
123 }
124 
125 
126 #endif
fml::cpumat
Matrix class for data held on a single CPU.
Definition: cpumat.hh:36
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::trace
REAL trace(const cpumat< REAL > &x)
Computes the trace, i.e. the sum of the diagonal.
Definition: linalg_misc.hh:109
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::unimat::ncols
len_t ncols() const
Number of columns.
Definition: unimat.hh:38
fml::linalg::det
void det(cpumat< REAL > &x, int &sign, REAL &modulus)
Computes the determinant in logarithmic form.
Definition: linalg_misc.hh:46
fml::unimat::data_ptr
REAL * data_ptr()
Pointer to the internal array.
Definition: unimat.hh:40
fml
Core namespace.
Definition: dimops.hh:10