Initial commit.

.gitignore

@@ -0,0 +1,39 @@
+### C++ ###
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+
+build*
+*json*
+# *fortran*
+
+.cache

CMakeLists.txt

@@ -0,0 +1,33 @@
+cmake_minimum_required(VERSION 3.14)
+project(cmc_ctm_matrices)
+
+set(CMAKE_CXX_STANDARD             20)
+set(CMAKE_MODULE_PATH              "${CMAKE_SOURCE_DIR}/tools/cmake")
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}")
+set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}")
+set(CMAKE_EXPORT_COMPILE_COMMANDS  ON)
+
+find_package(Catch REQUIRED)
+
+set(Clang_DIR "/usr/lib/llvm-16/lib/cmake/clang")
+set(LLVM_DIR "/usr/lib/llvm-16/lib/cmake/llvm")
+find_package(Clang CONFIG)
+
+include_directories(tools/commons)
+include(tools/cmake/TestSolution.cmake)
+include(tools/cmake/BuildFlags.cmake)
+
+if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/private)
+  add_subdirectory(private)
+endif()
+
+function(add_if_exists name)
+  if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${name})
+    add_subdirectory(${name})
+  endif()
+endfunction()
+
+set(BASIC_LA_PATH ${CMAKE_CURRENT_SOURCE_DIR}/basic_la)
+add_if_exists(basic_la)
+add_if_exists(qr_decomposition)
+add_if_exists(conjugate_gradient)

README.md

@@ -0,0 +1,50 @@
+# Basic matrix algorithms
+
+This repository is devoted to study course of computational linear algebra in  **MSU-CMC-CTM** (Moscow State University, Faculty of Computational Mathematics and Cybernetics, Department of Computational Technologies and Modeling).
+
+Here you can also find use-examples of some instruments like **conan package manager**, **google-test**, **google-benchmark** and maybe something else.
+
+## Project file tree
+
+`git ls-tree -r --name-only HEAD | tree --fromfile`
+
+```bash
+cmc-ctm-matrix-algos
+├── basic_la
+│   ├── CMakeLists.txt
+│   ├── helpers.h
+│   ├── matrix.cpp
+│   ├── matrix.h
+│   ├── vector.cpp
+│   └── vector.h
+├── CMakeLists.txt
+├── conjugate_gradient
+│   ├── bench_results.txt
+│   ├── CMakeLists.txt
+│   ├── conjugate_gradient.cpp
+│   ├── conjugate_gradient.h
+│   ├── README.md
+│   ├── run.cpp
+│   └── test.cpp
+├── .gitignore
+├── qr_decomposition
+│   ├── CMakeLists.txt
+│   ├── qr_decomposition.cpp
+│   ├── qr_decomposition.h
+│   ├── run.cpp
+│   └── test.cpp
+├── README.md
+└── tools
+    ├── catch_main.cpp
+    ├── cmake
+    │   ├── BuildFlags.cmake
+    │   ├── FindCatch.cmake
+    │   ├── FindFFTW.cmake
+    │   ├── Protobuf.cmake
+    │   └── TestSolution.cmake
+    └── commons
+        ├── run_channel.h
+        ├── runner.h
+        ├── test_channel.h
+        └── util.h
+```

basic_la/CMakeLists.txt

@@ -0,0 +1 @@
+add_library(BASIC_LA STATIC matrix.cpp common_matrix.cpp vector.cpp)

basic_la/common_container.cpp

@@ -0,0 +1 @@
+#include "common_container.h"

basic_la/common_container.h

@@ -0,0 +1,98 @@
+#pragma once
+#include "helpers.h"
+#include <cassert>
+#include <cstddef>
+#include <cstdlib>
+#include <ctime>
+#include <type_traits>
+#include <utility>
+
+template <typename ContainerType, typename T, bool Hold = true>
+struct CommonContainer {
+
+    CommonContainer() = default;
+
+    CommonContainer(size_t mem_size)
+        : mem_size_{mem_size}
+        , data_{new T[mem_size_]}
+    {
+    }
+
+    CommonContainer(const CommonContainer& other) {
+        mem_size_ = other.mem_size_;
+        data_ = new T[mem_size_];
+        for (size_t i = 0; i < mem_size_; ++i) {
+            data_[i] = other.data_[i];
+        }
+    }
+
+    CommonContainer(CommonContainer&& other) {
+        if (data_ && Hold) {
+            delete[] data_;
+            data_ = nullptr;
+        }
+        std::swap(data_, other.data_);
+        other.data_ = nullptr;
+        std::swap(mem_size_, other.mem_size_);
+        other.mem_size_ = 0;
+    }
+
+    CommonContainer& operator =(const CommonContainer& other) {
+        if (data_ && Hold) {
+            delete[] data_;
+            data_ = nullptr;
+        }
+        mem_size_ = other.mem_size_;
+        data_ = new T[mem_size_];
+        for (size_t i = 0; i < mem_size_; ++i) {
+            data_[i] = other.data_[i];
+        }
+        return *this;
+    }
+
+    CommonContainer& operator =(CommonContainer&& other) {
+        std::swap(data_, other.data_);
+        other.data_ = nullptr;
+        std::swap(mem_size_, other.mem_size_);
+        other.mem_size_ = 0;
+        return *this;
+    }
+
+    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_);
+        for (size_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 = 1) {
+        assert(x.mem_size_ == y.mem_size_);
+        if (data_ && Hold) {
+            delete[] data_;
+            data_ = nullptr;
+        }
+        mem_size_ = x.mem_size_;
+        data_ = new T[mem_size_];
+        for (size_t i = 0; i < mem_size_; ++i) {
+            data_[i] = x.data_[i] * alpha + y.data_[i] * beta;
+        }
+    }
+
+    template <typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
+    inline T Norm2() const noexcept {
+        T result = NHelpers::InnerProd(data_, data_, mem_size_);
+        return std::sqrt(result);
+    }
+
+    ~CommonContainer() {
+        if (data_ && Hold) {
+            delete[] data_;
+            data_ = nullptr;
+        }
+    }
+
+    size_t mem_size_;
+    T* data_{nullptr};
+};

