Skip to content

Commit

Permalink
Merge pull request #3 from ypodlesov/MatrixPowersMV
Browse files Browse the repository at this point in the history 
Add MatrixPowersMV.
  • Loading branch information
@ypodlesov
ypodlesov authored May 18, 2024
2 parents be218ef + d0e9d0f commit 633b70b
Show file tree
Hide file tree
Showing 27 changed files with 270 additions and 159 deletions.
12 changes: 11 additions & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,16 @@ function(add_if_exists name)
endfunction()

set(BASIC_LA_PATH ${CMAKE_CURRENT_SOURCE_DIR}/basic_la)
set(GKLIB_PATH ${CMAKE_CURRENT_SOURCE_DIR}/GKlib)
set(METIS_INCLUDE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/METIS/include)
set(METIS_COPTIONS "${METIS_COPTIONS} -DIDXTYPEWIDTH=64")
set(METIS_COPTIONS "${METIS_COPTIONS} -DREALTYPEWIDTH=64")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${METIS_COPTIONS}")

add_if_exists(basic_la)
add_if_exists(qr_decomposition)
add_if_exists(conjugate_gradient)
add_if_exists(conjugate_gradient)

add_if_exists(GKlib)
add_if_exists(METIS)
add_if_exists(matrix_powers_mv)
12 changes: 6 additions & 6 deletions basic_la/common_container.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ struct CommonContainer {

CommonContainer() = default;

CommonContainer(size_t mem_size)
CommonContainer(int64_t mem_size)
: mem_size_{mem_size}
, data_{new T[mem_size_]}
{
Expand All @@ -21,7 +21,7 @@ struct CommonContainer {
CommonContainer(const CommonContainer& other) {
mem_size_ = other.mem_size_;
data_ = new T[mem_size_];
for (size_t i = 0; i < mem_size_; ++i) {
for (int64_t i = 0; i < mem_size_; ++i) {
data_[i] = other.data_[i];
}
}
Expand All @@ -44,7 +44,7 @@ struct CommonContainer {
}
mem_size_ = other.mem_size_;
data_ = new T[mem_size_];
for (size_t i = 0; i < mem_size_; ++i) {
for (int64_t i = 0; i < mem_size_; ++i) {
data_[i] = other.data_[i];
}
return *this;
Expand All @@ -61,15 +61,15 @@ struct CommonContainer {
template <typename std::enable_if<std::is_arithmetic<T>::value>::type* = nullptr>
void PlusAX(const CommonContainer& x, const T alpha = 1) {
assert(mem_size_ == x.mem_size_ && data_ && x.data_);
for (size_t i = 0; i < mem_size_; ++i) {
for (int64_t i = 0; i < mem_size_; ++i) {
data_[i] += alpha * x.data_[i];
}
}

template <typename std::enable_if<std::is_arithmetic<T>::value>::type* = nullptr>
void AXPlusBY(const CommonContainer& x, const T alpha, const CommonContainer& y, const T beta) {
assert(mem_size_ == x.mem_size_ && data_ && x.mem_size_ == y.mem_size_);
for (size_t i = 0; i < mem_size_; ++i) {
for (int64_t i = 0; i < mem_size_; ++i) {
data_[i] = x.data_[i] * alpha + y.data_[i] * beta;
}
}
Expand All @@ -87,6 +87,6 @@ struct CommonContainer {
}
}

size_t mem_size_;
int64_t mem_size_;
T* data_{nullptr};
};
16 changes: 8 additions & 8 deletions basic_la/common_matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ struct CommonMatrix: public CommonContainer<CommonMatrix<SpecMatrix, T, Hold>, T

CommonMatrix() = default;

CommonMatrix(size_t row_cnt, size_t col_cnt, size_t mem_size)
CommonMatrix(int64_t row_cnt, int64_t col_cnt, int64_t mem_size)
: Base(mem_size)
, row_cnt_{row_cnt}
, col_cnt_{col_cnt}
Expand Down Expand Up @@ -47,17 +47,17 @@ struct CommonMatrix: public CommonContainer<CommonMatrix<SpecMatrix, T, Hold>, T
return *this;
}

inline T Get(size_t row, size_t col) const {
inline T Get(int64_t row, int64_t col) const {
T result;
static_cast<const SpecMatrix*>(this)->GetImpl(row, col, result);
return result;
}

inline void Get(size_t row, size_t col, T& result) const {
inline void Get(int64_t row, int64_t col, T& result) const {
static_cast<const SpecMatrix*>(this)->GetImpl(row, col, result);
}

inline size_t GetMemSize() const {
inline int64_t GetMemSize() const {
return static_cast<const SpecMatrix*>(this)->GetMemSizeImpl();
}

Expand All @@ -66,8 +66,8 @@ struct CommonMatrix: public CommonContainer<CommonMatrix<SpecMatrix, T, Hold>, T
bool flag = true;
const SpecMatrix* mat = static_cast<const SpecMatrix*>(this);
assert(mat);
for (size_t i = 0; i < row_cnt_; ++i) {
for (size_t j = (is_upper ? 0 : i + 1); j < (is_upper ? i : col_cnt_); ++j) {
for (int64_t i = 0; i < row_cnt_; ++i) {
for (int64_t j = (is_upper ? 0 : i + 1); j < (is_upper ? i : col_cnt_); ++j) {
flag &= NHelpers::RoughEq(mat->GetImpl(i, j), 0.0);
}
}
Expand All @@ -78,6 +78,6 @@ struct CommonMatrix: public CommonContainer<CommonMatrix<SpecMatrix, T, Hold>, T
row_cnt_ = col_cnt_ = 0;
}

size_t row_cnt_{};
size_t col_cnt_{};
int64_t row_cnt_{};
int64_t col_cnt_{};
};
62 changes: 31 additions & 31 deletions basic_la/matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ struct Matrix: public CommonMatrix<Matrix<T>, T, Hold> {

Matrix() = default;

Matrix(size_t row_cnt, size_t col_cnt)
Matrix(int64_t row_cnt, int64_t col_cnt)
: Base(row_cnt, col_cnt, row_cnt * col_cnt)
{
}
Expand All @@ -38,41 +38,41 @@ struct Matrix: public CommonMatrix<Matrix<T>, T, Hold> {
return *this;
}

inline T& operator ()(size_t row, size_t col) {
inline T& operator ()(int64_t row, int64_t col) {
return data_[col * row_cnt_ + row];
}

inline size_t GetMemSizeImpl() const {
inline int64_t GetMemSizeImpl() const {
return row_cnt_ * col_cnt_;
}

inline static size_t GetMemSize(const size_t& row_cnt, const size_t& col_cnt) {
inline static int64_t GetMemSize(const int64_t& row_cnt, const int64_t& col_cnt) {
return row_cnt * col_cnt;
}

inline T GetImpl(size_t row, size_t col) const {
inline T GetImpl(int64_t row, int64_t col) const {
assert(data_);
T result;
GetImpl(row, col, result);
return result;
}

inline void GetImpl(size_t row, size_t col, T& result) const {
inline void GetImpl(int64_t row, int64_t col, T& result) const {
assert(data_);
result = data_[col * row_cnt_ + row];
}

static void MMMult(const Matrix& a, const Matrix& b, Matrix& c) {
assert(a.col_cnt_ == b.row_cnt_);
size_t m = a.row_cnt_;
size_t n = b.col_cnt_;
int64_t m = a.row_cnt_;
int64_t n = b.col_cnt_;
c = Matrix(m, n);
for (size_t j = 0; j < n; ++j) {
for (int64_t j = 0; j < n; ++j) {
T* c_j_col = &c.data_[j * c.row_cnt_];
for (size_t i = 0; i < m; ++i) {
for (int64_t i = 0; i < m; ++i) {
T* b_j_col = &b.data_[j * b.row_cnt_];
c_j_col[i] = 0;
for (size_t k = 0; k < a.col_cnt_; ++k) {
for (int64_t k = 0; k < a.col_cnt_; ++k) {
c_j_col[i] += a.data_[k * m + i] * b_j_col[k];
}
}
Expand All @@ -82,10 +82,10 @@ struct Matrix: public CommonMatrix<Matrix<T>, T, Hold> {
void VecMult(const Vector<T>& x, Vector<T>& result) const {
assert(result.data_ && result.mem_size_ == row_cnt_ && col_cnt_ == x.mem_size_);
NHelpers::Nullify(result.data_, result.mem_size_);
for (size_t j = 0; j < col_cnt_; ++j) {
for (int64_t j = 0; j < col_cnt_; ++j) {
T x_j = x.data_[j];
T* j_col = &data_[j * row_cnt_];
for (size_t i = 0; i < row_cnt_; ++i) {
for (int64_t i = 0; i < row_cnt_; ++i) {
result.data_[i] += j_col[i] * x_j;
}
}
Expand All @@ -94,14 +94,14 @@ struct Matrix: public CommonMatrix<Matrix<T>, T, Hold> {

template <typename T, typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
bool operator ==(const Matrix<T>& a, const Matrix<T>& b) {
size_t m = a.row_cnt_;
size_t n = a.col_cnt_;
int64_t m = a.row_cnt_;
int64_t n = a.col_cnt_;
if (m != b.row_cnt_ || n != b.col_cnt_) {
return false;
}
double diff = 0;
for (size_t j = 0; j < n; ++j) {
for (size_t i = 0; i < m; ++i) {
for (int64_t j = 0; j < n; ++j) {
for (int64_t i = 0; i < m; ++i) {
double tmp = std::fabs(a.data_[j * m + i] - b.data_[j * m + i]);
diff += tmp * tmp;
}
Expand All @@ -126,13 +126,13 @@ namespace NHelpers {
std::random_device random_device;
std::mt19937 random_engine(random_device());
std::uniform_real_distribution<T> distribution(0, a.row_cnt_ * a.col_cnt_);
for (size_t i = 0; i < a.mem_size_; ++i) {
for (int64_t i = 0; i < a.mem_size_; ++i) {
a.data_[i] = is_int ? int(distribution(random_engine)) : distribution(random_engine);
}
}

template <typename T>
void GenRandomMatrix(Matrix<T>& a, size_t row_cnt, size_t col_cnt, bool is_int = false) {
void GenRandomMatrix(Matrix<T>& a, int64_t row_cnt, int64_t col_cnt, bool is_int = false) {
a = Matrix<T>(row_cnt, col_cnt);
FillRandomMatrix(a, is_int);
}
Expand All @@ -142,37 +142,37 @@ namespace NHelpers {
std::random_device random_device;
std::mt19937 random_engine(random_device());
std::uniform_real_distribution<T> distribution(0, a.row_cnt_ * a.col_cnt_);
for (size_t j = 0; j < a.col_cnt_; ++j) {
for (int64_t j = 0; j < a.col_cnt_; ++j) {
if (type == TriangularType::Upper) {
for (size_t i = 0; i <= j; ++i) {
for (int64_t i = 0; i <= j; ++i) {
a.data_[j * a.row_cnt_ + i] = is_int ? int(distribution(random_engine)) : distribution(random_engine);
}
} else {
for (size_t i = j; i < a.row_cnt_; ++i) {
for (int64_t i = j; i < a.row_cnt_; ++i) {
a.data_[j * a.row_cnt_ + i] = is_int ? int(distribution(random_engine)) : distribution(random_engine);
}
}
}
}

template <typename T>
void GenRandomTriangularMatrix(Matrix<T>& a, size_t row_cnt, size_t col_cnt, TriangularType type, bool is_int = false) {
void GenRandomTriangularMatrix(Matrix<T>& a, int64_t row_cnt, int64_t col_cnt, TriangularType type, bool is_int = false) {
a = Matrix<T>(row_cnt, col_cnt);
FillRandomTriangularMatrix(a, type, is_int);
}

template <typename T>
void FillRandomSymmetricMatrix(Matrix<T>& a, bool is_int = false) {
FillRandomTriangularMatrix(a, TriangularType::Upper, is_int);
for (size_t j = 0; j < a.col_cnt_; ++j) {
for (size_t i = j + 1; i < a.row_cnt_; ++i) {
for (int64_t j = 0; j < a.col_cnt_; ++j) {
for (int64_t i = j + 1; i < a.row_cnt_; ++i) {
a.data_[j * a.row_cnt_ + i] = a.data_[i * a.row_cnt_ + j];
}
}
}

template <typename T>
void GenRandomSymmetricMatrix(Matrix<T>& a, size_t row_cnt, size_t col_cnt, bool is_int = false) {
void GenRandomSymmetricMatrix(Matrix<T>& a, int64_t row_cnt, int64_t col_cnt, bool is_int = false) {
a = Matrix<T>(row_cnt, col_cnt);
FillRandomSymmetricMatrix(a, is_int);
}
Expand All @@ -183,15 +183,15 @@ namespace NHelpers {
std::random_device random_device;
std::mt19937 random_engine(random_device());
std::uniform_real_distribution<T> distribution(T{}, a.row_cnt_ * a.col_cnt_);
for (size_t i = 0; i < a.row_cnt_; ++i) {
for (int64_t i = 0; i < a.row_cnt_; ++i) {
a.data_[i * a.row_cnt_ + i] = definiteness == Definiteness::Positive
? std::fabs(a.data_[i * a.row_cnt_ + i]) + 1
: -std::fabs(a.data_[i * a.row_cnt_ + i]) - 1;
}
}

template <typename T>
void GenRandomDefiniteMatrix(Matrix<T>& a, size_t row_cnt, size_t col_cnt, Definiteness type, bool is_int = false) {
void GenRandomDefiniteMatrix(Matrix<T>& a, int64_t row_cnt, int64_t col_cnt, Definiteness type, bool is_int = false) {
a = Matrix<T>(row_cnt, col_cnt);
FillRandomDefiniteMatrix(a, type, is_int);
}
Expand All @@ -200,11 +200,11 @@ namespace NHelpers {

template <typename T>
std::istream & operator >>(std::istream &in, Matrix<T> &matrix) {
size_t sz1, sz2;
int64_t sz1, sz2;
in >> sz1 >> sz2;
matrix = Matrix<T>(sz1, sz2);
for (size_t i = 0; i < sz1; ++i) {
for (size_t j = 0; j < sz2; ++j) {
for (int64_t i = 0; i < sz1; ++i) {
for (int64_t j = 0; j < sz2; ++j) {
in >> matrix.data_[j * sz1 + i];
}
}
Expand Down
Loading

0 comments on commit 633b70b

Please sign in to comment.