矩阵及其特殊矩阵

之前我们已经用数组的数组方式写了一个二维数组，并实现了其转置，加法和乘法运算，现在我们将用行主描述方法重写一个矩阵及其运算。

template<class T>
class matrix
{
    friend ostream& operator<<(ostream&,const maxtrix<T>&);
    public:
    matrix(int theRows=0,int theColumns=0);
    matrix(const matrix<T>&);
    ~matrix(){delete [] element;}
    int rows() const {return theRows;}
    int columns() const {return theColumns;}
    T& operator() (int i,int j)const;
    matrix<T>& operator=(const matrix<T>&);
    matrix<T> operator+() const;
    matrix<T> operator+(const matrix<T>&) const;
    matrix<T> operator-() const;
    matrix<T> operator-(const matrix<T>&) const;
    matrix<T> operator*(const matrix<T>&)const;
    matrix<T>& operator+=(const T&);
    private:
    int theRows,theColumns;
    T *element; 
}

以上是类定义，下面是主要实现代码

template<class T>
matrix<T>::matrix(int theRows,int theColumns)
{
    if(theRows<0||theColumns<0)
        throw illegalParameterValue("Rows and Cols must be >=0");
    if((theRows==0||theColumns==0)&&if(theRows!=0||theColumns!=0))
        throw illegalParameterValue("Either both or neither rows and columns should be zero");
    this->theRows=theRows;
    this->theColumns=theColumns;
    element=new T[theRows*theColumns];
}

template<class T>
matrix<T>::matrix(const matrix<T>& m)
{
    theRows=m.theRows;
    theColumns=m.theColumns;
    element=new T[theColumns*theRows];

    copy(m.element,m.element+rows+theColumns,element);
}

以上是构造函数和复制构造函数

template<class T>
matrix<T>& matrix<T>::operator=(const matrix<T>& m)
{
    if(this!=&m)
    {
        delete [] element;
        theRows=m.theRows;
        theColumns=m.theColumns;
        element=new T[theRows*theColumns];
        copy(m.element,m.element+theRows*theColumns,element);
    }
    return *this;
}

赋值重载符与复制构造函数类似，就是将同类复制给自己。将其他类赋值给自己需要构造函数支持。

template<class T>
T& matrix<T>::operator()(int i,int j)const
{
    if(i<1||i>theRows||j<1||j>theColumns)
        throw matrixIndexOutOfBounds();
    return element[(i-1)*theColumns+i-1];
}

获取代码有一点需要注意，那就是对矩阵而言数据是从1开始的，而不是从0
但是我们存取的时候是从零开始的所以我们要在行列上减1.

template<class T>
matrix<T> matrix<T>::operator+(const matrix<T>& m)const
{
    if(theRows!=m.theRows||theColumns!=m.columns)
        throw matrixSizeMismatch();
    matrix<T>w(theRows,theColumns);
    for(int i=0;i<theRows*theColumns;i++)
        w.element[i]=element[i]+m.element[i];
    return w;
}

这个程序很简单，我就是想强调一下，在进行任何操作之前，都要对数据
进行检查，符合标准后在写接下来的操作，比如这里就要检查数据的形状
是否相同。

template<class T>
matrix<T> matrix<T>::operator*(const matirx<T>& m)const
{
    if(theColumns!=m.theRows)
        throw matrixSizeMismatch();
    matrix<T> w(theRow,m.theColumns);

    int ct=0,cm=0,cw=0;
    for(int i=1;i<=theRows;i++)
    {
        for(int j=1;j<=m.theColumns;j++)
        {
            T sum=element[ct]*m.element[cm];//这里之所以把行首和列首首先相乘是为了下面能用++进行自增操作
            for(int k=2;k<=theColumns;k++)
            {
                ct++;
                cm+=m.theColumns;
                sum+=element[ct]*m.element[cm];
            }
            w.element[cw++]=sum;
            ct-=theColumns-1;//上一个循环执行完后ct变为thecolumns-1,这时将其归零.
                             //之所以不直接等于0是考虑下次归为第二行的0个元素。
            cm=j;//从m的第j列开始乘。
        }
        ct+=theColumns;//变为第二行
        cm=0;//m从第0列开始乘。
    }
    return w;
}

