矩阵操作

#ifndef _SPARSE_SEQ_MATRIX
#define _SPARSE_SEQ_MATRIX

template<class T> class SparseSeqMatrix;

//稀疏序列矩阵
template<class T> class TripleType
{
public:
 friend class SparseSeqMatrix <T>;
public:
 int row, col; //非零元素所在行号/列号
 T value;  //非零元素的值
};

template<class T> class SparseSeqMatrix
{
public:
 SparseSeqMatrix()
 {
  this->Rows = 0;
  this->Cols = 0;
  this->Terms = 0;
  this->smArray = NULL;
 };
 SparseSeqMatrix(int Rows, int Cols, int Terms)
 {
  this->Rows = Rows;
  this->Cols = Cols;
  this->Terms = Terms;
  this->smArray = new TripleType<T>[Terms];
 }
 ~SparseSeqMatrix()
 {
  if (smArray != NULL)
  {
   delete[] smArray;
  }
 }

public:
 int Rows, Cols, Terms; //行/列/非零元素数
 TripleType<T> *smArray;
};

#endif

 

#ifndef _SPARSE_ORTH_MATRIX
#define _SPARSE_ORTH_MATRIX

template<class T> class OrthogonalListMatrix;

template<class T> class OLNode
{
public:
 friend class OrthogonalListMatrix <T>;
public:
 OLNode()
 {
  this->row = 0;
  this->col = 0;
  this->value = -1;
  this->right = NULL;
  this->down = NULL;
 }
 OLNode(int row, int col, T value)
 {
  this->row = row;
  this->col = col;
  this->value = value;
  this->right = NULL;
  this->down = NULL;
 }
public:
 int row,col; /*非零元素的行和列下标*/
 T value;
 OLNode<T> *right; /*非零元素所在行表、列表的后继链域*/
 OLNode<T> *down;
};

template<class T> class OrthogonalListMatrix
{
public:
 OrthogonalListMatrix()
 {
  this->Rows = 0;
  this->Cols = 0;
  this->Terms = 0;
  this->row_head = NULL;
  this->col_head = NULL;
 };
 OrthogonalListMatrix(int Rows, int Cols, int Terms)
 {
  this->Rows = Rows;
  this->Cols = Cols;
  this->Terms = Terms;
  this->row_head = new OLNode<T>*[Rows];
  this->col_head = new OLNode<T>*[Cols];
  memset(this->row_head, NULL, sizeof(OLNode<T>*)*Rows);
  memset(this->col_head, NULL, sizeof(OLNode<T>*)*Cols);
 }
 ~OrthogonalListMatrix()
 {
  for (int index=0; index<Rows; index++)
  {
   OLNode<T>*ptr = this->row_head[index];
   OLNode<T>*qtr = this->row_head[index];
   while (ptr != NULL)
   {
    qtr = ptr->right;
    delete ptr;
    ptr = qtr;
   }
  }
  delete[] this->row_head;
  delete[] this->col_head;
 }
public:
 int Rows, Cols, Terms; //行/列/非零元素数
 OLNode<T> **row_head;
 OLNode<T> **col_head;
};

#endif

 

#include <iostream>
#include <assert.h>  
#include "sparseseqmatrix.h"
#include "sparseorthmatrix.h"
using namespace std;

//特殊矩阵：非零元在矩阵中的分布有一定规则。例如: 三角矩阵，对角矩阵
//随机稀疏矩阵：非零元在矩阵中随机出现。存储方法：三元组顺序表，行逻辑联接的顺序表，十字链表。

