线性代数库

头文件

/*
 * linear_algebra.h
 *
 *  Created on: 2013年9月1日
 *      Author: sai
 * last modify: 2013年12月17日
 */

#ifndef LINEAR_ALGEBRA_H_
#define LINEAR_ALGEBRA_H_

#include <iostream>
#include <map>

class Matrix;
/******************************************************************************
 * 向量类
 * 虽然按行存,但默认是列向量
******************************************************************************/
class Vector
{
public:
    friend class Matrix;
    Vector();
    Vector(int n); //构造函数, 指定维数
    Vector(int n, double v); //用v填充
    Vector(int n, double* a, bool transfer=false); //构造函数,用数组构造
    virtual ~Vector(); //析构函数
    Vector(const Vector& vec); //拷贝构造函数
    Vector& operator = (const Vector& vec); //赋值函数
    Vector& operator = (const Matrix& m); //将1*n的矩阵赋给该向量
    void print(std::ostream& out) const;

    double& operator [] (int index); //使用
    double& operator [] (int index) const; //使用
    double& operator () (int index); //使用
    double& operator () (int index) const; //使用
    int size() const; //获取维数
    void resize(int n); //改变大小
    
    void random(double length=1, double from=0);
    void zero();
    double normalize();
    double max() const;
    double min() const;

    Vector& operator *= (double c); //数乘
    Vector operator * (double c) const; //数乘
    friend Vector operator * (double c, const Vector& vec); //数乘
    Vector& operator += (const Vector& vec); //向量+=
    Vector operator + (const Vector& vec) const; //向量加法
    Vector& operator -= (const Vector& vec); //向量-=
    Vector operator - (const Vector& vec) const; //向量减法
    double operator * (const Vector& vec) const; //内积
    Matrix operator * (const Matrix& vec) const; //向量*矩阵
    Matrix transpose() const; //转置为1行n列矩阵
public:
    int dim; //维数
    double* data; //数据
    int flag; //1需要释放但不传递,0需要释放传递给继任者,-1不需要释放
};

std::ostream& operator << (std::ostream& oo, const Vector& vec); //重载插入运算符


/******************************************************************************
 * 矩阵类
******************************************************************************/
class Matrix
{
public:
    friend class Vector;
    Matrix();
    Matrix(int n, int m); //构造函数,指定行数和列数
    Matrix(int n); //构造函数, 指定行数和列数均为n
    Matrix(const Matrix& m); //拷贝构造函数
    Matrix& operator = (const Matrix& m); //赋值函数
    virtual ~Matrix();
    void resize(int n, int m); //改变大小
    void print(std::ostream& out) const;
    Vector operator [] (int index); //使用, 返回该行向量的引用
    Vector operator [] (int index) const; //使用, 返回该行向量的引用
    double& operator () (int r, int c); //使用

    void unit(int n); //单位方阵
    void zero(); //0矩阵

    Matrix& operator *= (double c); //数乘
    Matrix operator * (double c) const; //数乘
    friend Matrix operator * (double c, const Matrix& m); //数乘
    Matrix operator + (const Matrix& mat) const; //矩阵加法
    Matrix& operator += (const Matrix& mat); //矩阵+=
    Matrix operator - (const Matrix& mat) const; //矩阵减法
    Matrix& operator -= (const Matrix& mat); //矩阵-=
    Vector operator * (const Vector& vec) const; //矩阵*向量
    Matrix operator * (const Matrix& mat) const; //矩阵乘法
    double trace() const; //迹
    Matrix transpose() const; //转置
    Matrix inverse() const; //逆
    double det() const; //行列式
    void QR(Matrix& Q, Matrix& R) const; //QR分解
    Vector solve(const Vector& b) const; //解方程
public:
    int row; //行数
    int col; //列数
    double** data; //数据
    int flag; //1需要释放但不传递,0需要释放传递给继任者,-1不需要释放
};