basic_la/common_matrix.cpp

@@ -0,0 +1 @@
+#include "common_matrix.h"

basic_la/common_matrix.h

@@ -0,0 +1,90 @@
+#pragma once
+#include "common_container.h"
+#include "helpers.h"
+#include "vector.h"
+#include <utility>
+
+template <typename SpecMatrix, typename T, bool Hold = true>
+struct CommonMatrix: public CommonContainer<CommonMatrix<SpecMatrix, T, Hold>, T, Hold> {
+    using Base = CommonContainer<CommonMatrix<SpecMatrix, T, Hold>, T, Hold>;
+    using Base::data_;
+
+    CommonMatrix() = default;
+
+    CommonMatrix(size_t row_cnt, size_t col_cnt, size_t mem_size)
+        : Base(mem_size)
+        , row_cnt_{row_cnt}
+        , col_cnt_{col_cnt}
+    {
+    }
+
+    CommonMatrix(CommonMatrix&& other)
+        : Base(std::move(other))
+        , row_cnt_{other.row_cnt_}
+        , col_cnt_{other.col_cnt_}
+    {
+    }
+
+    CommonMatrix& operator =(const CommonMatrix& other) {
+        Base::operator =(other);
+        row_cnt_ = other.row_cnt_;
+        col_cnt_ = other.col_cnt_;
+        return *this;
+    }
+
+    CommonMatrix& operator =(CommonMatrix&& other) {
+        Base::operator =(std::move(other));
+        std::swap(row_cnt_, other.row_cnt_);
+        std::swap(col_cnt_, other.col_cnt_);
+        return *this;
+    }
+
+    inline T Get(size_t row, size_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 {
+        static_cast<const SpecMatrix*>(this)->GetImpl(row, col, result);
+    }
+
+    inline size_t GetMemSize() const {
+        return static_cast<const SpecMatrix*>(this)->GetMemSizeImpl();
+    }
+
+    bool IsTriangular(NHelpers::ETriangularType type) const {
+        bool is_upper = type == NHelpers::ETriangularType::Upper;
+        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) {
+                flag &= NHelpers::RoughEq(mat->GetImpl(i, j), 0.0);
+            }
+        }
+        return flag;
+    }
+
+    void VecMult(const Vector<T>& x, Vector<T>& result) const {
+        if (!result.data_ || result.mem_size_ != row_cnt_) {
+            result = Vector<T>(row_cnt_);
+        }
+        NHelpers::Nullify(result.data_, result.mem_size_);
+        auto* original = static_cast<const SpecMatrix*>(this);
+        assert(original);
+        for (size_t j = 0; j < col_cnt_; ++j) {
+            T x_j = x.data_[j];
+            for (size_t i = 0; i < row_cnt_; ++i) {
+                result.data_[i] += original->GetImpl(i, j) * x_j;
+            }
+        }
+    }
+
+    ~CommonMatrix() {
+        row_cnt_ = col_cnt_ = 0;
+    }
+
+    size_t row_cnt_{};
+    size_t col_cnt_{};
+};

basic_la/helpers.h

@@ -0,0 +1,45 @@
+#pragma once
+#include <cmath>
+#include <ctime>
+
+namespace NHelpers {
+
+    constexpr double EPS = 1e-9;
+
+    template <class T1, class T2>
+    bool RoughEq(T1 lhs, T2 rhs, T2 epsilon = EPS) { 
+        return std::fabs(lhs - rhs) < epsilon;
+    }
+
+    template <class T1, class T2>
+    bool RoughLT(T1 lhs, T2 rhs, T2 epsilon = EPS) {
+        return rhs - lhs >= epsilon;
+    }
+
+    template <class T1, class T2>
+    bool RoughtLTE(T1 lhs, T2 rhs, T2 epsilon = EPS) {
+        return rhs - lhs > -epsilon;
+    }
+
+    template <typename T>
+    void Nullify(T* data_, size_t n) {
+        for (size_t i = 0; i < n; ++i) {
+            data_[i] = {};
+        }
+    }
+
+    template <typename T>
+    T InnerProd(T* data_1, T* data_2, size_t n) {
+        T result{};
+        for (size_t i = 0; i < n; ++i) {
+            result += data_1[i] * data_2[i];
+        }
+        return result;
+    }
+
+    enum class ETriangularType {
+        Upper,
+        Lower
+    };
+
+} // namespace NHelpers

basic_la/matrix.cpp

@@ -0,0 +1 @@
+#include "matrix.h"