template<class T>class CMatrix  
{  
private :
 CMatrix(){} // 必须定义, 且为private.  
 CMatrix(const CMatrix&);            // 不实现. 
 CMatrix& operator=(const CMatrix&); // 不实现.  
 ~CMatrix(){} // 可声明为public, 但这里声明为private没有错, 可被调用.

public :
 static CMatrix& GetInstance()
 {
  static CMatrix theSingleton;
  return theSingleton;
 }

public :

 //最简单的方法是把三元组表中的row与col的内容互换，然后再按照新的row中的行号对各三元组从小到大重新排序，
 //最后得到上图右半部分的三元组表，但是用这样的方法其时间复杂度为平方级的.

 //使用从列的0开始找，然后将行列互换放入新稀疏矩阵中
 //即：假设稀疏矩阵A有n列，相应地，需要针对它的三元组表中的col列进行n次扫描，
 //第k次扫描是在数组的col列中查找列号为k的三元组，若找到，
 //则意味着这个三元组所表示的矩阵元素在稀疏矩阵的第k列，需要把它存放到转置矩阵的第k行。
 void Matrix_Sparse_Transpose(T *data, int rows, int cols);
 //建立辅助数组 rowSize 和 rowStart，记录矩阵转置后各行非零元素个数和各行元素在转置三元组表中开始存放位置。
 //扫描矩阵三元组表，根据某项的列号，确定它转置后的行号，查 rowStart 表，按查到的位置直接将该项存入转置三元组表中。
 //通过设置相对位置，进行直接转置，而无需多次循环
 void Matrix_Sparse_FastTranspose(T *data, int rows, int cols);
 //稀疏矩阵相加
 void Matrix_Sparse_Add(T *data1, T *data2, int rows, int cols);
 //稀疏矩阵相减
 void Matrix_Sparse_Sub(T *data1, T *data2, int rows, int cols);
 //相乘A*B=C
 //首先，计算B中按行非0元素在B稀疏矩阵中的起始和结束为止
 //A中一行行处理
 //找到当前行中非0元素，获取列值，然后用该列值去B中对应行取值，相乘，一次找当前行处理
 void Matrix_Sparse_FastMulti(T *data1, int rows1, int cols1, T *data2, int rows2, int cols2);

 //将三对角矩阵压缩到一位数组中存储
 //以按行为主序的原则转存为一维数组M[k]中，则A[i,j]的对应关系为 k=2*i+j-2.(起始位置从1开始)
 //另一种计算方式为
 //  当i=j+1时k=3*i-3
 //  当i=j时k=3*i-2
 //  当j=i+1时k=3*i-1
 void Matrix_Tridiagonal(T *data, int rows, int cols);

 //对角线之上全是0，对应关系如下：
 //k = （i-1）*(2n-i+2)/2 + j - i;(下标从1开始)
 void Matrix_Up_Triangular(T *data, int rows, int cols);

 //对角线之下全是0，对应关系如下：
 //k = i*(i-1)/2 + j - 1;(下标从1开始)
 void Matrix_Down_Triangular(T *data, int rows, int cols);

 //将A：m*n 转换成 B：n*m
 void Matrix_Transpose(T *data, int rows, int cols);

 //对称矩阵,只保存上三角或下三角，剩下的通过对称计算推出其位置
 void Matrix_Symmetry(T *data, int rows, int cols);
private:
 //将稀疏矩阵保存到【行，列，值】为单元的数据中
 SparseSeqMatrix<T>* Matrix_SeqSparse(T *data, int rows, int cols);

 //十字链表（Orthogonal List)是有向图的另一种链式存储结构。
 //可以看成是将有向图的邻接表和逆邻接表结合起来得到的一种链表。在十字链表中，
 //对应于有向图中每一条弧都有一个结点，对应于每个定顶点也有一个结点。
 void Matrix_OrthogonalListSparse(T *data, int rows, int cols);

};

template<class T> SparseSeqMatrix<T>* CMatrix<T>::Matrix_SeqSparse(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Sparse ,data is NULL/n");
 assert(rows!=0,"In Matrix_Sparse ,rows is 0/n");
 assert(cols!=0,"In Matrix_Sparse ,cols is 0/n");

 int maxItems = 0;
 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   if (data[irow*cols + icol] != 0)
   {
    maxItems++;
   }
  }
 }
 if (maxItems == 0)
 {
  return NULL;
 }

 SparseSeqMatrix<T> *sparseMatrix = new SparseSeqMatrix<T>(rows, cols, maxItems);
 assert(sparseMatrix!=NULL,"In Matrix_Sparse ,sparseMatrix is NULL/n");

 sparseMatrix->Terms = 0;
 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   if (data[irow*cols + icol] != 0)
   {
    sparseMatrix->smArray[sparseMatrix->Terms].row = irow;
    sparseMatrix->smArray[sparseMatrix->Terms].col = icol;
    sparseMatrix->smArray[sparseMatrix->Terms].value = data[irow*cols + icol];
    sparseMatrix->Terms++;
   }
  }
 }

 cout<<sparseMatrix->Rows<<"/t"<<sparseMatrix->Cols<<"/t"<<maxItems<<endl;
 for (int index=0; index<maxItems; index++)
 {
  cout<<sparseMatrix->smArray[index].row<<"/t"
   <<sparseMatrix->smArray[index].col<<"/t"
   <<sparseMatrix->smArray[index].value<<endl;
 }
 cout<<endl;
 return sparseMatrix;
}

template<class T> void CMatrix<T>::Matrix_OrthogonalListSparse(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_OrthogonalList ,data is NULL/n");
 assert(rows!=0,"In Matrix_OrthogonalList ,rows is 0/n");
 assert(cols!=0,"In Matrix_OrthogonalList ,cols is 0/n");

 int maxItems = 0;
 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   if (data[irow*cols + icol] != 0)
   {
    maxItems++;
   }
  }
 }
 if (maxItems == 0)
 {
  return ;
 }

 OrthogonalListMatrix<T> *orthogonalListMatrix = new OrthogonalListMatrix<T>(rows, cols, maxItems);
 assert(orthogonalListMatrix!=NULL,"In Matrix_OrthogonalList ,orthogonalListMatrix is NULL/n");

 OLNode<T> *ptr = NULL;
 OLNode<T> *qtr = NULL;

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   if (data[irow*cols + icol] == 0)
   {
    continue;
   }
   ptr = new OLNode<T>(irow, icol, data[irow*cols + icol]);
   assert(ptr!=NULL,"In Matrix_OrthogonalList ,ptr is NULL/n");

   if (orthogonalListMatrix->row_head[irow] == NULL)
   {
    orthogonalListMatrix->row_head[irow] = ptr;
   }
   else
   {
    for (qtr=orthogonalListMatrix->row_head[irow]; qtr->right && qtr->right->col<icol; qtr=qtr->right);
    ptr->right = qtr->right;
    qtr->right = ptr;
   }

   if (orthogonalListMatrix->col_head[icol] == NULL)
   {
    orthogonalListMatrix->col_head[icol] = ptr;
   }
   else
   {
    for (qtr=orthogonalListMatrix->col_head[icol]; qtr->down && qtr->down->row<irow; qtr=qtr->down);
    ptr->down = qtr->down;
    qtr->down = ptr;
   }
  }
 }

 cout<<orthogonalListMatrix->Rows<<"/t"<<orthogonalListMatrix->Cols<<"/t"<<orthogonalListMatrix->Terms<<endl;
 for (int index=0; index<rows; index++)
 {
  ptr = orthogonalListMatrix->row_head[index];
  while (ptr != NULL)
  {
   cout<<ptr->row<<"/t"<<ptr->col<<"/t"<<ptr->value<<endl;
   ptr = ptr->right;
  }
 }
 cout<<endl;

 delete orthogonalListMatrix;
}

template<class T> void CMatrix<T>::Matrix_Sparse_Transpose(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Sparse ,data is NULL/n");
 assert(rows!=0,"In Matrix_Sparse ,rows is 0/n");
 assert(cols!=0,"In Matrix_Sparse ,cols is 0/n");

 SparseSeqMatrix<T> *sparseMatrix = Matrix_SeqSparse(data, rows, cols);
 assert(sparseMatrix!=NULL,"In Matrix_Sparse_Transpose ,sparseMatrix is NULL/n");

 SparseSeqMatrix<T>* transposeMatrix = new SparseSeqMatrix<T>(cols, rows, sparseMatrix->Terms);
 assert(transposeMatrix!=NULL,"In Matrix_Sparse_Transpose ,transposeMatrix is NULL/n");

 transposeMatrix->Terms = 0;
 for (int icol=0; icol<cols; icol++)
 {
  for (int iitem=0; iitem<sparseMatrix->Terms; iitem++)
  {
   if (icol == sparseMatrix->smArray[iitem].col)
   {
    transposeMatrix->smArray[transposeMatrix->Terms].col = sparseMatrix->smArray[iitem].row;
    transposeMatrix->smArray[transposeMatrix->Terms].row = sparseMatrix->smArray[iitem].col;
    transposeMatrix->smArray[transposeMatrix->Terms].value = sparseMatrix->smArray[iitem].value;
    transposeMatrix->Terms++;
   }
  }
 }

 cout<<transposeMatrix->Rows<<"/t"<<transposeMatrix->Cols<<"/t"<<transposeMatrix->Terms<<endl;
 for (int index=0; index<sparseMatrix->Terms; index++)
 {
  cout<<transposeMatrix->smArray[index].row<<"/t"
   <<transposeMatrix->smArray[index].col<<"/t"
   <<transposeMatrix->smArray[index].value<<endl;
 }
 cout<<endl;
 delete sparseMatrix;
 delete transposeMatrix;
}

template<class T> void CMatrix<T>::Matrix_Sparse_FastTranspose(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Sparse ,data is NULL/n");
 assert(rows!=0,"In Matrix_Sparse ,rows is 0/n");
 assert(cols!=0,"In Matrix_Sparse ,cols is 0/n");

 SparseSeqMatrix<T> *sparseMatrix = Matrix_SeqSparse(data, rows, cols);
 assert(sparseMatrix!=NULL,"In Matrix_Sparse_FastTranspose ,sparseMatrix is NULL/n");

 SparseSeqMatrix<T>* transposeMatrix = new SparseSeqMatrix<T>(cols, rows, sparseMatrix->Terms);
 assert(transposeMatrix!=NULL,"In Matrix_Sparse_Transpose ,transposeMatrix is NULL/n");

 int *rowSize = new int[cols];
 assert(rowSize!=NULL,"In Matrix_Sparse_FastTranspose ,rowSize is NULL/n");
 memset(rowSize, 0, sizeof(int)*(cols));

 int *rowStart = new int[cols];
 assert(rowStart!=NULL,"In Matrix_Sparse_FastTranspose ,rowStart is NULL/n");
 memset(rowStart, 0, sizeof(int)*(cols));

 for (int index=0; index<sparseMatrix->Terms; index++)
 {
  rowSize[sparseMatrix->smArray[index].col]++;
 }

 for (int index=1; index<cols; index++)
 {
  rowStart[index] = rowStart[index-1] + rowSize[index-1];
 }

 for (int index=0; index<sparseMatrix->Terms; index++)
 {
  int col = sparseMatrix->smArray[index].col;
  transposeMatrix->smArray[rowStart[col]].col = sparseMatrix->smArray[index].row;
  transposeMatrix->smArray[rowStart[col]].row = sparseMatrix->smArray[index].col;
  transposeMatrix->smArray[rowStart[col]].value = sparseMatrix->smArray[index].value;
  rowStart[col]++;
 }

 cout<<transposeMatrix->Rows<<"/t"<<transposeMatrix->Cols<<"/t"<<transposeMatrix->Terms<<endl;
 for (int index=0; index<sparseMatrix->Terms; index++)
 {
  cout<<transposeMatrix->smArray[index].row<<"/t"
   <<transposeMatrix->smArray[index].col<<"/t"
   <<transposeMatrix->smArray[index].value<<endl;
 }
 cout<<endl;

 delete sparseMatrix;
 delete transposeMatrix;
 delete[] rowSize;
 delete[] rowStart;
}