std::ostream& operator << (std::ostream& oo, const Matrix& m); //重载插入运算符


/******************************************************************************
 * 稀疏向量类
******************************************************************************/
class SparseVector
{
public:
    SparseVector();

    double& operator [] (int index); //使用
public:
    std::map<int, double> data;
};


/******************************************************************************
 * 稀疏矩阵类
******************************************************************************/
class SparseMatrix
{
private:
    struct Pos
    {
        int r;
        int c;
        Pos(int rr, int cc) : r(rr), c(cc) {}
        bool operator < (const Pos& t) const
        {
            return r < t.r || r==t.r && c<t.c;
        }
    };
public:
    SparseMatrix();

    double& operator () (int r, int c); //使用
public:
    std::map<Pos, double> data;
};



#endif /* LINEAR_ALGEBRA_H_ */

 

源文件

/*
 * linear_algebra.cpp
 *
 *  Created on: 2013年12月17日
 *      Author: sai
 */

#include "linalgebra.h"
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cfloat>
using namespace std;

const double ZERO = 1e-4;

Matrix null;

/******************************************************************************
 * 向量类
******************************************************************************/

//构造函数
Vector::Vector()
{
    dim = 0;
    data = NULL;
    flag = 1;
}

//构造函数
Vector::Vector(int n)
{
    if (n < 0)
    {
        cout<<"Vector构造函数参数错误.\n";
        return;
    }
    dim = n;
    data = new double[n];
    flag = 1;
}

Vector::Vector(int n, double v)
{
    if (n < 0)
    {
        cout<<"Vector构造函数参数错误.\n";
        return;
    }
    dim = n;
    data = new double[n];
    for (int i=0; i<dim; i++)
        data[i] = v;
    flag = 1;
}

//构造函数
Vector::Vector(int n, double* a, bool transfer)
{
    if (n < 0)
    {
        cout<<"Vector构造函数参数错误.\n";
        return;
    }
    dim = n;
    if (transfer) //浅拷贝
    {
        data = a;
        flag = -1;
    }
    else //深拷贝
    {
        data = new double[n];
        memcpy(data, a, n*sizeof(double));
        flag = 1;
    }
}

//析构函数
Vector::~Vector()
{
    if (data && flag > -1)
        delete []data;
    data = NULL;
}

//拷贝构造
Vector::Vector(const Vector& vec)
{
    dim = vec.dim;    
    if (vec.flag == 1)
    {
        data = new double[dim];
        memcpy(data, vec.data, dim*sizeof(double));
        flag = 1;
    }
    else
    {
        data = vec.data;
        if (vec.flag == 0)
        {
            (const_cast<Vector*>(&vec))->flag = -1;
            flag = 0;
        }
        else
            flag = -1;
    }
}

//赋值
Vector& Vector::operator = (const Vector& vec)
{
    if (data && flag > -1)
        delete []data;
    dim = vec.dim;
    if (vec.flag == 1)
    {
        data = new double[dim];
        memcpy(data, vec.data, dim*sizeof(double));
        flag = 1;
    }
    else if (vec.flag == 0) //需要传递
    {
        data = vec.data;
        flag = 1; //仅传递一次
        (const_cast<Vector*>(&vec))->flag = -1;
    }
    else //不需要传递的,将flag置为不需要释放
    {
        data = vec.data;
        flag = -1;
    }
    return *this;
}

//赋值
Vector& Vector::operator = (const Matrix& m)
{
    if (m.row != 1)
    {
        cerr<<"无法赋值, operator = (const Matrix& m)"<<endl;
        return *this;
    }
    if (data && flag > -1)
        delete []data;
    dim = m.col;
    data = new double[dim];
    memcpy(data, m.data[0], dim*sizeof(double));
    flag = 1;
    return *this;
}

void Vector::print(ostream& out) const
{
    for (int i=0; i<dim; i++)
        out<<data[i]<<' ';
}

