Commit 94968414 authored by Lucas Serrano's avatar Lucas Serrano
Browse files

Remove old files

parent 6bf99bd1
GTEST_DIR := /usr/src/googletest/googletest
CPATH := ${BOOST_SIMD_ROOT}
ifneq ($(CPATH),)
CPATH := -I $(CPATH)
endif
CXX ?= g++
CXXFLAGS = -O3 -march=native
all:
lib/gtest-all.o: $(GTEST_DIR)/src/gtest-all.cc
g++ -isystem $(GTEST_DIR)/include -I $(GTEST_DIR) -pthread -c $< -o $@
obj/unit_%.o: tests/unit_%.cpp fast_matrix.hpp
$(CXX) $(CPATH) $(CXXFLAGS) -isystem $(GTEST_DIR)/include -pthread $< -c -o $@
test: bin/tests
bin/tests: tests/run_all.cpp $(patsubst tests/unit_%.cpp,obj/unit_%.o,$(wildcard tests/unit_*.cpp)) fast_matrix.hpp lib/gtest-all.o
$(CXX) $(CPATH) $(CXXFLAGS) -isystem $(GTEST_DIR)/include -pthread $^ -o $@
$@
clean:
rm -f bin/*
rm -f lib/*
rm -f obj/*
#include<iostream>
#include<cstring>
#include "../dsl/fast_matrix.hpp"
#ifndef SIZE
#define SIZE 4
#endif
#ifndef VECSIZE
#define VECSIZE 8
#endif
template <typename M>
void chol_inv(M &a, M &i);
template<typename T, int MatrixSize, int MaxVecSize>
class InvertibleMatrix : public BaseMatrix<T, MatrixSize, MaxVecSize> {
protected:
using matrix_t = InvertibleMatrix<T, MatrixSize, MaxVecSize>;
using ElementType = T;
using pack_t = typename BaseMatrix<T, MatrixSize, MaxVecSize>::pack_t;
static constexpr int MatSize = MatrixSize;
static constexpr int VecSize = BaseMatrix<T, MatrixSize, MaxVecSize>::VecSize;
public:
InvertibleMatrix() = default;
InvertibleMatrix(T const a[])
: BaseMatrix<T, MatrixSize, MaxVecSize>(a)
{}
InvertibleMatrix(T const a)
: BaseMatrix<T, MatrixSize, MaxVecSize>(a)
{}
private:
friend void chol_inv<matrix_t>(matrix_t &a, matrix_t &i);
friend void matrix_mul_m_m<matrix_t>(matrix_t &a, matrix_t &b, matrix_t &c);
friend void matrix_mul_mt_m<matrix_t>(matrix_t &a, matrix_t &b, matrix_t &c);
friend void matrix_mul_m_mt<matrix_t>(matrix_t &a, matrix_t &b, matrix_t &c);
};
template <typename M>
void chol_inv(M &a, M &inv) {
M r;
using pack_t = typename M::pack_t;
// Cholesky decomposition
for (int i=0; i<M::MatSize; i++) {
pack_t sum_row = bs::Zero<pack_t>();
for (int j=0; j<i; j++) {
pack_t cur_row(&r.array[j*M::VecSize]);
pack_t factor {r.array[j*M::VecSize+i]};
sum_row += cur_row*factor;
}
pack_t res_row(&a.array[i*M::VecSize]);
res_row -= sum_row;
pack_t diviser {1/std::sqrt(res_row[i])};
res_row *= diviser;
bs::aligned_store(res_row, &r.array[i*M::VecSize]);
}
M tmp;
static M Id(1);
// Inversion of R by backward substitution
for (int i=(M::MatSize-1); i>=0; i--) {
pack_t sum_row = bs::Zero<pack_t>();
for (int j=(M::MatSize-1); j > i; j--) {
pack_t cur_row(&tmp.array[j*M::VecSize]);
pack_t factor {r.array[i*M::VecSize+j]};
sum_row += cur_row*factor;
}
pack_t res_row(&Id.array[i*M::VecSize]);
res_row -= sum_row;
pack_t diviser {1/r.array[i*M::VecSize+i]};
res_row *= diviser;
bs::aligned_store(res_row, &tmp.array[i*M::VecSize]);
}
// Compute inverse of initial matrix
matrix_mul_mt_m(tmp, tmp, inv);
}
int main(int argc, char *argv[]) {
alignas(32) float a[SIZE*SIZE] {0};
alignas(32) float b[SIZE*SIZE];
alignas(32) float p[SIZE*SIZE];
for (int i=0; i<SIZE; i++) {
for (int j=0; j<SIZE; j++) {
if(j>=i) a[i*SIZE+j] = j-i+1;
}
}
for (int i=0; i<SIZE; i++) {
for (int j=0; j<SIZE; j++) {
float val;
if (i == 0 && j == 1) val = 1;
else if (i == 1 && j == 0) val = 1;
else if (i > 1 && i == j && i % 2) val = 1;
else if (i > 1 && i == j && !(i % 2)) val = 1;
else val = 0;
b[i*SIZE+j] = val;
}
}
if (SIZE % 2 == 0) {
b[SIZE*SIZE-1] = 1;
}
using M = InvertibleMatrix<float, SIZE, VECSIZE>;
M A(a), P(b), I, Tmp;
// std::cout << pa << std::endl;
// std::cout << pb << std::endl;
std::cout << P << std::endl;
matrix_mul_mt_m(A, A, Tmp); //Symetric positive definite matrix
A = Tmp;
std::cout << A << std::endl;
// for (unsigned int i=0; i<20000000; i++) {
for (unsigned int k = 0; k< 1000000; k++) {
A = Tmp;
matrix_mul_m_m(P, A, Tmp);
matrix_mul_mt_m(Tmp, P, A);
Tmp = A;
for (unsigned int i=0; i<20; i++) {
chol_inv(A, I);
chol_inv(I, A);
}
}
std::cout << A << std::endl;
return 0;
}
#include<iostream>
#include"fast_matrix.hpp"
template <typename T>
void print_matrix(T* matrix, int n, int m) {
for (int i=0; i<n; i++) {
for (int j=0; j<m; j++) {
std::cout << matrix[i*m+j] << (j == m-1 ? "" : ", ");
}
std::cout << std::endl;
}
}
void test_matrix_add() {
alignas(32) float input_a[3*3] = {
1, 2, 0,
2, 5, 1,
1, 1, 0.1,
};
alignas(32) float input_b[3*3] = {
2, 2, 4,
1, 0, 0.5,
1, 1, 2,
};
alignas(32) float truth[3*3] = {
3, 4, 4,
3, 5, 1.5,
2, 2, 2.1,
};
BaseMatrix<float, 3, 3, 8> A(input_a), B(input_b), C;
matrix_add(A, B, C);
float* output = C.dump_array();
std::cout << "Testing: Matrix Add" << std::endl;
std::cout << "\tOutput: " << std::endl;
print_matrix<float>(output, 3, 3);
std::cout << std::endl;
std::cout << "\tTruth: " << std::endl;
print_matrix<float>(truth, 3, 3);
std::cout << std::endl;
}
void test_matrix_sub() {
alignas(32) float input_a[3*3] = {
1, 2, 0,
2, 5, 1,
1, 1, 0.1,
};
alignas(32) float input_b[3*3] = {
2, 2, 4,
1, 0, 0.5,
1, 1, 2,
};
alignas(32) float truth[3*3] = {
-1, 0, -4,
1, 5, 0.5,
0, 0, -1.9,
};
BaseMatrix<float, 3, 3, 8> A(input_a), B(input_b), C;
matrix_sub(A, B, C);
float* output = C.dump_array();
std::cout << "Testing: Matrix sub" << std::endl;
std::cout << "\tOutput: " << std::endl;
print_matrix<float>(output, 3, 3);
std::cout << std::endl;
std::cout << "\tTruth: " << std::endl;
print_matrix<float>(truth, 3, 3);
std::cout << std::endl;
}
void test_matrix_mul_m_m() {
alignas(32) float input_a[3*3] = {
1, 2, 0,
2, 5, 1,
1, 1, 0.1,
};
alignas(32) float input_b[3*3] = {
2, 2, 4,
1, 0, 0.5,
1, 1, 2,
};
alignas(32) float truth[3*3] = {
4, 2, 5,
10, 5, 12.5,
3.1, 2.1, 4.7,
};
BaseMatrix<float, 3, 3, 8> A(input_a), B(input_b), C;
matrix_mul_m_m(A, B, C);
float* output = C.dump_array();
std::cout << "Testing: Matrix Mul M M" << std::endl;
std::cout << "\tOutput: " << std::endl;
print_matrix<float>(output, 3, 3);
std::cout << std::endl;
std::cout << "\tTruth: " << std::endl;
print_matrix<float>(truth, 3, 3);
std::cout << std::endl;
}
void test_identity() {
//IdentityMatrix<float, 5, 8> A;
BaseMatrix<float, 5, 5, 8> A(1);
std::cout << "Testing identity matrix" << std::endl;
std::cout << A;
}
void test_matrix_inv() {
alignas(32) float input[3*3] = {
1, 0, 0,
0, 2, 0,
0, 0, 4,
};
alignas(32) float truth[3*3] = {
1, 0, 0,
0, 0.5, 0,
0, 0, 0.25,
};
BaseMatrix<float, 3, 3, 4> A(input), I;
Inverse<float, 3, 3, 4>::inverse(A, I);
std::cout << "Testing Matrix inversion" << std::endl;
std::cout << "\tOutput: " << std::endl;
print_matrix<float>(I.dump_array(), 3, 3);
std::cout << "\tTruth: " << std::endl;
print_matrix<float>(truth, 3, 3);
std::cout << std::endl;
}
int main(int argc, char **argv) {
test_matrix_add();
test_matrix_sub();
test_matrix_mul_m_m();
test_identity();
test_matrix_inv();
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment