fml  0.1-0
Fused Matrix Library
det.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_DET_H
6 #define FML_CPU_LINALG_DET_H
7 #pragma once
8 
9 
10 #include <cmath>
11 
12 #include "../cpuvec.hh"
13 #include "../cpumat.hh"
14 
15 #include "lu.hh"
16 
17 
18 namespace fml
19 {
20 namespace linalg
21 {
40  template <typename REAL>
41  void det(cpumat<REAL> &x, int &sign, REAL &modulus)
42  {
43  const len_t m = x.nrows();
44  if (!x.is_square())
45  throw std::runtime_error("'x' must be a square matrix");
46 
47  cpuvec<int> p;
48  int info;
49  lu(x, p, info);
50 
51  if (info != 0)
52  {
53  if (info > 0)
54  {
55  sign = 1;
56  modulus = -INFINITY;
57  return;
58  }
59  else
60  return;
61  }
62 
63 
64  // get determinant
65  REAL mod = 0.0;
66  int sgn = 1;
67 
68  const int *ipiv = p.data_ptr();
69  for (int i=0; i<m; i++)
70  {
71  if (ipiv[i] != (i + 1))
72  sgn = -sgn;
73  }
74 
75  const REAL *a = x.data_ptr();
76 
77  #pragma omp parallel for reduction(+:mod) reduction(*:sgn)
78  for (int i=0; i<m; i++)
79  {
80  const REAL d = a[i + m*i];
81  if (d < 0)
82  {
83  mod += log(-d);
84  sgn *= -1;
85  }
86  else
87  mod += log(d);
88  }
89 
90  modulus = mod;
91  sign = sgn;
92  }
93 }
94 }
95 
96 
97 #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::lu
void lu(cpumat< REAL > &x, cpuvec< int > &p, int &info)
Computes the PLU factorization with partial pivoting.
Definition: lu.hh:48
fml::cpuvec
Vector class for data held on a single CPU.
Definition: cpuvec.hh:31
fml::linalg::det
void det(cpumat< REAL > &x, int &sign, REAL &modulus)
Computes the determinant in logarithmic form.
Definition: det.hh:41
fml::unimat::data_ptr
REAL * data_ptr()
Pointer to the internal array.
Definition: unimat.hh:40
fml
Core namespace.
Definition: dimops.hh:10