//重载插入运算符
ostream& operator << (ostream& oo, const Vector& vec)
{
    vec.print(oo);
    return oo;
}

//数组下标引用
double& Vector::operator [] (int index)
{
    return data[index];
}

//数组下标引用
double& Vector::operator [] (int index) const
{
    return data[index];
}

//括号下标引用
double& Vector::operator () (int index)
{
    return data[index];
}

//括号下标引用
double& Vector::operator () (int index) const
{
    return data[index];
}

//获得维数
int Vector::size() const
{
    return dim;
}

//改变维数
void Vector::resize(int n)
{
    if (n != dim)
    {
        if (data && flag > -1)
            delete []data;
        dim = n;
        data = new double[dim];
        flag = 1;
    }
}

//赋随机值
void Vector::random(double length, double from)
{
    for (int i=0; i<dim; i++)
        data[i] = length * rand() / RAND_MAX + from;
}

//全部为0
void Vector::zero()
{
    memset(data, 0, dim*sizeof(double));
}

double Vector::normalize()
{
    double norm = 0;
    for (int i=0; i<dim; i++)
        norm += data[i] * data[i];
    norm = sqrt(norm);
    for (int i=0; i<dim; i++)
        data[i] /= norm;
    return norm;
}

double Vector::max() const
{
    double ans = -DBL_MAX;
    for (int i=0; i<dim; i++)
        ans = data[i] > ans ? data[i] : ans;
    return ans;
}

double Vector::min() const
{
    double ans = DBL_MAX;
    for (int i=0; i<dim; i++)
        ans = data[i] < ans ? data[i] : ans;
    return ans;
}

//数乘
Vector& Vector::operator *= (double c)
{
    for (int i=0; i<dim; i++)
        data[i] *= c;
    return *this;
}

//数乘
Vector Vector::operator *(double c) const
{
    Vector ans(dim);
    ans.flag = 0;
    for (int i=0; i<dim; i++)
        ans.data[i] = c * data[i];
    return ans;
}

//数乘
Vector operator * (double c, const Vector& vec)
{
    Vector ans(vec.dim);
    ans.flag = 0;
    for (int i=0; i<vec.dim; i++)
        ans.data[i] = c * vec.data[i];
    return ans;
}

//+=重载
Vector& Vector::operator +=(const Vector& vec)
{
    for (int i=0; i<dim; i++)
        data[i] += vec.data[i];
    return *this;
}

//+重载
Vector Vector::operator +(const Vector& vec) const
{
    Vector ans(dim);
    ans.flag = 0;
    for (int i=0; i<dim; i++)
        ans.data[i] = data[i] + vec.data[i];
    return ans;
}

//-=重载
Vector& Vector::operator -=(const Vector& vec)
{
    for (int i=0; i<dim; i++)
        data[i] -= vec.data[i];
    return *this;
}

//-重载
Vector Vector::operator -(const Vector& vec) const
{
    Vector ans(dim);
    ans.flag = 0;
    for (int i=0; i<dim; i++)
        ans.data[i] = data[i] - vec.data[i];
    return ans;
}

//内积
double Vector::operator * (const Vector& vec) const
{
    double ans = 0;
    for (int i=0; i<dim; i++)
        ans += data[i] * vec.data[i];
    return ans;
}

//向量*矩阵, (n,1) * (1,n)
Matrix Vector::operator * (const Matrix& m) const
{
    if (m.row != 1)
    {
        cerr<<"无法相乘, Vector::operator * (const Matrix& m)"<<endl;
        return Matrix();
    }
    Matrix ans(dim, m.col);
    ans.flag = 0;
    for (int i=0; i<dim; i++)
        for (int j=0; j<m.col; j++)
            ans.data[i][j] = data[i] * m.data[0][j];
    return ans;
}

//转置为1*n矩阵
Matrix Vector::transpose() const
{
    Matrix ans(1, dim);
    ans.flag = 0;
    memcpy(ans.data[0], data, dim*sizeof(double));
    return ans;
}




/******************************************************************************
 * 稀疏向量类
******************************************************************************/
SparseVector::SparseVector()
{
}

double& SparseVector::operator [] (int index)
{
    return data[index];
}





/******************************************************************************
 * 矩阵类
******************************************************************************/
Matrix::Matrix()
{
    row = 0;
    col = 0;
    data = NULL;
    flag = 1;
}

Matrix::Matrix(int n, int m)
{
    if (n <= 0 || m <= 0)
    {
        cout<<"Matrix构造函数参数错误\n";
        return;
    }
    row = n;
    col = m;
    data = new double*[n];
    for (int i=0; i<n; i++)
        data[i] = new double[m];
    flag = 1;
}

Matrix::Matrix(int n)
{
    if (n <= 0)
    {
        cout<<"Matrix构造函数参数错误\n";
        return;
    }
    row = n;
    col = n;
    data = new double*[n];
    for (int i=0; i<n; i++)
        data[i] = new double[n];
    flag = 1;
}

Matrix::~Matrix()
{
    if (data != NULL && flag > -1)
    {
        for (int i=0; i<row; i++)
            delete []data[i];
        delete []data;
    }
    data = NULL;
}

//拷贝构造
Matrix::Matrix(const Matrix& m)
{
    row = m.row;
    col = m.col;
    if (m.flag == 1)
    {
        data = new double*[row];
        for (int i=0; i<row; i++)
        {
            data[i] = new double[col];
            memcpy(data[i], m.data[i], col*sizeof(double));
        }
    }
    else
    {
        data = m.data;
        if (m.flag == 0)
        {
            (const_cast<Matrix*>(&m))->flag = -1;
            flag = 0;
        }
        else
            flag = -1;
    }
}

//赋值
Matrix& Matrix::operator = (const Matrix& m)
{
    if (data && flag > -1)
    {
        for (int i=0; i<row; i++)
            delete []data[i];
        delete []data;
    }
    row = m.row;
    col = m.col;
    if (m.flag == 1)
    {
        data = new double*[row];
        for (int i=0; i<row; i++)
        {
            data[i] = new double[col];
            memcpy(data[i], m.data[i], col*sizeof(double));
        }
        flag = 1;
    }
    else if (m.flag == 0) //若原Matrix需要释放,拷贝构造将释放传递给新Matrix
    {
        data = m.data;
        flag = 1;
        (const_cast<Matrix*>(&m))->flag = -1;
    }
    else
    {
        data = m.data;
        flag = -1;
    }
    return *this;
}

void Matrix::print(ostream& out) const
{
    for(int i=0; i<row; i++)
    {
        for(int j=0; j<col; j++)
            out<<data[i][j]<<' ';
        out<<endl;
    }
}

//重载插入运算符
ostream& operator << (ostream& oo, const Matrix& m)
{
    m.print(oo);
    return oo;
}


Vector Matrix::operator [] (int index)
{
    Vector ans(col, data[index]);
    ans.flag = -1;
    return ans;
}

Vector Matrix::operator [] (int index) const
{
    Vector ans(col, data[index]);
    ans.flag = -1;
    return ans;
}

double& Matrix::operator () (int r, int c)
{
    return data[r][c];
}

void Matrix::resize(int n, int m)
{
    if (row == n && col == m)
        return;
    if (data != NULL)
    {
        for (int i=0; i<row; i++)
            delete []data[i];
        delete []data;
    }
    row = n;
    col = m;
    data = new double*[row];
    for (int i=0; i<row; i++)
        data[i] = new double[col];
    flag = 1;
}

void Matrix::unit(int n)
{
    resize(n, n);
    for (int i=0; i<n; i++)
    {
        data[i][i] = 1;
        for (int j=0; j<i; j++)
        {
            data[i][j] = 0;
            data[j][i] = 0;
        }
    }
}