template<class T> void CMatrix<T>::Matrix_Sparse_Add(T *data1, T *data2, int rows, int cols)
{
 assert(data1!=NULL,"In Matrix_Sparse ,data1 is NULL/n");
 assert(data2!=NULL,"In Matrix_Sparse ,data2 is NULL/n");
 assert(rows!=0,"In Matrix_Sparse ,rows2 is 0/n");
 assert(cols!=0,"In Matrix_Sparse ,cols2 is 0/n");

 SparseSeqMatrix<T> *sparseMatrix1 = Matrix_SeqSparse(data1, rows, cols);
 assert(sparseMatrix1!=NULL,"In Matrix_Sparse_Add ,sparseMatrix1 is NULL/n");

 SparseSeqMatrix<T> *sparseMatrix2 = Matrix_SeqSparse(data2, rows, cols);
 assert(sparseMatrix2!=NULL,"In Matrix_Sparse_Add ,sparseMatrix2 is NULL/n");

 SparseSeqMatrix<T> *addMatrix = new SparseSeqMatrix<T>(rows, cols, sparseMatrix1->Terms + sparseMatrix2->Terms);
 assert(addMatrix!=NULL,"In Matrix_Sparse_Add ,addMatrix is NULL/n");

 int value1,value2;
 addMatrix->Terms = 0;
 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   value1 = value2 = 0;
   for (int iitem1=0; iitem1<sparseMatrix1->Terms; iitem1++)
   {
    if (irow == sparseMatrix1->smArray[iitem1].row &&
     icol == sparseMatrix1->smArray[iitem1].col)
    {
     value1 = sparseMatrix1->smArray[iitem1].value;
     break;
    }
   }

   for (int iitem2=0; iitem2<sparseMatrix2->Terms; iitem2++)
   {
    if (irow == sparseMatrix2->smArray[iitem2].row &&
     icol == sparseMatrix2->smArray[iitem2].col)
    {
     value2 = sparseMatrix2->smArray[iitem2].value;
     break;
    }
   }
   if (value1 + value2 != 0)
   {
    addMatrix->smArray[addMatrix->Terms].col = icol;
    addMatrix->smArray[addMatrix->Terms].row = irow;
    addMatrix->smArray[addMatrix->Terms].value = value1 + value2;
    addMatrix->Terms++;
   }
  }
 }

 cout<<addMatrix->Rows<<"/t"<<addMatrix->Cols<<"/t"<<addMatrix->Terms<<endl;
 for (int index=0; index<addMatrix->Terms; index++)
 {
  cout<<addMatrix->smArray[index].row<<"/t"
   <<addMatrix->smArray[index].col<<"/t"
   <<addMatrix->smArray[index].value<<endl;
 }
 cout<<endl;

 delete sparseMatrix1;
 delete sparseMatrix2;
 delete addMatrix;
}

template<class T> void CMatrix<T>::Matrix_Sparse_Sub(T *data1, T *data2, int rows, int cols)
{
 assert(data1!=NULL,"In Matrix_Sparse ,data1 is NULL/n");
 assert(data2!=NULL,"In Matrix_Sparse ,data2 is NULL/n");
 assert(rows!=0,"In Matrix_Sparse ,rows2 is 0/n");
 assert(rows!=0,"In Matrix_Sparse ,cols2 is 0/n");

 SparseSeqMatrix<T> *sparseMatrix1 = Matrix_SeqSparse(data1, rows, cols);
 assert(sparseMatrix1!=NULL,"In Matrix_Sparse_Sub ,sparseMatrix1 is NULL/n");

 SparseSeqMatrix<T> *sparseMatrix2 = Matrix_SeqSparse(data2, rows, cols);
 assert(sparseMatrix2!=NULL,"In Matrix_Sparse_Sub ,sparseMatrix2 is NULL/n");

 SparseSeqMatrix<T> *subMatrix = new SparseSeqMatrix<T>(rows, cols, sparseMatrix1->Terms + sparseMatrix2->Terms);
 assert(subMatrix!=NULL,"In Matrix_Sparse_Sub ,addMatrix is NULL/n");

 int value1,value2;
 subMatrix->Terms = 0;
 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   value1 = value2 = 0;
   for (int iitem1=0; iitem1<sparseMatrix1->Terms; iitem1++)
   {
    if (irow == sparseMatrix1->smArray[iitem1].row &&
     icol == sparseMatrix1->smArray[iitem1].col)
    {
     value1 = sparseMatrix1->smArray[iitem1].value;
     break;
    }
   }

   for (int iitem2=0; iitem2<sparseMatrix2->Terms; iitem2++)
   {
    if (irow == sparseMatrix2->smArray[iitem2].row &&
     icol == sparseMatrix2->smArray[iitem2].col)
    {
     value2 = sparseMatrix2->smArray[iitem2].value;
     break;
    }
   }
   if (value1 - value2 != 0)
   {
    subMatrix->smArray[subMatrix->Terms].col = icol;
    subMatrix->smArray[subMatrix->Terms].row = irow;
    subMatrix->smArray[subMatrix->Terms].value = value1 - value2;
    subMatrix->Terms++;
   }
  }
 }

 cout<<subMatrix->Rows<<"/t"<<subMatrix->Cols<<"/t"<<subMatrix->Terms<<endl;
 for (int index=0; index<subMatrix->Terms; index++)
 {
  cout<<subMatrix->smArray[index].row<<"/t"
   <<subMatrix->smArray[index].col<<"/t"
   <<subMatrix->smArray[index].value<<endl;
 }
 cout<<endl;

 delete sparseMatrix1;
 delete sparseMatrix2;
 delete subMatrix;
}

