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_MPI_LINALG_DET_H
6 #define FML_MPI_LINALG_DET_H
7 #pragma once
8 
9 
10 #include <stdexcept>
11 
12 #include "../../cpu/cpuvec.hh"
13 
14 #include "../mpimat.hh"
15 
16 #include "lu.hh"
17 
18 
19 namespace fml
20 {
21 namespace linalg
22 {
41  template <typename REAL>
42  void det(mpimat<REAL> &x, int &sign, REAL &modulus)
43  {
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 len_t m_local = x.nrows_local();
69  const len_t n_local = x.ncols_local();
70 
71  const int *ipiv = p.data_ptr();
72  const REAL *a = x.data_ptr();
73  const grid g = x.get_grid();
74 
75  for (len_t i=0; i<m_local; i++)
76  {
77  len_t gi = fml::bcutils::l2g(i, x.bf_rows(), g.nprow(), g.myrow());
78 
79  if (ipiv[i] != (gi + 1))
80  sgn = -sgn;
81  }
82 
83  for (len_t j=0; j<n_local; j++)
84  {
85  for (len_t i=0; i<m_local; i++)
86  {
87  len_t gi = fml::bcutils::l2g(i, x.bf_rows(), g.nprow(), g.myrow());
88  len_t gj = fml::bcutils::l2g(j, x.bf_cols(), g.npcol(), g.mycol());
89 
90  if (gi == gj)
91  {
92  const REAL d = a[i + m_local*j];
93  if (d < 0)
94  {
95  mod += log(-d);
96  sgn *= -1;
97  }
98  else
99  mod += log(d);
100  }
101  }
102  }
103 
104  g.allreduce(1, 1, &mod);
105 
106  sgn = (sgn<0 ? 1 : 0);
107  g.allreduce(1, 1, &sgn, 'C');
108  sgn = (sgn%2==0 ? 1 : -1);
109 
110  modulus = mod;
111  sign = sgn;
112  }
113 }
114 }
115 
116 
117 #endif
fml::grid::allreduce
void allreduce(const int m, const int n, int *x, const char scope='A', const blacsops op=BLACS_SUM) const
Sum reduce operation across all processes in the grid.
Definition: grid.hh:420
fml::grid
2-dimensional MPI process grid.
Definition: grid.hh:70
fml::mpimat
Matrix class for data distributed over MPI in the 2-d block cyclic format.
Definition: mpimat.hh:40
fml::grid::mycol
int mycol() const
The process column (0-based index) of the calling process.
Definition: grid.hh:129
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::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
fml::grid::npcol
int npcol() const
The number of processes columns in the BLACS context.
Definition: grid.hh:125
fml::grid::myrow
int myrow() const
The process row (0-based index) of the calling process.
Definition: grid.hh:127
fml::grid::nprow
int nprow() const
The number of processes rows in the BLACS context.
Definition: grid.hh:123