void Matrix::zero()
{
    for (int i=0; i<row; i++)
        memset(data[i], 0, col*sizeof(double));
}

//数乘
Matrix& Matrix::operator *= (double c)
{
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            data[i][j] *= c;
    return *this;
}

//数乘
Matrix Matrix::operator *(double c) const
{
    Matrix ans(row, col);
    ans.flag = 0;
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            ans.data[i][j] = c * data[i][j];
    return ans;
}

//数乘
Matrix operator * (double c, const Matrix& m)
{
    Matrix ans = m;
    ans.flag = 0;
    for (int i=0; i<m.row; i++)
        for (int j=0; j<m.col; j++)
            ans.data[i][j] *= c;
    return ans;
}

Matrix Matrix::operator + (const Matrix& mat) const
{
    Matrix ans(row, col);
    ans.flag = 0;
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            ans.data[i][j] = data[i][j] + mat.data[i][j];
    return ans;
}

Matrix& Matrix::operator +=(const Matrix& mat)
{
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            data[i][j] += mat.data[i][j];
    return *this;
}

Matrix Matrix::operator - (const Matrix& mat) const
{
    Matrix ans(row, col);
    ans.flag = 0;
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            ans.data[i][j] = data[i][j] - mat.data[i][j];
    return ans;
}

Matrix& Matrix::operator -=(const Matrix& mat)
{
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            data[i][j] -= mat.data[i][j];
    return *this;
}

Vector Matrix::operator * (const Vector& vec) const
{
    if (col != vec.dim)
    {
        cout<<"维数不一致,无法相乘\n";
    }

    Vector ans(row);
    ans.flag = 0;
    for (int i=0; i<row; i++)
    {
        ans.data[i] = 0;
        for (int j=0; j<col; j++)
            ans.data[i] += data[i][j] * vec.data[j];
    }
    return ans;
}

Matrix Matrix::operator * (const Matrix& mat) const
{
    if (col != mat.row)
    {
        cout<<"维数不一致,无法相乘\n";
    }
    Matrix ans(row, mat.col);
    ans.flag = 0;
    for (int i=0; i<row; i++)
        for (int j=0; j<mat.col; j++)
        {
            ans.data[i][j] = 0;
            for (int k=0; k<col; k++)
                ans.data[i][j] += data[i][k] * mat.data[k][j];
        }
    return ans;
}

double Matrix::trace() const
{
    if (row != col)
    {
        cout<<"不是方阵,无法求迹\n";
        return 0;
    }
    double ans = 0;
    for (int i=0; i<row; i++)
        ans += data[i][i];
    return ans;
}

Matrix Matrix::transpose() const
{
    Matrix ans(col, row);
    ans.flag = 0;
    for (int i=0; i<row; i++)
        for (int j=0; j<col; j++)
            ans.data[i][j] = data[j][i];
    return ans;
}

//矩阵逆,Gauss-Jordan列主元法
Matrix Matrix::inverse() const
{
    if(row != col)
    {
        cout<<"不是方阵,无法求逆\n";
        return null;
    }
    Matrix ans;
    ans.flag = 0;
    ans.unit(row);

    for(int i=0; i<row; i++)
    {
        //选主元
        int max = i;
        for(int j=i; j<row; j++)
            if(fabs(data[j][i]) > fabs(data[max][i]))
                max = j;
        if( fabs(data[max][i]) < ZERO )
            return null;

        //换主元
        if(max != i)
        {
            for(int j=i; j<row; j++)
            {
                double temp = data[i][j];
                data[i][j] = data[max][j];
                data[max][j] = temp;
            }
            for(int j=0; j<row; j++)
            {
                double temp = ans.data[i][j];
                ans.data[i][j] = ans.data[max][j];
                ans.data[max][j] = temp;
            }
        }

        //对角线元素化为1
        for(int j=i+1; j<row; j++)
            data[i][j] /= data[i][i];
        for(int j=0; j<row; j++)
            ans.data[i][j] /= data[i][i];

        //消0
        for(int j=i+1; j<row; j++)
        {
            if(j == i)
                continue;
            for(int k=i+1; k<row; k++)
                data[j][k] -= data[i][k] * data[j][i];
            for(int k=0; k<row; k++)
                ans.data[j][k] -= ans.data[i][k] * data[j][i];
        }
    }

    for(int i=row-1; i>0; i--)
        for(int j=0; j<i; j++)
            for(int k=0; k<row; k++)
                ans.data[j][k] -= ans.data[i][k] * data[j][i];
    return ans;
}