template<class T> void CMatrix<T>::Matrix_Sparse_FastMulti(T *data1, int rows1, int cols1, T *data2, int rows2, int cols2)
{
 assert(data1!=NULL,"In Matrix_Sparse_Multi ,data1 is NULL/n");
 assert(data2!=NULL,"In Matrix_Sparse_Multi ,data2 is NULL/n");
 assert(rows1!=0,"In Matrix_Sparse_Multi ,rows1 is 0/n");
 assert(cols1!=0,"In Matrix_Sparse_Multi ,cols1 is 0/n");
 assert(rows2!=0,"In Matrix_Sparse_Multi ,rows2 is 0/n");
 assert(cols2!=0,"In Matrix_Sparse_Multi ,cols2 is 0/n");

 assert(cols1 == rows2,"In Matrix_Sparse_Multi ,cols2 is 0/n");

 SparseSeqMatrix<T> *sparseMatrix1 = Matrix_SeqSparse(data1, rows1, cols1);
 assert(sparseMatrix1!=NULL,"In Matrix_Sparse_Multi ,sparseMatrix1 is NULL/n");

 SparseSeqMatrix<T> *sparseMatrix2 = Matrix_SeqSparse(data2, rows2, cols2);
 assert(sparseMatrix2!=NULL,"In Matrix_Sparse_Multi ,sparseMatrix2 is NULL/n");

 SparseSeqMatrix<T> *multiMatrix = new SparseSeqMatrix<T>(rows1, cols2, rows1*cols2);
 assert(multiMatrix!=NULL,"In Matrix_Sparse_Multi ,multiMatrix2 is NULL/n");
 multiMatrix->Terms = 0;

 int *rowSize = new int[rows2];
 assert(rowSize!=NULL,"In Matrix_Sparse_Multi ,rowSize is NULL/n");
 memset(rowSize, 0, sizeof(int)*(rows2));

 //记录起始点和结束点，所以最后需要多加一个
 int *rowStart = new int[rows2+1];
 assert(rowStart!=NULL,"In Matrix_Sparse_Multi ,rowSize is NULL/n");
 memset(rowStart, 0, sizeof(int)*(rows2+1));

 //保存C中当前计算出的结果行
 int *rowTmp = new int[multiMatrix->Cols];
 assert(rowTmp!=NULL,"In Matrix_Sparse_Multi ,rowTmp is NULL/n");

 //将B矩阵按行算出稀疏矩阵非0元素的相对位置
 for (int index=0; index<sparseMatrix2->Terms; index++)
 {
  rowSize[sparseMatrix2->smArray[index].row]++;
 }
 for (int index=1; index<=rows2; index++)
 {
  rowStart[index] = rowStart[index-1] + rowSize[index-1];
 }

 int itemIndex = 0;
 int curARow = 0;

 //从A矩阵行开始，从上开始，内部每次处理一行
 while (itemIndex<sparseMatrix1->Terms)
 {
  //记录A当前处理行
  curARow = sparseMatrix1->smArray[itemIndex].row;
  memset(rowTmp, 0, sizeof(int)*(multiMatrix->Cols));

  //每次A处理属于同一行的非0元素
  while (itemIndex<sparseMatrix1->Terms && curARow == sparseMatrix1->smArray[itemIndex].row)
  {
   int curACol = sparseMatrix1->smArray[itemIndex].col;//A的列对应B的行
   for(int start = rowStart[curACol]; start<rowStart[curACol+1]; start++)
   {
    //A行中N列不为0的元素和B列中所有不为0的元素相乘，放入对如列中
    int curBCol = sparseMatrix2->smArray[start].col;
    rowTmp[curBCol] = rowTmp[curBCol] + sparseMatrix1->smArray[itemIndex].value*sparseMatrix2->smArray[start].value;
   }
   itemIndex++;
  }

  //将A计算完的行放入C中
  for (int index=0; index<sparseMatrix2->Cols; index++)
  {
   if (rowTmp[index] != 0)
   {
    multiMatrix->smArray[multiMatrix->Terms].row = curARow;
    multiMatrix->smArray[multiMatrix->Terms].col = index;
    multiMatrix->smArray[multiMatrix->Terms].value = rowTmp[index];
    multiMatrix->Terms++;
   }
  }
 }

 cout<<multiMatrix->Rows<<"/t"<<multiMatrix->Cols<<"/t"<<multiMatrix->Terms<<endl;
 for (int index=0; index<multiMatrix->Terms; index++)
 {
  cout<<multiMatrix->smArray[index].row<<"/t"
   <<multiMatrix->smArray[index].col<<"/t"
   <<multiMatrix->smArray[index].value<<endl;
 }
 cout<<endl;

 delete sparseMatrix1;
 delete sparseMatrix2;
 delete multiMatrix;

 delete[] rowStart;
 delete[] rowSize;
 delete[] rowTmp;
}

template<class T> void CMatrix<T>::Matrix_Tridiagonal(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Tridiagonal ,data is NULL/n");
 assert(rows!=0,"In Matrix_Tridiagonal ,rows is 0/n");
 assert(cols!=0,"In Matrix_Tridiagonal ,cols is 0/n");

 T* tridata = new T[rows*cols];
 assert(tridata!= NULL,"In Matrix_Tridiagonal ,tridata is NULL/n");
 memset(tridata, 0, sizeof(T)*rows*cols);

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   if (irow == icol || irow == icol+1 || icol == irow+1)
   {
    tridata[2*irow+icol] = data[irow*cols + icol];
   }
  }
 }
 for (int index=0; index<rows*cols; index++)
 {
  cout<<tridata[index]<<" ";
 }
 cout<<endl;

 delete[] tridata;
}

template<class T> void CMatrix<T>::Matrix_Down_Triangular(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Down_Triangular ,data is NULL/n");
 assert(rows!=0,"In Matrix_Down_Triangular ,rows is 0/n");
 assert(cols!=0,"In Matrix_Down_Triangular ,cols is 0/n");
 assert(cols==rows,"In Matrix_Down_Triangular ,cols is 0/n");

 T* tridata = new T[rows*(rows+1)/2];
 assert(tridata!= NULL,"In Matrix_Down_Triangular ,tridata is NULL/n");
 memset(tridata, 0, sizeof(T)*(rows*(rows+1)/2));

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<=irow; icol++)
  {
   tridata[irow*(irow+1)/2 + icol] = data[irow*cols + icol];
  }
 }
 for (int index=0; index<rows*(rows+1)/2; index++)
 {
  cout<<tridata[index]<<" ";
 }
 cout<<endl;

 delete[] tridata;
}

