AxioDL/metaforce
2
0
Fork 0
mirror of https://github.com/AxioDL/metaforce.git synced 2025-12-08 18:24:55 +00:00
Code Issues Projects Releases Wiki Activity

Windows fixes

Browse Source
This commit is contained in:
Jack Andersen
2018-10-14 10:16:21 -10:00
parent 1559163f98
commit 1851308021
68 changed files with 1343 additions and 968 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
gmm/gmm_dense_lu.h
View File

@@ -70,10 +70,50 @@
#define GMM_DENSE_LU_H
#include "gmm_dense_Householder.h"
#include "gmm_opt.h"
namespace gmm {
/* ********************************************************************** */
/* IPVT structure. */
/* ********************************************************************** */
// For compatibility with lapack version with 64 or 32 bit integer.
// Should be replaced by std::vector<size_type> if 32 bit integer version
// of lapack is not used anymore (and lapack_ipvt_int set to size_type)
// Do not use iterators of this interface container
class lapack_ipvt : public std::vector<size_type> {
bool is_int64;
size_type &operator[](size_type i)
{ return std::vector<size_type>::operator[](i); }
size_type operator[] (size_type i) const
{ return std::vector<size_type>::operator[](i); }
void begin(void) const {}
void begin(void) {}
void end(void) const {}
void end(void) {}
public:
void set_to_int32() { is_int64 = false; }
const size_type *pfirst() const
{ return &(*(std::vector<size_type>::begin())); }
size_type *pfirst() { return &(*(std::vector<size_type>::begin())); }
lapack_ipvt(size_type n) : std::vector<size_type>(n), is_int64(true) {}
size_type get(size_type i) const {
const size_type *p = pfirst();
return is_int64 ? p[i] : size_type(((const int *)(p))[i]);
}
void set(size_type i, size_type val) {
size_type *p = pfirst();
if (is_int64) p[i] = val; else ((int *)(p))[i] = int(val);
}
};
}
#include "gmm_opt.h"
namespace gmm {
/** LU Factorization of a general (dense) matrix (real or complex).
@@ -85,24 +125,23 @@ namespace gmm {
The pivot indices in ipvt are indexed starting from 1
so that this is compatible with LAPACK (Fortran).
*/
template <typename DenseMatrix, typename Pvector>
size_type lu_factor(DenseMatrix& A, Pvector& ipvt) {
template <typename DenseMatrix>
size_type lu_factor(DenseMatrix& A, lapack_ipvt& ipvt) {
typedef typename linalg_traits<DenseMatrix>::value_type T;
typedef typename linalg_traits<Pvector>::value_type int_T;
typedef typename number_traits<T>::magnitude_type R;
size_type info(0), i, j, jp, M(mat_nrows(A)), N(mat_ncols(A));
size_type NN = std::min(M, N);
std::vector<T> c(M), r(N);
GMM_ASSERT2(ipvt.size()+1 >= NN, "IPVT too small");
for (i = 0; i+1 < NN; ++i) ipvt[i] = int_T(i);
for (i = 0; i+1 < NN; ++i) ipvt.set(i, i);
if (M || N) {
for (j = 0; j+1 < NN; ++j) {
R max = gmm::abs(A(j,j)); jp = j;
for (i = j+1; i < M; ++i) /* find pivot. */
if (gmm::abs(A(i,j)) > max) { jp = i; max = gmm::abs(A(i,j)); }
ipvt[j] = int_T(jp + 1);
ipvt.set(j, jp + 1);
if (max == R(0)) { info = j + 1; break; }
if (jp != j) for (i = 0; i < N; ++i) std::swap(A(jp, i), A(j, i));
@@ -112,7 +151,7 @@ namespace gmm {
rank_one_update(sub_matrix(A, sub_interval(j+1, M-j-1),
sub_interval(j+1, N-j-1)), c, conjugated(r));
}
ipvt[NN-1] = int_T(NN);
ipvt.set(NN-1, NN);
}
return info;
}
@@ -126,7 +165,7 @@ namespace gmm {
typedef typename linalg_traits<DenseMatrix>::value_type T;
copy(b, x);
for(size_type i = 0; i < pvector.size(); ++i) {
size_type perm = pvector[i]-1; // permutations stored in 1's offset
size_type perm = pvector.get(i)-1; // permutations stored in 1's offset
if(i != perm) { T aux = x[i]; x[i] = x[perm]; x[perm] = aux; }
}
/* solve Ax = b -> LUx = b -> Ux = L^-1 b. */
@@ -138,7 +177,7 @@ namespace gmm {
void lu_solve(const DenseMatrix &A, VectorX &x, const VectorB &b) {
typedef typename linalg_traits<DenseMatrix>::value_type T;
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
std::vector<int> ipvt(mat_nrows(A));
lapack_ipvt ipvt(mat_nrows(A));
gmm::copy(A, B);
size_type info = lu_factor(B, ipvt);
GMM_ASSERT1(!info, "Singular system, pivot = " << info);
@@ -154,10 +193,9 @@ namespace gmm {
lower_tri_solve(transposed(LU), x, false);
upper_tri_solve(transposed(LU), x, true);
for(size_type i = pvector.size(); i > 0; --i) {
size_type perm = pvector[i-1]-1; // permutations stored in 1's offset
size_type perm = pvector.get(i-1)-1; // permutations stored in 1's offset
if(i-1 != perm) { T aux = x[i-1]; x[i-1] = x[perm]; x[perm] = aux; }
}
}
@@ -212,7 +250,7 @@ namespace gmm {
typedef typename linalg_traits<DenseMatrix>::value_type T;
DenseMatrix& A = const_cast<DenseMatrix&>(A_);
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
std::vector<int> ipvt(mat_nrows(A));
lapack_ipvt ipvt(mat_nrows(A));
gmm::copy(A, B);
size_type info = lu_factor(B, ipvt);
if (doassert) GMM_ASSERT1(!info, "Non invertible matrix, pivot = "<<info);
@@ -229,7 +267,7 @@ namespace gmm {
for (size_type j = 0; j < std::min(mat_nrows(LU), mat_ncols(LU)); ++j)
det *= LU(j,j);
for(size_type i = 0; i < pvector.size(); ++i)
if (i != size_type(pvector[i]-1)) { det = -det; }
if (i != size_type(pvector.get(i)-1)) { det = -det; }
return det;
}
@@ -238,7 +276,7 @@ namespace gmm {
lu_det(const DenseMatrix& A) {
typedef typename linalg_traits<DenseMatrix>::value_type T;
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
std::vector<int> ipvt(mat_nrows(A));
lapack_ipvt ipvt(mat_nrows(A));
gmm::copy(A, B);
lu_factor(B, ipvt);
return lu_det(B, ipvt);