//行列式
double Matrix::det() const
{
    if(row != col)
    {
        cout<<"不是方阵,无法求行列式\n";
        return 0;
    }
    bool sgn = true;
    for(int i=0; i<row; i++)
    {
        //选主元
        int maxi = i, maxj = i;
        for(int j=i; j<row; j++)
            for(int k=i; k<col; k++)
                if(fabs(data[j][k]) > fabs(data[maxi][maxj]))
                {
                    maxi = j;
                    maxj = k;
                }
        if(fabs(data[maxi][maxj]) < ZERO)
            return 0;
        //换主元
        if(maxi != i)
        {
            sgn = !sgn;
            for(int j=i; j<col; j++)
            {
                double temp = data[i][j];
                data[i][j] = data[maxi][j];
                data[maxi][j] = temp;
            }
        }
        if(maxj != i)
        {
            sgn = !sgn;
            for(int j=i; j<row; j++)
            {
                double temp = data[j][i];
                data[j][i] = data[j][maxj];
                data[j][maxj] = temp;
            }
        }
        //消0
        for(int j=i+1; j<row; j++)
            for(int k=i+1; k<col; k++)
                data[j][k] -= data[i][k] * data[j][i] / data[i][i];
    }

    double ans = data[0][0];
    for(int i=1; i<row; i++)
        ans *= data[i][i];
    if(sgn)
        return ans;
    else
        return -ans;
}

//QR分解,Schimidt正交化法
void Matrix::QR(Matrix& Q, Matrix& R) const
{
    if(row != col)
    {
        cout<<"不是方阵,无法QR分解.\n";
        return;
    }

    Q = *this;
    R.unit(row);

    Vector lambda(row);
    double mu;

    for(int i=1; i<row; i++)
    {
        lambda.data[i-1] = 0;
        for(int j=0; j<row; j++)
            lambda.data[i-1] += Q.data[j][i-1] * Q.data[j][i-1];
        if(fabs(lambda.data[i-1]) < ZERO)
        {
            Q.resize(0,0);
            R.resize(0,0);
            return;
        }
        for(int j=0; j<i; j++)
        {
            double mu = 0;
            for(int k=0; k<row; k++)
                mu += Q.data[k][i] * Q.data[k][j];
            for(int k=0; k<row; k++)
            {
                Q.data[k][i] -= mu / lambda.data[j] * Q.data[k][j];
                R.data[j][k] += mu / lambda.data[j] * R.data[i][k];
            }
        }
    }
    for(int i=0; i<row; i++)
    {
        mu = 0;
        for(int j=0; j<row; j++)
            mu += Q.data[j][i] * Q.data[j][i];
        mu = sqrt(mu);
        for(int j=0; j<row; j++)
        {
            Q.data[j][i] /= mu;
            R.data[i][j] *= mu;
        }
    }
}

Vector Matrix::solve(const Vector& b) const
{
    Vector ans(b.dim);
    ans.flag = 0;
    return ans;
}



/******************************************************************************
 * 稀疏矩阵类
******************************************************************************/
double& SparseMatrix::operator () (int r, int c)
{
    return data[Pos(r,c)];
}