Clase Matriz con template

#include <iostream>

using namespace std;

#include <iomanip>
#include <cstdlib>
#include <cassert>

template<typename T> class Matriz;
template<typename T> ostream& operator << (ostream &salida, Matriz<T> &z);
template<typename T> istream& operator >> (istream &entrada, Matriz<T> &z);

template<typename T> Matriz<T> operator + (Matriz<T> x, Matriz<T> y);
template<typename T> Matriz<T> operator - (Matriz<T> x, Matriz<T> y);
template<typename T> Matriz<T> operator * (Matriz<T> x, Matriz<T> y);
template<typename T> Matriz<T> operator * (T x, Matriz<T> y);

template <class T>

class Matriz {

private:

    T **a;
    int rows, cols;

public:

    Matriz(int m, int n);
    ~Matriz();  // Destructor de la clase

    void initMatriz(T **x);
    void inverMat(T **a);

    T deterMat();

    friend ostream& operator << <> (ostream &salida, Matriz<T> &z); //sobrecarga del operador <<
    friend istream& operator >> <> (istream &entrada, Matriz<T> &z); //sobrecarga del operador >>

    friend Matriz   operator  + <> (Matriz<T> x, Matriz<T> y); //sobrecarga al operador + para sumar matrices
    friend Matriz   operator  - <> (Matriz<T> x, Matriz<T> y); //sobrecarga al operador - para sumar matrices
    friend Matriz   operator  * <> (Matriz<T> x, Matriz<T> y); //sobrecarga al operador * para multiplicar matrices
    friend Matriz   operator  * <> (T x, Matriz<T> y); //sobrecarga al operador * para multiplicar escalar*matriz


};

//Constructor por parametros

template<typename T>

Matriz<T>::Matriz(int m , int n){

    int i;

    rows = m;

    cols = n;

    a = new T*[rows];

    for(i = 0; i< rows; i++){

        a[i] = new T[cols];

    }

}

template <typename T>

Matriz<T>::~Matriz() {}

template<typename T>

void Matriz<T>::initMatriz(T **x){

    int i,j;

    for(i=0; i < rows; i++){

            for(j=0; j < cols; j++){

                a[i][j] = x[i][j];

        }
    }


}

template<typename T>

istream &operator >> (istream &entrada, Matriz<T> &z){//sobrecarga del operador >>

    cout << "\nIntroduzca los valores a traves del teclado:\n\n";

    int i, j;

    for (i = 0; i < z.rows; i++) {

        for (j = 0; j < z.cols; j++) {

            cout << "A[" << i+1 <<"][" << j+1 << "] = ? ";

            entrada >> z.a[i][j];

        }

    }

    return entrada;

}

template<typename T>

ostream &operator << (ostream &salida, Matriz<T> &z){//sobrecarga del operador <<

    salida.setf(ios::fixed);
    salida.precision(4);

    int i, j;

    for (i = 0; i < z.rows; i++) {

        for (j = 0; j < z.cols; j++) {

            salida << setw(10) << z.a[i][j];

        }

        salida << "\n";

    }

    return salida;

}

template<typename T>

Matriz<T> operator +(Matriz<T> x, Matriz<T> y){//sobrecarga al operador + para sumar matrices

    assert(x.rows == y.rows || x.cols == y.cols);

    Matriz<T> sum_mat(x.rows, x.cols);

    int i,j;

    for(i=0; i < x.rows; i++){

        for(j=0; j < x.cols; j++){

            sum_mat.a[i][j] = x.a[i][j] + y.a[i][j];

        }
    }

    return sum_mat;

}

template<typename T>

Matriz<T> operator -(Matriz<T> x, Matriz<T> y){//sobrecarga al operador + para restar matrices

    assert(x.rows == y.rows || x.cols == y.cols);

    Matriz<T> rest_mat(x.rows, x.cols);

    int i,j;

    for(i=0; i < x.rows; i++){

        for(j=0; j < x.cols; j++){

            rest_mat.a[i][j] = x.a[i][j] + y.a[i][j];

        }
    }

    return rest_mat;

}

template<typename T>

Matriz<T> operator *(Matriz<T> x, Matriz<T> y){//sobrecarga al operador * para multiplicar matrices

    assert(x.cols == y.rows);

    Matriz<T> mult_mat(x.rows, y.cols);

    int i,j,k;

    for(i=0; i < x.rows; i++){

        for(j=0; j < x.cols; j++){

            mult_mat.a[i][j] = 0;

            for(k=0; k < x.cols; k++){

                mult_mat.a[i][j] += x.a[i][k] * y.a[k][j];

            }

        }
    }

    return mult_mat;

}

