我写了一个矩阵类，以下是矩阵类的头文件，在源文件中调试时，重载矩阵乘法那一块儿出了问题，相关的部分代码如下：

Matrix operator*(Matrix &M1,Matrix &M2)
{
    Matrix M3(M1.row, M2.col);
    if (M1.row==M2.col)        
        {
            for(int i=0;i<M1.row;++i)
                for(int j=0;j<M2.col;++j)
                {
                    for(int k=0;k<M1.row;++k)
                    M3.mat[i][j]+=M1.mat[i][k]*M2.mat[k][j];
                }
        }
    return M3;
}
其中构造函数作用是产生一个全0矩阵，两个参数是其行列，不知道是哪儿除了问题，请大家 帮帮忙，谢谢！
整个头文件如下：

//Mymat.h

#include   <iostream.h>   
#include   <fstream>   
#include   <iomanip>


class   Matrix                                                              
{
protected:
    int row;   
    int col;   
    int** mat;
public:
    Matrix():row(0),col(0),mat(NULL){}
    Matrix(const Matrix& rhs);
    Matrix(const Matrix& rhs,int r,int c,int height,int width);
    Matrix(int row,int col);
    ~Matrix();
    void Create(int r,int c);
    void Clear();
    int GetRow() const {return row;}   
    int GetCol() const {return col;}  
    int Get(int i,int j) const {return mat[i][j];}
    void Set(int i,int j,int val) {mat[i][j]=val;}  
    friend istream& operator>>(istream& in,Matrix &M);
    friend ostream& operator<<(ostream& out,Matrix &M);
    friend Matrix operator+(Matrix &M1,Matrix &M2);
    Matrix& operator=(const Matrix& M);
    friend Matrix operator*(Matrix &M1,Matrix &M2);
};


void Matrix::Create(int r,int c)   
  {   
  Clear();   
  row = r;   
  col = c;   
  mat = new int*[row];   
  for(int i=0;i<row;++i)   
  {   
  mat[i] = new int[col];   
  for(int j=0;j<col;++j)   
  mat[i][j]=0;   
  }   
  } 
  
void Matrix::Clear()   
  {   
  if(mat)   
  {   
  for(int i=0;i<row;++i)   
  delete[] mat[i];   
  delete[] mat;   
  row=col=0;   
  }   
  }   

Matrix::Matrix(const Matrix& rhs)   
    {
        mat=NULL;   
        Create(rhs.GetRow(),rhs.GetCol());   
        for(int i=0; i<row;++i)   
        for(int j=0; j<col;++j)   
        mat[i][j]=rhs.Get(i,j);   
    }   

 Matrix::Matrix(const Matrix& rhs,int r,int c,int height,int width)   
  {   
  mat=NULL;   
  Create(height,width);   
  for(int i=0;i<height;++i)   
  for(int j=0;j<width;++j)   
  mat[i][j]=rhs.Get(r+i,c+j);   
  }   

Matrix::Matrix(int row,int col)   
  {   
    mat=NULL;   
    Create(row,   col);   
  }   

Matrix::~Matrix()   
  {   
  Clear();   
  } 

istream& operator>>(istream& in,Matrix &M)
{
    for(int i=0;i<M.row;++i)
    {
        for(int j=0;j<M.col;++j)
            in>>M.mat[i][j];
        }
    return in;
}

ostream& operator<<(ostream& out,Matrix &M)
{
    for(int i=0;i<M.row;++i)
    {
        for(int j=0;j<M.col;++j)
            out<<"   "<<M.mat[i][j];
        out<<endl;
    }
    return out;
}

Matrix& Matrix::operator=(const Matrix& M)   
  {   
  Clear();   
  row=M.GetRow();   
  col=M.GetCol();   
  Create(row,col);   
  for(int i=0;i<row;++i)   
  for(int j=0;j<col;++j)   
  mat[i][j]=M.Get(i,j);   
  return   *this;   
  }   

Matrix operator+(Matrix &M1,Matrix &M2)
{
    Matrix M3(M1.row,M1.col);
    if ((M1.col==M2.col)&&(M1.row==M2.row))
    {
        M3.col=M1.col;
        M3.row=M1.row;
        {
            for(int i=0;i<M1.row;++i)
                for(int j=0;j<M1.col;++j)
                    M3.mat[i][j]=M1.mat[i][j]+M2.mat[i][j];
        }
    }
    return M3;
}

Matrix operator*(Matrix &M1,Matrix &M2)
{
    Matrix M3(M1.row, M2.col);
    if (M1.row==M2.col)        
        {
            for(int i=0;i<M1.row;++i)
                for(int j=0;j<M2.col;++j)
                {
                    for(int k=0;k<M1.row;++k)
                    M3.mat[i][j]+=M1.mat[i][k]*M2.mat[k][j];
                }
        }
    return M3;
}