From f55824cdd36f58a5dd2af332e77178d0d6bb2be9 Mon Sep 17 00:00:00 2001
From: chn <chn@chn.moe>
Date: Mon, 9 Oct 2023 19:07:13 +0800
Subject: [PATCH] =?UTF-8?q?=E6=9A=82=E5=AD=98?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 include/ufo/unfold.hpp |  4 ++++
 src/fold.cpp           |  3 ++-
 src/unfold.cpp         | 41 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 47 insertions(+), 1 deletion(-)

diff --git a/include/ufo/unfold.hpp b/include/ufo/unfold.hpp
index a652e29..f6466ea 100644
--- a/include/ufo/unfold.hpp
+++ b/include/ufo/unfold.hpp
@@ -133,6 +133,10 @@ namespace ufo
       UnfoldSolver(std::string config_file);
       UnfoldSolver& operator()() override;
 
+      // 分解 SuperCellMultiplier
+      static std::tuple<Eigen::Matrix3i, Eigen::Vector3i> decompose_super_cell_multiplier
+        (Eigen::Matrix3i super_cell_multiplier);
+
       // 构建基
       // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。
       // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置。
diff --git a/src/fold.cpp b/src/fold.cpp
index bb6958a..706c6d1 100644
--- a/src/fold.cpp
+++ b/src/fold.cpp
@@ -40,7 +40,8 @@ namespace ufo
       for (auto& qpoint : Input_.Qpoints)
         Output_->Qpoints.push_back(fold
         (
-          qpoint, Input_.SuperCellMultiplier,
+          qpoint,
+          Input_.SuperCellMultiplier,
           Input_.SuperCellDeformation
         ));
     }
diff --git a/src/unfold.cpp b/src/unfold.cpp
index 3f9bbc5..fecb3f7 100644
--- a/src/unfold.cpp
+++ b/src/unfold.cpp
@@ -325,6 +325,47 @@ namespace ufo
     return *this;
   }
 
+  std::tuple<Eigen::Matrix3i, Eigen::Vector3i> UnfoldSolver::decompose_super_cell_multiplier
+    (Eigen::Matrix3i super_cell_multiplier)
+  {
+    struct
+    {
+      Eigen::Matrix3i SuperCellDeformation = Eigen::Matrix3i::Identity();
+      Eigen::Vector3i SuperCellSizeMultiplier;
+    } result;
+    // 生成交换两行的初等矩阵
+    auto elementary_exchange = [](int nrow1, int nrow2) -> Eigen::Matrix3i
+    {
+      Eigen::Matrix3i result = Eigen::Matrix3i::Identity();
+      result.row(nrow1).swap(result.row(nrow2));
+      return result;
+    };
+    // 生成将一行乘以一个整数加到另一行的初等矩阵及其逆矩阵
+    auto elementary_add = [](int nrow1, int nrow2, int multiplier)
+      -> std::pair<Eigen::Matrix3i, Eigen::Matrix3i>
+    {
+      Eigen::Matrix3i result1 = Eigen::Matrix3i::Identity(), result2 = Eigen::Matrix3i::Identity();
+      result1(nrow2, nrow1) = multiplier;
+      result2(nrow2, nrow1) = -multiplier;
+      return { result1, result2 };
+    };
+    for (int ncol = 0; ncol < 2; ncol++)
+      while (true)
+      {
+        // 将第 ncol 列从第 ncol 列开始绝对值最小的非零值调整到第 ncol 行
+        for (int nrow = ncol + 1; nrow < 3; nrow++)
+          if (std::abs(super_cell_multiplier(nrow, ncol)) < std::abs(super_cell_multiplier(ncol, ncol)))
+          {
+            super_cell_multiplier.row(ncol).swap(super_cell_multiplier.row(nrow));
+            result.SuperCellDeformation.col(ncol).swap(result.SuperCellDeformation.col(nrow));
+          }
+        // 如果第 ncol 列第 ncol + 1 行开始的元素为零, 则退出循环
+        if (super_cell_multiplier.block(ncol + 1, ncol, 2 - ncol, 1).isZero())
+          break;
+        // 尝试用
+      }
+  }
+
   UnfoldSolver::BasisType UnfoldSolver::construct_basis
   (
     const decltype(InputType::PrimativeCell)& primative_cell,