template<typename T>

Matriz<T> operator *(T x, Matriz<T> y){//sobrecarga al operador * para multiplicar escalar*matriz

    Matriz<T> esc_mat(y.rows, y.cols);

    int i,j;

    for(i=0; i < y.rows; i++){

        for(j=0; j < y.cols; j++){

            esc_mat.a[i][j] = x*y.a[i][j];

        }
    }

    return esc_mat;

}

template<typename T>

void Matriz<T>::inverMat(T **z) {

// Algoritmo para la eliminación simple de Gauss

    int i, j, k, n= rows, m = cols;

    assert (rows == cols);

    double factor;

    double **L, **ainv, *D, *X;

    X = new double [n]; D = new double [n];

    L = new double* [n];    ainv = new double* [n];

    for (j = 0; j < n; j++) {

        L[j] = new double [n]; ainv[j] = new double [n];

    }

    for (k = 0; k < n - 1; k++) {

        for (i = k+1; i < n;  i++) {

            factor = z[i][k]/z[k][k];

            for (j = k+1; j < n + 1; j++) {

                z[i][j] = z[i][j] - factor * z[k][j];

            }



        }

    }

// Cálculo del determinante

        double determ;

        determ = 1.;

    for (i = 0; i < n; i++) {

        determ = determ * z[i][i];

    }

    assert(determ != 0);

// Rutina para determinar las matrices L (inferior) y U (superior) de la
// descomposición LU


        for (i = 0; i < n; i++) {

               for (j = 0; j < n; j++) {

                  if (i > j) {

                      L[i][j] = z[i][j]/z[j][j];

                      z[i][j] = 0;

                   }


               }

        }


       for (i = 0; i < n; i++) {

            for (j = 0; j < n; j++) {

                  L[j][j] = 1;

            }

        }


// Implementación de la rutina para el cálculo de la inversa


 for (k = 0; k < n; k++) {


// Esta rutina inicializa los L[i][n] para ser utilizados con la matriz L


          for (i = 0; i < n; i++) {

               if (i == k) L[i][n] = 1;

               else  L[i][n] = 0;

           }

// Esta función implementa la sustitución hacia adelante con la matriz L y los L[i][n]
// que produce la rutina anterior

  double sum;

  D[0] = L[0][n];

  for (i = 1; i < n; i++) {

       sum = 0;

       for (j = 0; j < i; j++) {

            sum = sum + L[i][j]*D[j];

       }

        D[i] = L[i][n] - sum;

   }



// Esta rutina asigna los D[i] que produce forward para ser utilizados con la matriz U

  for (i = 0; i < n; i++) {

          z[i][n] = D[i];

  }

// Rutina que aplica la sustitución hacia atras


 X[n-1] = z[n-1][n]/z[n-1][n-1];

 // Determinación de las raíces restantes


  for (i = n - 2; i > -1; i--) {

        sum = 0;

        for (j = i+1; j < n; j++) {

              sum = sum + z[i][j]*X[j];

         }

         X[i] = (z[i][n] - sum)/z[i][i];

   }


// Esta rutina asigna los X[i] que produce Sustituir como los elementos de la matriz inversa

  for (i = 0; i < n; i++) {

         ainv[i][k] = X[i];

  }

 }   // llave de cierre del for para k

    for(i=0; i < n; i++){

            for(j=0; j < n; j++){

                a[i][j] = ainv[i][j];

        }
    }

    delete L, D, X, ainv;

}

template<typename T>

T Matriz<T>::deterMat() {

// Algoritmo para la eliminación simple de Gauss

    int i, j, k, n=rows, m=cols;

    assert (rows == cols);

    double **b;

    b = new double* [n];

    for(j=0; j<n; j++){

        b[j] = new double [n];

    }

    for(i=0; i<n; i++){

        for(j=0; j<n; j++){

            b[i][j] = a[i][j];

        }

    }

    double factor;

    for (k = 0; k < n - 1; k++) {

        for (i = k+1; i < n;  i++) {

            factor = b[i][k]/b[k][k];

            for (j = k+1; j < n + 1; j++) {

                b[i][j] = b[i][j] - factor * b[k][j];

            }



        }

    }

// Cálculo del determinante

    double determ = 1.;

    for (i = 0; i < n; i++) {

        determ = determ * b[i][i];

    }


    delete [] b;

    return determ;

}