一个头文件告诉你 C++ 中矩阵、向量及其相关运算(矩阵乘向量等)的定义

如何优雅地写好科学计算的 C 代码

动机

我有个同学叫xxx。那天和xx一起去上课,他坐在我旁边,抱着 Y7000P,在 Qt 里敲着计算流体力学的 C 代码,不一会儿,洋洋洒洒码出近百行的代码。定睛一看,密密麻麻地,三层、四层、五层,全是 for 循环,除了 for 循环,似乎没有别的东西了。一个搞高性能计算的精神小伙,对于写 for 循环,乐此不疲。

我们做科学计算的,免不了在程序中大批次大规模大概率地和矩阵和向量打交道,向量内积,外积,加减乘除,矩阵乘向量,矩阵乘矩阵,求范数等等。比起 MATLAB,C 如果不调线性代数的库的话,向量和矩阵的操作要写不少的 for 循环。既然如此,我们可以提前定义好矩阵和向量的数据结构,重载四则运算和下标索引,模仿 MATLAB 定义好 A*x 等操作 ,提高可复用性,让代码看起来更简洁优雅,让写出来的程序带的 bug 更少。

说干就干。把双手从 for 循环中解放出来。

实现的功能

  • 任意类型,任意长度的数组
  • double 类型的向量
  • double 类型的矩阵
  • 各个类型的构造和初始化
  • 迭代器的实现
  • 向量[] 索引取值
  • 返回size
  • 比较大小
  • 返回变量地址
  • += -= *= /=+ - * /的重载
  • 范数及其平方
  • 质心、凸包、开方、绝对值
  • 内积、叉积、外积
  • 标准输出
  • 向量标准归一化
  • 标准基、对角矩阵、矩阵转置
  • 矩阵()索引取值
  • 矩阵取某一列
  • 返回行数、列数
  • 计算 A T A A^TA ATA A T v A^Tv ATv x T A y x^TAy xTAy A x Ax Ax等。
  • 矩阵的迹
  • ……

矩阵、向量定义的 h 文件

看代码,不解释。

#pragma once
/// \file container.h
/// \brief Define some useful vector or matrix container for scientific computation
/// \author LSEC: Song Lu, lusong@lsec.cc.ac.cn


#include <vector>
#include <list>
#include <cmath>
#include <iostream>
#include <valarray>
#include <stdexcept>

enum InitStateT { Uninitialized, Initialized };
typedef unsigned int Uint;
template <class T, Uint _Size>
class SArrayCL;
template <Uint _Size>
class SVectorCL;
template <Uint _Rows, Uint _Cols>
class SMatrixCL;
/// Stores 2D coordinates
typedef SVectorCL<2> Point2DCL;
/// Stores 3D coordinates
typedef SVectorCL<3> Point3DCL;
/// Stores 4D coordinates
typedef SVectorCL<4> Point4DCL;
/// Stores barycentric coordinates of a space-time simplex (pentatope)
typedef SVectorCL<5> STBaryCoordCL;
template <class T, Uint _Size>
inline bool
operator==(const SArrayCL<T, _Size>&, const SArrayCL<T, _Size>&);
template <class T, Uint _Size>
inline bool
operator<(const SArrayCL<T, _Size>&, const SArrayCL<T, _Size>&);


//**************************************************************************
// Class:    SArrayCL                                                      *
// Purpose:  an array that remembers its size                              *
// Remarks:  All functions are inline, should be as fast as a "bare" array *
//**************************************************************************
template <class T, Uint _Size>
class SArrayCL
{
private:
    T Array[_Size];

public:
    typedef       T* iterator;
    typedef const T* const_iterator;
    typedef       T& reference;
    typedef const T& const_reference;
    typedef       T  value_type;

    //    SArrayCL() {}
    explicit           SArrayCL(T val = T()) { std::fill_n(Array + 0, _Size, val); }
    /*uninitialized memory; mainly for faster SVectorCL-math*/
    explicit SArrayCL(InitStateT) {}
    template<class In> explicit SArrayCL(In start) { std::copy(start, start + _Size, Array + 0); }
    template<class In> SArrayCL(In start, In end) { std::copy(start, end, Array + 0); }
    // Default copy-ctor, assignment operator, dtor

