Ошибка сегментации, перегрузка ostream (‹‹)

Итак, я практикую кодирование на С++ и пытаюсь написать класс для матриц (хранящихся в виде массивов) с соответствующими перегруженными операциями.

Мне удалось определить класс и попытаться перегрузить оператор ‹‹, но мой текущий код вызывает ошибку сегментации (я компилирую с помощью g++ в Ubuntu). Я посмотрел в Интернете, и люди с похожими проблемами, как правило, всегда забывают возвращать os в своей функции перегрузки, но я сделал это, поэтому понятия не имею, в чем может быть моя проблема. Кроме того, мой перегруженный оператор сработает несколько раз, прежде чем вызовет ошибку сегментации.

Любая помощь будет принята с благодарностью.

Вот мой код:

#include<iostream>
#include<stdlib.h> // for c style exit
using namespace std;

class matrix
{
  // Friends
  friend ostream & operator<<(ostream &os, const matrix &mat);
  friend istream & operator>>(istream &is, matrix &mat);

private:
  double *mdata;
  int rows,columns;
public:
  // Default constructor
  matrix(){mdata=0; rows=columns=0;}
  // Parameterized constructor
  matrix(int m, int n){mdata = new double[ m*n ]; rows = m; columns = n;}
  // Copy constructor
  matrix(matrix &mat)
  // Destructor
  ~matrix(){delete[] mdata; cout<<"Destructing array."<<endl;}
  // Access functions
  int getrows() const {return rows;} // Return number of rows
  int getcols() const {return columns;} // Return number of columns
  int index(int m, int n) const // Return position in array of element (m,n)
  {
    if(m>0 && m<=rows && n>0 && n<=columns) return (n-1)+(m-1)*columns;
    else {cout<<"Error: out of range"<<endl; exit(1);}
  }
  double & operator()(int m, int n)const {return mdata[index(m,n)];}
  // Other access functions go here
  double & operator[](int i) {return mdata[i];}
  // Other functions 
  // Copy  Assignment operator
  matrix & operator=(matrix &mat);
};

// Member functions defined outside class
matrix::matrix(matrix &mat){
  rows = mat.getrows();
  columns = mat.getcols();
  for(int j = 0; j<rows*columns; j++){mdata[j] = mat[j];}
 }

matrix & matrix::operator=(matrix &mat){
  if (&mat == this) return *this;

  delete[] mdata; rows = 0; columns = 0;

  rows = mat.getrows(); columns = mat.getcols();
  if(rows>0&&columns>0){
    mdata = new double[(columns-1) + (rows-1)*columns + 1];
    for(int j = 0; j<rows*columns; j++){mdata[j] = mat[j];}
  }
  return *this;
}


// Overload insertion to output stream for matrices
ostream & operator<<(ostream &os, const matrix &mat){
  for(int j = 0;j<mat.rows;j++){
    for(int k = 0;k<mat.columns;k++){
      os << mat(j+1,k+1) << " ";
   }
    os << endl;
 }
  return os;
}

// Main program

int main(){

  // Demonstrate default constructor
  matrix a1;
  cout<<a1;

  // Parameterized constructor
  const int m(2),n(2);
  matrix a2(m,n);
  // Set values for a2 here
  a2[0] = 1; a2[1] = 2; a2[2] = 3; a2[3] = 4;
  // Print matrix a2
  cout<<a2;


  // Deep copy by assignment: define new matrix a3 then copy from a2 to a3
  matrix a3(m,n);
  cout<<a3;
  a3=a2;
  cout<<a3;
  // Modify contents of original matrix and show assigned matrix is unchanged here
  a2[0] = 5;
  cout<<a2;
  cout<<a3; //here is where segmentation fault occurs
  return 0;
}