first commit

.gitignore

@@ -0,0 +1,2 @@
+test/*
+!test/*.cpp

README

@@ -0,0 +1,30 @@
+matris - a simple matrix class with some MATLAB-like functions
+
+base functions:
+
+include/lu.hpp   LU分解
+include/inv.hpp  求逆
+include/tr.hpp   转置
+include/det.hpp  求行列式的值
+include/rref.hpp 将矩阵化为行阶梯型
+include/rank.hpp 矩阵的秩
+include/norm.hpp 范数
+include/cond.hpp 条件数
+include/diag.hpp 取矩阵对角元素的值
+
+extend functions:
+
+include/guassian_elimination.hpp
+用高斯消元部分选主元法求解线性方程组
+include/qr.hpp
+基于豪斯霍尔德变换的QR分解
+include/givens_rotation.hpp
+基于吉文斯变换的QR分解
+include/power_method.hpp
+幂法求最大特征值及其特征向量
+include/qr_method.hpp
+QR算法求一般矩阵的所有特征值
+include/least_squares.hpp
+基于QR分解的最小二乘问题解法
+include/gram_schmidt_modified.hpp
+基于改进的格拉姆-斯密特正交化过程的QR分解

include/cond.hpp

@@ -0,0 +1,13 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef COND_HPP
+#define COND_HPP
+
+template <typename T>
+T cond( const Matrix<T> &A, string p="2" ) {
+    T cd = norm(A,p)*norm(inv(A),p);
+    return cd;
+}
+
+#endif

include/det.hpp

@@ -0,0 +1,36 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef DET_HPP
+#define DET_HPP
+
+template <typename T>
+T det( const Matrix<T> &A ) {
+    // argument must be a square matrix
+    assert(A.rows() == A.cols());        
+
+    int n = A.rows();
+    vector< Matrix<T> > V;
+
+    V = lu(A);
+
+    if(V.empty()) return 0;
+
+    Matrix<T> U = V.at(1);
+    Matrix<T> F = V.at(3);
+
+    int count = F[0][0];
+
+    T value = U[0][0];
+
+    for (int i = 1; i < n; ++i) {
+        value *= U[i][i];
+    }
+
+    if(count%2) value = -value;
+    if(!value) value = 0;
+
+    return value;
+}
+
+#endif

include/diag.hpp

@@ -0,0 +1,18 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef DIAG_HPP
+#define DIAG_HPP
+
+template <typename T>
+Matrix<T> diag( const Matrix<T> &A ) {
+    int n = A.rows() < A.cols() ? A.rows() : A.cols();
+    Matrix<T> d(n,1);
+    for (int i = 0; i < n; ++i) {
+        d[i][0] = A[i][i];
+    }
+
+    return d;
+}
+
+#endif

include/givens_rotation.hpp

@@ -0,0 +1,44 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef GIVENS_ROTATION_HPP
+#define GIVENS_ROTATION_HPP
+
+template <typename T>
+vector< Matrix<T> > givens_rotation( const Matrix<T> &A ) {
+    // argument must be a square matrix
+    assert(A.rows() == A.cols());
+
+    int n = A.rows();
+
+    Matrix<T> R(A);
+    Matrix<T> Q(n,n);
+    Q = 1;
+
+    for (int i = 0; i < n; ++i) {
+        for (int j = i+1; j < n; ++j) {
+            Matrix<T> G(n,n);
+            G = 1;
+
+            T r = sqrt( R[i][i]*R[i][i] + R[j][i]*R[j][i] );
+            T c = R[i][i]/r;
+            T s = R[j][i]/r;
+
+            G[i][i] = c;
+            G[i][j] = s;
+            G[j][i] = s;
+            G[j][j] = -c;
+
+            R = G*R;
+            Q = G*Q;
+        }
+    }
+
+    vector< Matrix<T> > V;
+    V.push_back(Q);
+    V.push_back(R);
+
+    return V;
+}
+
+#endif

include/gram_schmidt_modified.hpp

@@ -0,0 +1,48 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef GRAM_SCHMIDT_MODIFIED_HPP
+#define GRAM_SCHMIDT_MODIFIED_HPP
+
+template <typename T>
+vector< Matrix<T> > gram_schmidt_modified( const Matrix<T> &A ) {
+    int n = A.cols();
+    int m = A.rows();
+
+    vector< Matrix<T> > a,q;
+
+    Matrix<T> Q(m,n);
+    Matrix<T> R(n,n);
+
+    for (int i = 0; i < n; ++i) {
+        Matrix<T> tmp = A.cols(i);
+        a.push_back(tmp);
+        q.push_back(tmp);
+    }
+
+    for (int k = 0; k < n; ++k) {
+        R[k][k] = norm(a.at(k));
+        q.at(k) = (1/R[k][k])*a.at(k);
+
+        for (int j = k+1; j < n; ++j) {
+            Matrix<T> tmp = tr(q.at(k))*a.at(j);
+            R[k][j] = tmp[0][0];
+            a.at(j) = a.at(j) - R[k][j]*q.at(k);
+        }
+    }
+
+    for (int i = 0; i < n; ++i) {
+        Matrix<T> tmp = q.at(i);
+        for (int j = 0; j < m; ++j) {
+            Q[j][i] = tmp[j][0];
+        }
+    }
+
+    vector< Matrix<T> > V;
+    V.push_back(Q);
+    V.push_back(R);
+
+    return V;
+}
+
+#endif

include/guassian_elimination.hpp

@@ -0,0 +1,68 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef GUASSIAN_ELIMINATION_HPP
+#define GUASSIAN_ELIMINATION_HPP
+
+template <typename T>
+Matrix<T> guassian_elimination(const Matrix<T> &A, const Matrix<T> &b) {
+    int n = A.rows();
+    double eps = 1.0e-10;
+
+    Matrix<T> B(A);
+    Matrix<T> a(b);
+
+    for (int i = 0; i < n-1; ++i) {
+        T maxi = 0, tmpi;
+        int nth;
+
+        for (int j = i; j< n; ++j) {
+            if( (tmpi = abs(B[j][i])) > maxi ) {
+                maxi = tmpi;
+                nth = j;
+            }
+        }
+
+        if(maxi < eps) {
+            return Matrix<T>(n,1);
+        }
+
+        if(nth > i) {
+            swap(rows(B)[nth],rows(B)[i]);
+            swap(rows(a)[nth],rows(a)[i]);
+        }
+
+        for (int j = i+1; j < n; ++j) {
+            T tmp = B[j][i]/B[i][i];
+
+            for (int k = i+1; k < n; ++k) {
+                B[j][k] -= tmp*B[i][k];
+            }
+
+            for (int k = 0; k < (int)a.cols(); ++k) {
+                a[j][k] -= tmp*a[i][k];
+            }
+
+        }
+    }
+
+
+    for (int i = n-1; i >= 0; --i) {
+        for ( int j = i+1; j < n; ++j) {
+            for (int k = 0; k < (int)a.cols(); ++k) {
+                a[i][k] -= B[i][j]*a[j][k];
+            }
+        }
+
+        for (int j = 0; j < (int)a.cols(); ++j) {
+            a[i][j] /= B[i][i];
+            if (abs(a[i][j]) < eps) {
+                a[i][j] = 0;
+            }
+        }
+    }
+
+    return a;
+}
+
+#endif

include/inv.hpp

@@ -0,0 +1,84 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef INV_HPP
+#define INV_HPP
+
+template <typename T>
+Matrix<T> inv( const Matrix<T> &A ) {
+    // argument must be a square matrix
+    assert(A.rows() == A.cols());
+
+    int n = A.rows();
+    double eps = 1.0e-10;
+
+    Matrix<T> B(A);
+    Matrix<T> I(n,n);
+    I = 1;
+
+    for (int i = 0; i < n-1; ++i) {
+        T maxi = 0, tmpi;
+        int nth;
+
+        for (int j = i; j< n; ++j) {
+            if( (tmpi = abs(B[j][i])) > maxi ) {
+                maxi = tmpi;
+                nth = j;
+            }
+        }
+
+        if(maxi < eps) {
+            return Matrix<T>(n,n);
+        }
+
+        if(nth > i) {
+            swap(rows(B)[nth],rows(B)[i]);
+            swap(rows(I)[nth],rows(I)[i]);
+        }
+
+        for (int j = i+1; j < n; ++j) {
+            T tmp = B[j][i]/B[i][i];
+            B[j][i] = 0;
+            for (int k = i+1; k < n; ++k) {
+                B[j][k] -= tmp*B[i][k];
+                if(abs(B[j][k]) < eps) B[j][k] = 0;
+            }
+
+            for (int k = 0; k < n; ++k) {
+                I[j][k] -= tmp*I[i][k];
+                if(abs(I[j][k]) < eps) I[j][k] = 0;
+            }
+        }
+    }
+
+
+    for (int i = 1; i < n; ++i) {
+        for (int j = i-1; j >= 0; --j) {
+            T tmp = B[j][i]/B[i][i];
+            B[j][i] = 0;
+            for (int k = i+1; k < n; ++k) {
+                B[j][k] -= tmp*B[i][k];
+                if(abs(B[j][k]) < eps) B[j][k] = 0;
+            }
+
+            for (int k = 0; k < n; ++k) {
+                I[j][k] -= tmp*I[i][k];
+                if(abs(I[j][k]) < eps) I[j][k] = 0;
+            }
+        }
+    }
+
+
+    for (int i = 0; i < n; ++i) {
+        T tmp = B[i][i];
+        B[i][i] = 1;
+
+        for (int j = 0; j < n; ++j) {
+            I[i][j] /= tmp;
+        }
+    }
+
+    return I;
+}
+
+#endif

include/least_squares.hpp

@@ -0,0 +1,41 @@
+// Author:   timebug
+// License:  BSD
+
+#ifndef LEAST_SQUARES_HPP
+#define LEAST_SQUARES_HPP
+
+template <typename T>
+Matrix<T> least_squares( const Matrix<T> &A, const Matrix<T> &b ) {
+    int m = A.rows();
+    int n = A.cols();
+
+    vector< Matrix<T> > V;
+
+    V = qr(A);
+
+    Matrix<T> Q = V.at(0); // mxm
+    Matrix<T> R = V.at(1); // mxn
+
+    Matrix<T> Q1(m,n);
+    Matrix<T> R1(n,m);
+
+    for (int i = 0; i < m; ++i) {
+        for (int j = 0; j < n; ++j) {
+            Q1[i][j] = Q[i][j];
+        }
+    }
+
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < m; ++j) {
+            R1[i][j] = R[i][j];
+        }
+    }
+
+    Matrix<T> c1 = tr(Q1)*b; // nx1 = (nxm)*(mx1)
+
+    Matrix<T> x = guassian_elimination(R1,c1);
+
+    return x;
+}
+
+#endif