    iterator       begin() { return static_cast<T*>(Array); }
    const_iterator begin()       const { return static_cast<const T*>(Array); }
    iterator       end() { return Array + _Size; }
    const_iterator end()       const { return Array + _Size; }
    reference      operator[](Uint i) {
        if (!(i < _Size))
            throw std::exception("SArrayCL::operator[]: wrong index");
        //Assert(i < _Size, DROPSErrCL("SArrayCL::operator[]: wrong index"), DebugContainerC);
        return Array[i];
    }
    value_type     operator[](Uint i) const {
        if (!(i < _Size))
            throw std::exception("SArrayCL::operator[]: wrong index");
        //Assert(i < _Size, DROPSErrCL("SArrayCL::operator[]: wrong index"), DebugContainerC);
        return Array[i];
    }
    Uint           size()       const { return _Size; }

    friend bool operator==<>(const SArrayCL&, const SArrayCL&); // Component-wise equality
    friend bool operator< <>(const SArrayCL&, const SArrayCL&); // lexicographic ordering
};

template <class T, Uint _Size>
inline bool
operator==(const SArrayCL<T, _Size>& a0, const SArrayCL<T, _Size>& a1)
{
    for (Uint i = 0; i < _Size; ++i)
        if (a0[i] != a1[i]) return false;
    return true;
}



template <class T, Uint _Size>
inline bool
operator<(const SArrayCL<T, _Size>& a0, const SArrayCL<T, _Size>& a1)
{
    for (Uint i = 0; i < _Size; ++i)
        if (a0[i] < a1[i]) return true;
        else if (a0[i] > a1[i]) return false;
    return false;
}

template <class T, Uint _Size>
inline const T*
Addr(const SArrayCL<T, _Size>& a)
{
    return a.begin();
}

template <class T, Uint _Size>
inline T*
Addr(SArrayCL<T, _Size>& a)
{
    return a.begin();
}


//**************************************************************************
// Class:    SVectorCL                                                     *
// Purpose:  primitive vector class with templates for short vectors       *
// Remarks:  Many optimizations are possible.                              *
//**************************************************************************

template <Uint _Size>
class SVectorCL : public SArrayCL<double, _Size>
{
public:
    typedef SArrayCL<double, _Size> base_type;

    SVectorCL() {}
    explicit           SVectorCL(InitStateT i) : base_type(i) {}
    explicit           SVectorCL(double val) : base_type(val) {}
    // constructors for Point3D/Point4D/BaryCoordCL 
    // (constructor throws Exception if number of arg. does not coincide with _Size)
    explicit SVectorCL(double arg1, double arg2, double arg3);
    explicit SVectorCL(double arg1, double arg2, double arg3, double arg4);
    template<Uint _Size2>
    explicit SVectorCL(SVectorCL<_Size2> arg1, double arg2);

    template<class In> explicit SVectorCL(In start) : base_type(start) {}
    template<class In> SVectorCL(In start, In end) : base_type(start, end) {}

    SVectorCL& operator+=(const SVectorCL&);
    SVectorCL& operator-=(const SVectorCL&);
    SVectorCL& operator*=(double);
    SVectorCL& operator/=(double);

    // komponentenweise Operatoren
    SVectorCL& operator*=(const SVectorCL&);
    SVectorCL& operator/=(const SVectorCL&);

    double norm_sq() const;
    double norm()    const { return std::sqrt(norm_sq()); }
};


template <Uint _Size>
template <Uint _Size2>
inline SVectorCL<_Size>::SVectorCL(SVectorCL<_Size2> arg1, double arg2) : base_type(Uninitialized)
{
    if (_Size != _Size2 + 1) // <- this is optimized away
    {
         throw "SVectorCL-classes have wrong size for constructor";
        //throw DROPSErrCL("SVectorCL-classes have wrong size for constructor");
    }
    for (Uint i = 0; i < _Size2; i++)
        (*this)[i] = arg1[i];
    (*this)[_Size2] = arg2;
}