template<class T> void CMatrix<T>::Matrix_Up_Triangular(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Up_Triangular ,data is NULL/n");
 assert(rows!=0,"In Matrix_Up_Triangular ,rows is 0/n");
 assert(cols!=0,"In Matrix_Up_Triangular ,cols is 0/n");
 assert(cols==rows,"In Matrix_Up_Triangular ,cols is 0/n");

 T* tridata = new T[rows*(rows+1)/2];
 assert(tridata!= NULL,"In Matrix_Up_Triangular ,tridata is NULL/n");
 memset(tridata, 0, sizeof(T)*(rows*(rows+1)/2));

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=irow; icol<cols; icol++)
  {
   tridata[rows*irow - irow*(irow-1)/2 + icol - irow] = data[irow*cols + icol];
  }
 }
 for (int index=0; index<rows*(rows+1)/2; index++)
 {
  cout<<tridata[index]<<" ";
 }
 cout<<endl;

 delete[] tridata;
}

template<class T> void CMatrix<T>::Matrix_Transpose(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Transpose ,data is NULL/n");
 assert(rows!=0,"In Matrix_Transpose ,rows is 0/n");
 assert(cols!=0,"In Matrix_Transpose ,cols is 0/n");

 T* tridata = new T[rows*cols];
 assert(tridata!= NULL,"In Matrix_Transpose ,tridata is NULL/n");
 memset(tridata, 0, sizeof(T)*(rows*cols));

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   cout<<data[irow*cols+icol]<<"/t";
  }
  cout<<endl;
 }
 cout<<endl;

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   tridata[icol*rows+irow] = data[irow*cols+icol];
  }
 }

 for (int irow=0; irow<cols; irow++)
 {
  for (int icol=0; icol<rows; icol++)
  {
   cout<<tridata[irow*rows+icol]<<"/t";
  }
  cout<<endl;
 }

 delete[] tridata;
}

template<class T> void CMatrix<T>::Matrix_Symmetry(T *data, int rows, int cols)
{
 assert(data!=NULL,"In Matrix_Symmetry ,data is NULL/n");
 assert(rows!=0,"In Matrix_Symmetry ,rows is 0/n");
 assert(cols!=0,"In Matrix_Symmetry ,cols is 0/n");
 assert(rows==cols,"In Matrix_Symmetry ,cols is 0/n");

 T* tridata = new T[rows*cols];
 assert(tridata!= NULL,"In Matrix_Symmetry ,tridata is NULL/n");
 memset(tridata, 0, sizeof(T)*(rows*cols));

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   cout<<data[irow*cols+icol]<<"/t";
  }
  cout<<endl;
 }
 cout<<endl;

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<=irow; icol++)
  {
   tridata[irow*(irow+1)/2 + icol] = data[irow*cols + icol];
  }
 }
 for (int index=0; index<rows*(rows+1)/2; index++)
 {
  cout<<tridata[index]<<" ";
 }
 cout<<endl;

 T* symdata = new T[rows*cols];
 assert(symdata!= NULL,"In Matrix_Symmetry ,tridata is NULL/n");
 memset(symdata, 0, sizeof(T)*(rows*cols));

 for (int irow=0; irow<rows; irow++)
 {
  for (int icol=0; icol<cols; icol++)
  {
   if (irow >= icol)
   {
    cout<<tridata[irow*(irow+1)/2 + icol]<<"/t";
   }
   else
   {
    cout<<tridata[icol*(icol+1)/2 + irow]<<"/t";
   }
  }
  cout<<endl;
 }
}