这个代码的思想和之前写的乘法思想是一模一样的，三层循环，难的地方在于映射公式的处理
比如，下一行的第一个元素要在第一行第一个元素的基础上要加theColumns等。
不过其实在我们重写了（）符号后，我们可以调用这个符号对元素进行访问，那么代码就和
之前的数组的数组来描述二维数组的矩阵乘法代码差不多了。

下面介绍几个特殊矩阵的实现
1，对角矩阵

template<class T>
class diagonalMatrix
{
    public:
    diagonalMatrix(int theN=10);
    ~diagonalMatrix(){delete [] element}
    T get(int,int)const;
    void set(int,int,const T&);
    private:
    int n;
    T*element;
}
template<class T>
diagonalMatrix<T>::diagonalMatrix(int theN)
{
    if(theN<1)
        throw illegalParameterValue("Matrix size must be >0");
    n=theN;
    element=new T[n];
}

template<class T>
T diagonalMatrix<T>::get(int i,int j)const
{
    if(i<1||j<1||i>n)
        throw matrixIndexOutOfBounds();
    if(i==j)
        return element[i-1];
    else
        return 0;   
}

template<class T>
void diagonalMatrix<T>::set(int i,int j,const T& newValue)
{
    if(i<1||j<1||i>n||j>n)
        throw matrixIndexOutOfBounds();
    if(i==j);
    element[i-1=newValue];
    else
        if(newValue!=0)
            throw illegalParameterValue("nondiagonal elements must be zero ");
}

很简单且节约内存，这就是为什么对焦矩阵在很多软件中都单独定义的原因。

1. 三对角矩阵
   要理解三对角矩阵，最关键的是理解其索引的规律，低于对角矩阵，我们的规律是i=j,对于三对角矩阵，我们有如下的规律：
   (1）主对角线——i=j;
   (2) 主对角线之下的对角线（称低对角线）——i=j+1;
   (3) 主对角线之上的对角线（称高对角线）——i=j-1;

   这三条对角线上的元素总数为3*rows-2。可以用一个容量为3*rows-2的一维数组element来描述三对角矩阵，然后就是确定元素的排列方式了，有逐行映射，逐列映射和逐条映射，逐条映射就是先从低对角线开始存储，然后对角戏，最后高对角线。这里我们采用逐条映射。
   下面只给出get方法，set方法类推：

template<class T>
T tridiagonalMatrix<T>::get(int i,int j)const
{
    if(i<1||j<1||i>n||j>n)
        throw matrixIndexOutOfBound();
    switch(i-j)
    {
        case 1:
        return element[i-2]
        case 0:
        return element[n+i-2];
        case -1:
        return element[2*n+i-2]
        default:
        return 0;
    }
}

上式中n代表矩阵的阶数（注意对角矩阵和三对角矩阵是方阵）第一条即低对角线的个数为n-1,而矩阵是从第一行第一列开始的而不是0行0列，所以对于下对角线，其第一个元素的索引为（2，1）即i从2开始，所以其对应的索引为i-2，中对角线第一个元素索引为(1,1) i从1开始，而前面有n-1个低对角线元素，所以其索引为n-1+i-1=n+i-2,同理上对角线的索引为2*n+i-2.

1. 三角矩阵
   在一个n行的下三角矩阵中，非0区域的第一行有1个元素，第二行有2个元素，……第n行有n各元素，上三角矩阵则正好与其相反。
   这两种三角矩阵都可以用一个大小为n(n+1)/2的一维数组来表示。
   同样，其排列顺序也有行主和列主之分。这里使用行主放大，且只给出set方法。

template<class T>
T lowerTridiagonalMatrix<T>::set(int i,int j,const T& newValue)
{
    if(i<1||j<1||i>n||j>n)
        throw matrixIndexOutOfBound();

    if(i>=j)
        element[i*(i-1)/2+j-1]=newValue;
    else
        if(newValue!=0)
            throw illegalParameterValue("element not in lower triangle must be zero");      
}

对于索引为(i,j)的下三角型元素，其上面有(i-1)(i-1+1)/2个非零元素，第第i行的第j个元素的索引为就j-1 因为矩阵索引从1开始计算，而我们存储时从0开始，所以加起来就是i*(i-1)/2+j-1。

4.一个n*n对称矩阵，可以视为下三角或上三角矩阵，用三角矩阵的表示方法，用一个大小为n(n+1)/2的一维数组来表示。未存储的元素可以用存储的元素来计算。