template <Uint _Size>
inline SVectorCL<_Size>::SVectorCL(double arg1, double arg2, double arg3) : base_type(Uninitialized)
{
    if (_Size != 3) // <- this is optimized away
        throw "SVectorCL-class has wrong size for constructor, number of args : 3";
        //throw DROPSErrCL("SVectorCL-class has wrong size for constructor, number of args : 3");
    (*this)[0] = arg1;
    (*this)[1] = arg2;
    (*this)[2] = arg3;
}

template <Uint _Size>
inline SVectorCL<_Size>::SVectorCL(double arg1, double arg2, double arg3, double arg4) : base_type(Uninitialized)
{
    if (_Size != 4) // <- this is optimized away
        throw "SVectorCL-class has wrong size for constructor, number of args : 4";
        //throw DROPSErrCL("SVectorCL-class has wrong size for constructor, number of args : 4");
    (*this)[0] = arg1;
    (*this)[1] = arg2;
    (*this)[2] = arg3;
    (*this)[3] = arg4;
}



template <Uint _Size>
SVectorCL<_Size>&
SVectorCL<_Size>::operator+=(const SVectorCL& v)
{
    for (Uint i = 0; i != _Size; ++i) (*this)[i] += v[i];
    return *this;
}

template <Uint _Size>
SVectorCL<_Size>&
SVectorCL<_Size>::operator-=(const SVectorCL<_Size>& v)
{
    for (Uint i = 0; i != _Size; ++i) (*this)[i] -= v[i];
    return *this;
}

template <Uint _Size>
SVectorCL<_Size>&
SVectorCL<_Size>::operator*=(const SVectorCL<_Size>& v)
{
    for (Uint i = 0; i != _Size; ++i) (*this)[i] *= v[i];
    return *this;
}

template <Uint _Size>
SVectorCL<_Size>&
SVectorCL<_Size>::operator/=(const SVectorCL<_Size>& v)
{
    for (Uint i = 0; i != _Size; ++i) (*this)[i] /= v[i];
    return *this;
}

template <Uint _Size>
SVectorCL<_Size>&
SVectorCL<_Size>::operator*=(double s)
{
    for (Uint i = 0; i != _Size; ++i) (*this)[i] *= s;
    return *this;
}

template <Uint _Size>
SVectorCL<_Size>&
SVectorCL<_Size>::operator/=(double s)
{
    for (Uint i = 0; i != _Size; ++i) (*this)[i] /= s;
    return *this;
}

template <Uint _Size>
double SVectorCL<_Size>::norm_sq() const
{
    double x = 0.0;
    for (Uint i = 0; i < _Size; ++i) x += (*this)[i] * (*this)[i];
    return x;
}

