更新README并测试遍历和填充

README.cpp

@@ -4,7 +4,12 @@
 将每个main_xxx.cpp文件改为main.cpp，然后编译运行即可。
 
 目前测试了
-- BiCGStab.cpp
+- BiCGStab(一种迭代求解器)
+- Cholesky(一种分解求解器)
+- fillSparseMatrix(填充稀疏矩阵)
+- traverseSparseMatrix(遍历稀疏矩阵)
+- matrixFree(无矩阵迭代求解器)
+- SparseMatrix1(稀疏矩阵求解Laplace问题的案例, 采用SparseCholesky)
 
 ## data
 data文件夹下面的是matrix market格式的矩阵数据。

main-fillSparseMatrix.cpp

@@ -55,7 +55,7 @@ SparseMatrix<double> fillMatrix4() {
 
 int main()
 {
-    auto mat = fillMatrix2();
+    auto mat = fillMatrix4();
     std::cout<<"nonzero elements: "<<mat.nonZeros()<<std::endl;
 
 }

main-traverseSparseMatrix.cpp

@@ -1,33 +1,61 @@
-// main-traverseSparseMatrix.cpp
+// main-fillSparseMatrix.cpp
 #include <Eigen/Sparse>
+#include <Eigen/Dense> //支持dense matrix
+#include <unsupported/Eigen/SparseExtra> //支持从matrix market文件读取
 #include <iostream>
 using namespace Eigen;
 
-//遍历非零元素（如果要输出所有元素，直接cout<<A即可,但不建议）
-void traverseSparseMatrix(const SparseMatrix<double>& mat) {
-    for (int k = 0; k < mat.outerSize(); ++k)
-        for (SparseMatrix<double>::InnerIterator it(mat, k); it; ++it)
-        {
-            it.value();
-            it.row();   // row index
-            it.col();   // col index (here it is equal to k)
-            it.index(); // inner index, here it is equal to it.row()
-
-            static int step = 0;
-            std::cout << "step = " << step++ << std::endl;
-            std::cout << "it.value() = " << it.value() << std::endl;
-            // std::cout << "it.row() = " << it.row() << std::endl;
-            // std::cout << "it.col() = " << it.col() << std::endl;
-            // std::cout << "it.index() = " << it.index() << std::endl;
-        }
+const int N = 4, M = 4;
+
+//方法1：使用coeffRef()方法
+SparseMatrix<double> fillMatrix1() {
+  SparseMatrix<double> mat(N,M);
+  // reserve()方法用于预先分配空间，避免后续的内存分配
+  mat.reserve(VectorXi::Constant(M, 4)); // 4: estimated number of non-zero enties per column
+  mat.coeffRef(0,0) = 1;
+  mat.coeffRef(0,1) = 2.;
+  mat.coeffRef(1,1) = 3.;
+  mat.coeffRef(2,2) = 4.;
+  mat.coeffRef(2,3) = 5.;
+  mat.coeffRef(3,2) = 6.;
+  mat.coeffRef(3,3) = 7.;
+  mat.makeCompressed();
+  return mat;
+}
+
+//方法2：从dense matrix转换而来
+SparseMatrix<double> fillMatrix2() {
+    MatrixXd mat_dense = MatrixXd::Random(N,M);
+    Eigen::SparseMatrix<double> mat = mat_dense.sparseView(0.5, 1);
+    return mat;
+}
+
+//方法3：从matrix market文件读取
+SparseMatrix<double> fillMatrix3() {
+    SparseMatrix<double> mat;
+    Eigen::loadMarket(mat, "e05r0500.mtx");
+    return mat;
+}
+
+//方法4：使用setFromTriplets()方法(官方推荐)
+SparseMatrix<double> fillMatrix4() {
+  std::vector<Eigen::Triplet<double>> coefficients;
+
+  coefficients.push_back(Eigen::Triplet<double>(0,0,2));//row, col, value
+  coefficients.push_back(Eigen::Triplet<double>(0,1,6));
+  coefficients.push_back(Eigen::Triplet<double>(1,0,6));
+  coefficients.push_back(Eigen::Triplet<double>(1,1,1));
+  coefficients.push_back(Eigen::Triplet<double>(2,2,5));
+
+  Eigen::SparseMatrix<double> mat(N,M);
+  mat.setFromTriplets(coefficients.begin(), coefficients.end());
+  return mat;
 }
 
 
 int main()
 {
-    int n = 10, m = 5;
-    MatrixXd mat_dense = MatrixXd::Random(n,m);
-    Eigen::SparseMatrix<double> mat = mat_dense.sparseView(0.5, 1);
+    auto mat = fillMatrix4();
     std::cout<<"nonzero elements: "<<mat.nonZeros()<<std::endl;
-    TraverseSparseMatrix(mat);
+    std::cout<<mat<<std::endl;
 }

main.cpp

@@ -1,84 +1,4 @@
-//SparseMatrix1（去掉Qt的部分）
-
-#include <Eigen/Sparse>
-#include <vector>
-#include <iostream>
-#include <fstream>
-#include "utils.h"
-
-/* -------------------------------------------------------------------------- */
-/*                            构建稀疏矩阵，从Laplace方程构建                    */
-/* -------------------------------------------------------------------------- */
-#define _USE_MATH_DEFINES
-#include <math.h>
-
-typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
-typedef Eigen::Triplet<double> T;
-
-
-void insertCoefficient(int id, int i, int j, double w, std::vector<T>& coeffs,
-                       Eigen::VectorXd& b, const Eigen::VectorXd& boundary)
+int main()
 {
-  int n = int(boundary.size());
-  int id1 = i+j*n;
-
-        if(i==-1 || i==n) b(id) -= w * boundary(j); // constrained coefficient
-  else  if(j==-1 || j==n) b(id) -= w * boundary(i); // constrained coefficient
-  else  coeffs.push_back(T(id,id1,w));              // unknown coefficient
-}
-
-void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n)
-{
-  b.setZero();
-  Eigen::ArrayXd boundary = Eigen::ArrayXd::LinSpaced(n, 0,M_PI).sin().pow(2);
-  for(int j=0; j<n; ++j)
-  {
-    for(int i=0; i<n; ++i)
-    {
-      int id = i+j*n;
-      insertCoefficient(id, i-1,j, -1, coefficients, b, boundary);
-      insertCoefficient(id, i+1,j, -1, coefficients, b, boundary);
-      insertCoefficient(id, i,j-1, -1, coefficients, b, boundary);
-      insertCoefficient(id, i,j+1, -1, coefficients, b, boundary);
-      insertCoefficient(id, i,j,    4, coefficients, b, boundary);
-    }
-  }
-}
-
-
-/* -------------------------------------------------------------------------- */
-/*                                    求解矩阵                                 */
-/* -------------------------------------------------------------------------- */
-
-int main(int argc, char** argv)
-{
-  int n = 300;  // size of the image
-  int m = n*n;  // number of unknowns (=number of pixels)
-
-  // Assembly:
-  std::vector<T> coefficients;            // list of non-zeros coefficients
-  Eigen::VectorXd b(m);                   // the right hand side-vector resulting from the constraints
-  buildProblem(coefficients, b, n);
-
-  SpMat A(m,m);
-  A.setFromTriplets(coefficients.begin(), coefficients.end());
-
-  Profiler p;
-  p.start();
-  // Solving:
-  Eigen::SimplicialCholesky<SpMat> chol(A);  // performs a Cholesky factorization of A
-  Eigen::VectorXd x = chol.solve(b);         // use the factorization to solve for the given right hand side
-  p.end();
-
-  if (chol.info()==Eigen::Success)
-    std::cout << "Solve Success" << std::endl;
-  else
-    std::cout << "Fail" << std::endl;
-
-  std::ofstream myfile;
-  myfile.open ("SparseMatrix1.txt");
-  myfile << x;
-  myfile.close();
-
-  return 0;
+    return 0;
 }