template <Uint _Size>
SVectorCL<_Size> BaryCenter(const SVectorCL<_Size>& v1, const SVectorCL<_Size>& v2)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = .5 * (v1[i] + v2[i]);
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> ConvexComb(double a,
    const SVectorCL<_Size>& v1,
    const SVectorCL<_Size>& v2)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = (1.0 - a) * v1[i] + a * v2[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator+(const SVectorCL<_Size>& v1,
    const SVectorCL<_Size>& v2)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = v1[i] + v2[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator-(const SVectorCL<_Size>& v1,
    const SVectorCL<_Size>& v2)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = v1[i] - v2[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator-(const SVectorCL<_Size>& v1)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = -v1[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator*(double d, const SVectorCL<_Size>& v)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = d * v[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator*(const SVectorCL<_Size>& v, double d)
{
    return d * v;
}

template <Uint _Size>
double inner_prod(const SVectorCL<_Size>& v1, const SVectorCL<_Size>& v2)
{
    double ret = 0.0;
    for (Uint i = 0; i < _Size; ++i) ret += v1[i] * v2[i];
    return ret;
}

inline double
inner_prod(const SVectorCL<3u>& v1, const SVectorCL<3u>& v2)
{
    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}

template <Uint _Size>
SVectorCL<_Size> operator/(const SVectorCL<_Size>& v, double d)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = v[i] / d;
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator*(const SVectorCL<_Size>& v1,
    const SVectorCL<_Size>& v2)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = v1[i] * v2[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> operator/(const SVectorCL<_Size>& v1,
    const SVectorCL<_Size>& v2)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = v1[i] / v2[i];
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> sqrt(const SVectorCL<_Size>& v)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = std::sqrt(v[i]);
    return tempv;
}

template <Uint _Size>
SVectorCL<_Size> fabs(const SVectorCL<_Size>& v)
{
    SVectorCL<_Size> tempv(Uninitialized);
    for (Uint i = 0; i < _Size; ++i) tempv[i] = std::fabs(v[i]);
    return tempv;
}


template <Uint _Size>
std::ostream& operator<<(std::ostream& os, const SVectorCL<_Size>& v)
{
    //    os << v.size() << "    ";
    for (Uint i = 0; i < v.size(); ++i)
        os << v[i] << ' ';
    return os;
}


inline void
cross_product(Point3DCL& res, const Point3DCL& v0, const Point3DCL& v1)
// res= v0 x v1
{
    res[0] = v0[1] * v1[2] - v0[2] * v1[1];
    res[1] = v0[2] * v1[0] - v0[0] * v1[2];
    res[2] = v0[0] * v1[1] - v0[1] * v1[0];
}


//taken from thesis hollasch
inline void
cross_product(Point4DCL& res, const Point4DCL& u, const Point4DCL& v, const Point4DCL& w)
// res= X( v1, v2, v3); 
{
    const double a = v[0] * w[1] - v[1] * w[0];
    const double b = v[0] * w[2] - v[2] * w[0];
    const double c = v[0] * w[3] - v[3] * w[0];
    const double d = v[1] * w[2] - v[2] * w[1];
    const double e = v[1] * w[3] - v[3] * w[1];
    const double f = v[2] * w[3] - v[3] * w[2];
    res[0] = u[1] * f - u[2] * e + u[3] * d;
    res[1] = -u[0] * f + u[2] * c - u[3] * b;
    res[2] = u[0] * e - u[1] * c + u[3] * a;
    res[3] = -u[0] * d + u[1] * b - u[2] * a;
}

//taken from thesis hollasch
inline Point4DCL
cross_product(const Point4DCL& u, const Point4DCL& v, const Point4DCL& w)
// res= X( v1, v2, v3); 
{
    Point4DCL res;
    cross_product(res, u, v, w);
    return res;
}


inline double
l2norm_of_3D_vector(Point4DCL& in)
{
    return std::sqrt(in[0] * in[0] + in[1] * in[1] + in[2] * in[2]);
}

template <Uint _Size>
inline double
normalize_l2norm(SVectorCL<_Size>& in)
{
    double ret = in.norm();
    in /= ret;
    return ret;
}


// std_basis<n>(0)==Null, std_basis<n>(i)[j]==Delta_i-1_j
template <Uint _Size>
inline SVectorCL<_Size> std_basis(Uint i)
{
    SVectorCL<_Size> ret(0.);
    if (i > 0) ret[i - 1] = 1.;
    return ret;
}

template<class T>
SArrayCL<T, 2>
MakeSArray(T a, T b)
{
    SArrayCL<T, 2> ret(Uninitialized);
    ret[0] = a; ret[1] = b;
    return ret;
}

template<class T>
SArrayCL<T, 3>
MakeSArray(T a, T b, T c)
{
    SArrayCL<T, 3> ret(Uninitialized);
    ret[0] = a; ret[1] = b; ret[2] = c;
    return ret;
}

template<class T>
SArrayCL<T, 4>
MakeSArray(T a, T b, T c, T d)
{
    SArrayCL<T, 4> ret(Uninitialized);
    ret[0] = a; ret[1] = b; ret[2] = c; ret[3] = d;
    return ret;
}

template<class T>
SArrayCL<T, 5>
MakeSArray(T a, T b, T c, T d, T e)
{
    SArrayCL<T, 5> ret(Uninitialized);
    ret[0] = a; ret[1] = b; ret[2] = c; ret[3] = d; ret[4] = e;
    return ret;
}

template <class T, Uint _Size>
std::ostream& operator<<(std::ostream& os, const SArrayCL<T, _Size>& a)
{
    //    os << v.size() << "    ";
    for (Uint i = 0; i < a.size(); ++i)
        os << a[i] << ' ';
    return os;
}

//**************************************************************************
// Class:    SMatrixCL                                                      *
// Purpose:  an matrix that remembers its size                              *
// Remarks:  All functions are inline, should be as fast as a "bare" matrix *
//**************************************************************************
template <Uint _Rows, Uint _Cols>
class SMatrixCL : public SVectorCL<_Rows* _Cols>
{
public:
    typedef SVectorCL<_Rows* _Cols> _vec_base;

    SMatrixCL() {}
    explicit           SMatrixCL(InitStateT i) : _vec_base(i) {}
    explicit           SMatrixCL(double val) : _vec_base(val) {}
    template<class In> explicit SMatrixCL(In start) : _vec_base(start) {}
    template<class In> SMatrixCL(In start, In end) : _vec_base(start, end) {}

    // Schreib- & Lesezugriff
    double& operator() (int row, int col) { return (*this)[row * _Cols + col]; }// Matrix(i,j)
    double  operator() (int row, int col) const { return (*this)[row * _Cols + col]; }

    SVectorCL<_Rows> col(int) const;
    void             col(int, const SVectorCL<_Rows>&);

    // Zuweisung & Co.
    SMatrixCL& operator+=(const SMatrixCL&);                // Matrix=Matrix+Matrix'
    SMatrixCL& operator-=(const SMatrixCL&);                // Matrix=Matrix-Matrix'
    SMatrixCL& operator*=(double s);                               // Matrix = c * Matrix
    SMatrixCL& operator/=(double s);                               // Matrix = Matrix'/c

// Dimensionen feststellen
    Uint num_rows() const { return _Rows; }                        // Zeilenzahl
    Uint num_cols() const { return _Cols; }                        // Spaltenzahl
};

template<Uint _Rows, Uint _Cols>
SVectorCL<_Rows>
SMatrixCL<_Rows, _Cols>::col(int c) const
{
    SVectorCL<_Rows> ret(Uninitialized);
    for (Uint i = 0; i != _Rows; ++i, c += _Cols)
        ret[i] = (*this)[c];
    return ret;
}

template<Uint _Rows, Uint _Cols>
void
SMatrixCL<_Rows, _Cols>::col(int c, const SVectorCL<_Rows>& v)
{
    for (Uint i = 0; i != _Rows; ++i)
        (*this)(i, c) = v[i];
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>&
SMatrixCL<_Rows, _Cols>::operator+=(const SMatrixCL<_Rows, _Cols>& m)
{
    *static_cast<_vec_base*>(this) += *static_cast<const _vec_base*>(&m);
    return *this;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>&
SMatrixCL<_Rows, _Cols>::operator-=(const SMatrixCL<_Rows, _Cols>& m)
{
    *static_cast<_vec_base*>(this) -= *static_cast<const _vec_base*>(&m);
    return *this;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>&
SMatrixCL<_Rows, _Cols>::operator*=(double d)
{
    *static_cast<_vec_base*>(this) *= d;
    return *this;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>&
SMatrixCL<_Rows, _Cols>::operator/=(double d)
{
    *static_cast<_vec_base*>(this) /= d;
    return *this;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>
operator+(const SMatrixCL<_Rows, _Cols>& m1, const SMatrixCL<_Rows, _Cols>& m2)
{
    SMatrixCL<_Rows, _Cols> ret(Uninitialized);
    *static_cast<typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&ret)
        = *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m1)
        + *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m2);
    return ret;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>
operator-(const SMatrixCL<_Rows, _Cols>& m1, const SMatrixCL<_Rows, _Cols>& m2)
{
    SMatrixCL<_Rows, _Cols> ret(Uninitialized);
    *static_cast<typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&ret)
        = *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m1)
        - *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m2);
    return ret;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>
operator-(const SMatrixCL<_Rows, _Cols>& m)
{
    SMatrixCL<_Rows, _Cols> ret(Uninitialized);
    *static_cast<typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&ret)
        = -*static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m);
    return ret;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>
operator*(double d, const SMatrixCL<_Rows, _Cols>& m)
{
    SMatrixCL<_Rows, _Cols> ret(Uninitialized);
    *static_cast<typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&ret)
        = d * *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m);
    return ret;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>
operator*(const SMatrixCL<_Rows, _Cols>& m, double d)
{
    SMatrixCL<_Rows, _Cols> ret(Uninitialized);
    *static_cast<typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&ret)
        = *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m) * d;
    return ret;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Rows, _Cols>
operator/(const SMatrixCL<_Rows, _Cols>& m, double d)
{
    SMatrixCL<_Rows, _Cols> ret(Uninitialized);
    *static_cast<typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&ret)
        = *static_cast<const typename SMatrixCL<_Rows, _Cols>::_vec_base*>(&m) / d;
    return ret;
}

template<Uint _RowsL, Uint _ColsR, Uint _Dim>
SMatrixCL<_RowsL, _ColsR>
operator*(const SMatrixCL<_RowsL, _Dim>& m1, const SMatrixCL<_Dim, _ColsR>& m2)
{
    SMatrixCL<_RowsL, _ColsR> ret(0.0);
    for (Uint row = 0; row != _RowsL; ++row)
        for (Uint col = 0; col != _ColsR; ++col)
            for (Uint i = 0; i != _Dim; ++i)
                ret(row, col) += m1(row, i) * m2(i, col);
    return ret;
}

template<Uint _Rows, Uint _Cols>
SMatrixCL<_Cols, _Cols>
GramMatrix(const SMatrixCL<_Rows, _Cols>& m)
/// Computes m^T*m
{
    SMatrixCL<_Cols, _Cols> ret(0.0);
    for (Uint row = 0; row != _Cols; ++row) {
        for (Uint col = 0; col < row; ++col) {
            for (Uint i = 0; i != _Rows; ++i)
                ret(row, col) += m(i, row) * m(i, col);
            ret(col, row) = ret(row, col);
        }
        for (Uint i = 0; i != _Rows; ++i)
            ret(row, row) += std::pow(m(i, row), 2);
    }
    return ret;
}

template<Uint _Rows, Uint _Cols>
SVectorCL<_Cols>
transp_mul(const SMatrixCL<_Rows, _Cols>& m, const SVectorCL<_Rows>& v)
{
    SVectorCL<_Cols> ret(0.0);
    for (Uint col = 0; col != _Cols; ++col)
        for (Uint i = 0; i != _Rows; ++i)
            ret[col] += m(i, col) * v[i];
    return ret;
}

template<Uint _Rows, Uint _Cols>
SVectorCL<_Rows>
operator*(const SMatrixCL<_Rows, _Cols>& m, const SVectorCL<_Cols>& v)
{
    SVectorCL<_Rows> ret(0.0);
    for (Uint row = 0; row != _Rows; ++row)
        for (Uint i = 0; i != _Cols; ++i)
            ret[row] += m(row, i) * v[i];
    return ret;
}

inline SVectorCL<3>
operator*(const SMatrixCL<3, 3>& m, const SVectorCL<3>& v)
{
    SVectorCL<3> ret(Uninitialized);
    const double* const a = m.begin();
    ret[0] = a[0] * v[0] + a[1] * v[1] + a[2] * v[2];
    ret[1] = a[3] * v[0] + a[4] * v[1] + a[5] * v[2];
    ret[2] = a[6] * v[0] + a[7] * v[1] + a[8] * v[2];
    return ret;
}

inline SVectorCL<4>
operator*(const SMatrixCL<4, 4>& m, const SVectorCL<4>& v)
{
    SVectorCL<4> ret(Uninitialized);
    const double* const a = m.begin();
    ret[0] = a[0] * v[0] + a[1] * v[1] + a[2] * v[2] + a[3] * v[3];
    ret[1] = a[4] * v[0] + a[5] * v[1] + a[6] * v[2] + a[7] * v[3];
    ret[2] = a[8] * v[0] + a[9] * v[1] + a[10] * v[2] + a[11] * v[3];
    ret[3] = a[12] * v[0] + a[13] * v[1] + a[14] * v[2] + a[15] * v[3];
    return ret;
}


template <Uint Rows, Uint Cols>
double inner_prod(const SVectorCL<Rows>& v1, const SMatrixCL<Rows, Cols>& m, const SVectorCL<Cols>& v2)
{
    double ret = 0.0;
    for (Uint i = 0; i < Rows; ++i)
        for (Uint j = 0; j < Cols; ++j)
            ret += v1[i] * m(i, j) * v2[j];
    return ret;
}

template <Uint _Rows, Uint _Cols>
std::ostream& operator << (std::ostream& os, const SMatrixCL<_Rows, _Cols>& m)
{
    const Uint M = m.num_rows();
    const Uint N = m.num_cols();

    os << M << ' ' << N << '\n';

    for (Uint row = 0; row < M; ++row)
    {
        os << "  ";
        for (Uint col = 0; col < N; ++col)
            os << m(row, col) << ' ';
        os << '\n';
    }
    return os;
}

template <Uint _Rows>
inline SMatrixCL<_Rows, _Rows>
outer_product(const SVectorCL<_Rows>& a, const SVectorCL<_Rows>& b)
{
    SMatrixCL<_Rows, _Rows> ret(Uninitialized);
    for (Uint i = 0; i < _Rows; ++i)
        for (Uint j = 0; j < _Rows; ++j)
            ret(i, j) = a[i] * b[j];
    return ret;
}

template <Uint _Rows>
inline double
frobenius_norm_sq(const SMatrixCL<_Rows, _Rows>& a)
{
    double ret = 0;
    for (Uint i = 0; i < _Rows * _Rows; ++i)
        ret += a[i] * a[i];
    return ret;
}

template <Uint _Rows>
inline double
trace(const SMatrixCL<_Rows, _Rows>& a)
{
    double ret = 0.;
    for (Uint i = 0; i < _Rows; ++i)
        ret += a(i, i);
    return ret;
}

template <Uint _Rows>
inline SMatrixCL<_Rows, _Rows>&
assign_transpose(SMatrixCL<_Rows, _Rows>& out, const SMatrixCL<_Rows, _Rows>& in)
{
    for (Uint i = 0; i < _Rows; ++i) {
        for (Uint j = 0; j < i; ++j) {
            out(i, j) = in(j, i);
            out(j, i) = in(i, j);
        }
        out(i, i) = in(i, i);
    }
    return out;
}

template <Uint _Rows, Uint _Cols>
inline SMatrixCL<_Rows, _Cols> eye()
{
    SMatrixCL<_Rows, _Cols> ret;
    const Uint m = std::min(_Rows, _Cols);
    for (Uint i = 0; i < m; ++i)
        ret(i, i) = 1.;
    return ret;
}


一个主程序告诉你方法怎么用

不解释,看代码。

#pragma once
/// \file container.h
/// \brief Test for container.h
/// \author LSEC: Song Lu, lusong@lsec.cc.ac.cn

#include "container.h" 
int main()
{
	/// <summary>
	/// SArrayCL 测试
	/// </summary>
	SArrayCL<double,3> myArray0(1.0);//初始化方式1,全部赋值为一个数
	std::vector<double> vec = { 1.0,2.0,3.0 };
	SArrayCL<double, 3> myArray1(vec.begin());//初始化方式2,通过已有量赋值
	SArrayCL<double, 3> myArray2(vec.begin(),vec.end());//初始化方式2,通过已有量赋值
	myArray1.begin();//得到开始位置的迭代器
	myArray1.end();//得到末尾位置的迭代器
	std::cout << myArray1[1] << std::endl;//可以通过中括号取值
	std::cout << myArray1.size() << std::endl;//返回大小
	bool ifEqual = (myArray1 == myArray2);//通过等于号判断是否相等
	bool ifSmall= myArray0 < myArray2;//通过小于号比较大小
	Addr(myArray0);//Addr 函数可以返回地址,即myArray0.begin();
	SArrayCL<double, 2> myArray3 = MakeSArray(double(1), double(2));//生成数组,支持长度为2、3、4、5的生成
#if 0
	myArray1[100];//下标异常测试
#endif


	/// <summary>
	/// SVectorCL 测试,SVectorCL 继承自 SArrayCL,自测试独有方法
	/// </summary>
	SVectorCL<2> myVector0(1.0);
	SVectorCL<3> myVector1(1.0,2.0,3.0);//直接构造,只支持个数为3或者4
	SVectorCL<4> myVector2(myVector1, 1.0);//支持向量加一个数的扩充构造
	SVectorCL<3> myVector3(vec.begin());//同样支持拷贝构造
	SVectorCL<3> myVector4(vec.begin(), vec.end());//同样支持拷贝构造
	myVector3 += myVector4;//重载的 += 的使用,-=、*=、/= 类似
	myVector3 *= 3.0;//*= 和 /= 还支持和一个数的运算 v = 3*v
	myVector3.norm();//求2范数
	myVector3.norm_sq();//求2范数的平方
	BaryCenter(myVector3, myVector4);//求两点中心位置
	ConvexComb(1/3,myVector3, myVector4);//求凸包(1.0 - a) * v1 + a * v2;
	SVectorCL<3> myVector5 = myVector4 + myVector3;//加法+ -类似
	SVectorCL<3> myVector6 = myVector5 * 1.0;//和数的乘法,这里的数可以放前面,也可以放后面
	double ret = inner_prod(myVector5, myVector6);//两个的内积
	SVectorCL<3> myVector7 = sqrt(myVector6);//对每个分量开根号
	SVectorCL<3> myVector8 = fabs(myVector6);//对每个分量求绝对值
	std::cout << myVector6 << std::endl;//重载了<< 标准输出
	Point3DCL res;
	cross_product(res, myVector3, myVector4);//三维向量的叉积
	Point4DCL res2;
	cross_product(res2, myVector2, myVector2, myVector2);//四维向量的叉积
	res2 = cross_product(myVector2, myVector2, myVector2);//四维向量叉积表达方式2
	double l2norm = l2norm_of_3D_vector(myVector2);//四维向量的 l2 模
	double l2norm2 = normalize_l2norm(myVector2);//myVector2 归一化,即把长度变成1,并返回模
	SVectorCL<5> myVector9 = std_basis<5>(2);//生成第2个位置是1,别的位置是0的五维向量
	SMatrixCL<3, 3> Mop = outer_product(myVector1, myVector1);//计算外积,即 x^T*x,变成矩阵


	/// <summary>
	/// SMatrixCL 测试,继承自 SVectorCL
	/// </summary>
	/// <returns></returns>
	SMatrixCL<2, 2> M0(1.0);//初始化方式1
	std::vector<double> vecM = { 1.0,2.0,3.0,4.0};
	SMatrixCL<2, 2> M1(vecM.begin());//初始化方式2,按行优先,即一行一行往下排,这和MATLAB是不太一样的
	SMatrixCL<2, 2> M2(vecM.begin(),vecM.end());//初始化方式3,按行优先
	std::cout << M2(1, 0) << std::endl;//矩阵可以使用圆括号取值
	SVectorCL<2> v = M1.col(1);//取第2列
	M1.col(1, v);//取第2列赋值给v
	Uint rowNum = M2.num_rows();//返回行数
	Uint colNum = M2.num_cols();//返回列数
	SMatrixCL<2, 2> G = GramMatrix(M1);//返回A^TA
	SMatrixCL<2, 2> Mtrans;
	assign_transpose(Mtrans, M2);//求矩阵M2的转置
	SMatrixCL<2, 2> D = eye<2, 2>();//生成对角矩阵,如果不是方的,别的部分为0
	std::cout << M1 << std::endl;//一行一行输出,重载了<<
	double r_fnsq = frobenius_norm_sq(M2);//计算frobenius范数的平方
	double tr = trace(M0);//求矩阵的迹,即对角线的和
	M2 += M1;//-= 含义类似
	M2 *= 0.3;// /= 含义类似
	M0 = M1 - M2;//+ 含义类似
	M0 =  - M2;//取反操作
	M0 = 2.0 * M1;//乘积,反过来 M1*2.0 也行,除法类似
	M0 = M1 * M1;//就是一般的矩阵乘法,C=AB



	/// <summary>
	/// 矩阵和向量之间的运算测试
	/// </summary>
	SVectorCL<2> v2 = transp_mul(M1, myVector0);// 返回A^T*v的结果,即v和A的每一列分别内积得到的向量,也可以理解为v^TA,向量乘矩阵
	double r = inner_prod(myVector0, M0, myVector0);//放回x^TAx的结果
	Mop(1, 0) = 1;
	SVectorCL<3> v3 = Mop * myVector1;//矩阵乘向量Ax
	

	/// <returns></returns>
	return 0;
}
相关推荐
©️2020 CSDN 皮肤主题: 酷酷鲨 设计师:CSDN官方博客 返回首页