I have Gentoo linux 1.4 (not sure about rc), running on a Pentium 200 MMX, built from stage 1, with GCC 3.2.2 (which has built my entire system, including Gnome, and everything works fine). I tried compiling my own C++ program with "g++ Matrix.cpp -o matrixtest", and got the following: Matrix.h: In static member function `static Matrix<T> Matrix<T>::Householder(Matrix<T>) [with T = float]': Matrix.cpp:702: instantiated from here Matrix.h:264: Internal compiler error in c_expand_expr, at c-common.c:3714 Please submit a full bug report, with preprocessed source if appropriate. See <URL:http://bugs.gentoo.org/> for instructions. I eventually managed to get the program compiled by moving the declaration of the class "SizeMismatch" further up in the program (almost to the beginning), and changing "throw SizeMismatch;" to "throw SizeMismatch();". But still, the compiler should tell you "SizeMismatch undeclared", not "internal error"... Matrix.h and Matrix.cpp follow. Matrix.h contents: (cut here)------------------------------------------------------ typedef unsigned long u32; typedef signed long s32; // Forward declaration template<typename T> class Matrix; // Expression template classes: template<typename T, typename Left, typename Right> class MatrixAddPair { const Left &left; const Right &right; public: MatrixAddPair(const Left &l, const Right &r) : left(l), right(r) {} T operator()(u32 i, u32 j) const { return left(i,j) + right(i,j); } u32 Columns() const {return left.Columns();} u32 Rows() const {return left.Rows();} }; template<typename T, typename Left, typename Right> class MatrixSubPair { const Left &left; const Right &right; public: MatrixSubPair(const Left &l, const Right &r) : left(l), right(r) {} T operator()(u32 i, u32 j) const { return left(i,j) - right(i,j); } u32 Columns() const {return left.Columns();} u32 Rows() const {return left.Rows();} }; template<typename T, typename Left, typename Right> class MatrixMulPair { const Left &left; const Right &right; public: MatrixMulPair(const Left &l, const Right &r) : left(l), right(r) {} T operator()(u32 i, u32 j) const { T sum = 0; const u32 rows = right.Rows(); for (u32 k = 0; k < rows; k++) sum += left(i,k)*right(k,j); return sum; } u32 Columns() const {return right.Columns();} u32 Rows() const {return left.Rows();} }; template<typename T, typename M, typename Unused=int> class MatrixMulScalarPair { const M &matrix; const T scalar; public: MatrixMulScalarPair(const M &m, T t) : matrix(m), scalar(t) {} T operator()(u32 i, u32 j) const { return matrix(i,j)*scalar; } u32 Columns() const {return matrix.Columns();} u32 Rows() const {return matrix.Rows();} }; // Operators //Addition template<typename T, typename L1, typename L2, typename R1, typename R2, template<typename,typename,typename> class Left, template<typename,typename,typename> class Right> MatrixAddPair<T, Left<T,L1,L2>, Right<T,R1,R2> > inline operator+(const Left<T,L1,L2> &l, const Right<T,R1,R2> &r) { return MatrixAddPair<T, Left<T,L1,L2>, Right<T,R1,R2> > (l,r); } template<typename T, typename L1, typename L2, template<typename,typename,typename> class Left> MatrixAddPair<T, Left<T,L1,L2>, Matrix<T> > inline operator+(const Left<T,L1,L2> &l, const Matrix<T> &r) { return MatrixAddPair<T, Left<T,L1,L2>, Matrix<T> > (l,r); } template<typename T, typename R1, typename R2, template<typename,typename,typename> class Right> MatrixAddPair<T, Matrix<T>, Right<T,R1,R2> > inline operator+(const Matrix<T> &l, const Right<T,R1,R2> &r) { return MatrixAddPair<T, Matrix<T>, Right<T,R1,R2> > (l,r); } template<typename T> MatrixAddPair<T, Matrix<T>, Matrix<T> > inline operator+(const Matrix<T> &l, const Matrix<T> &r) { return MatrixAddPair<T, Matrix<T>, Matrix<T> > (l,r); } // Subtraction template<typename T, typename L1, typename L2, typename R1, typename R2, template<typename,typename,typename> class Left, template<typename,typename,typename> class Right> MatrixSubPair<T, Left<T,L1,L2>, Right<T,R1,R2> > inline operator-(const Left<T,L1,L2> &l, const Right<T,R1,R2> &r) { return MatrixSubPair<T, Left<T,L1,L2>, Right<T,R1,R2> > (l,r); } template<typename T, typename L1, typename L2, template<typename,typename,typename> class Left> MatrixSubPair<T, Left<T,L1,L2>, Matrix<T> > inline operator-(const Left<T,L1,L2> &l, const Matrix<T> &r) { return MatrixSubPair<T, Left<T,L1,L2>, Matrix<T> > (l,r); } template<typename T, typename R1, typename R2, template<typename,typename,typename> class Right> MatrixSubPair<T, Matrix<T>, Right<T,R1,R2> > inline operator-(const Matrix<T> &l, const Right<T,R1,R2> &r) { return MatrixSubPair<T, Matrix<T>, Right<T,R1,R2> > (l,r); } template<typename T> MatrixSubPair<T, Matrix<T>, Matrix<T> > inline operator-(const Matrix<T> &l, const Matrix<T> &r) { return MatrixSubPair<T, Matrix<T>, Matrix<T> > (l,r); } // Multiplication template<typename T, typename L1, typename L2, typename R1, typename R2, template<typename,typename,typename> class Left, template<typename,typename,typename> class Right> MatrixMulPair<T, Left<T,L1,L2>, Right<T,R1,R2> > inline operator*(const Left<T,L1,L2> &l, const Right<T,R1,R2> &r) { return MatrixMulPair<T, Left<T,L1,L2>, Right<T,R1,R2> > (l,r); } template<typename T, typename L1, typename L2, template<typename,typename,typename> class Left> MatrixMulPair<T, Left<T,L1,L2>, Matrix<T> > inline operator*(const Left<T,L1,L2> &l, const Matrix<T> &r) { return MatrixMulPair<T, Left<T,L1,L2>, Matrix<T> > (l,r); } template<typename T, typename R1, typename R2, template<typename,typename,typename> class Right> MatrixMulPair<T, Matrix<T>, Right<T,R1,R2> > inline operator*(const Matrix<T> &l, const Right<T,R1,R2> &r) { return MatrixMulPair<T, Matrix<T>, Right<T,R1,R2> > (l,r); } template<typename T> MatrixMulPair<T, Matrix<T>, Matrix<T> > inline operator*(const Matrix<T> &l, const Matrix<T> &r) { return MatrixMulPair<T, Matrix<T>, Matrix<T> > (l,r); } // Matrix * scalar template<typename T, typename L1, typename L2, template<typename,typename,typename> class Left> MatrixMulScalarPair<T, Left<T,L1,L2> > inline operator*(const Left<T,L1,L2> &l, T r) { return MatrixMulScalarPair<T, Left<T,L1,L2> > (l,r); } template<typename T> MatrixMulScalarPair<T, Matrix<T> > inline operator*(const Matrix<T> &l, T r) { return MatrixMulScalarPair<T, Matrix<T> > (l,r); } // Scalar * matrix template<typename T, typename R1, typename R2, template<typename,typename,typename> class Right> MatrixMulScalarPair<T, Right<T,R1,R2> > inline operator*(T r, const Right<T,R2,R2> &l) { return MatrixMulScalarPair<T, Right<T,R1,R2> > (l,r); } template<typename T> MatrixMulScalarPair<T, Matrix<T> > inline operator*(T r, const Matrix<T> &l) { return MatrixMulScalarPair<T, Matrix<T> > (l,r); } // The matrix class template<typename T> class Matrix { u32 rows; u32 columns; T **entry; public: // Constructors Matrix(const Matrix&); Matrix(u32 newRows, u32 newColumns); // Comparison and assignment template<typename RHS> Matrix& operator=(const RHS&); template<typename RHS> bool operator==(const RHS&) const; // Mutating operators on this matrix Matrix& operator+=(const Matrix&); Matrix& operator-=(const Matrix&); Matrix& operator*=(const T); // Inlines T operator()(u32 i, u32 j) const { return entry[i][j]; } T& operator()(u32 i, u32 j) { return entry[i][j]; } u32 Columns() const {return columns;} u32 Rows() const {return rows;} // Static methods // Householder matrix for column vector x static Matrix Householder(Matrix x) { if (x.Columns() != 1) throw SizeMismatch; T alpha = sqrt((x.Transposed() * x)(0,0)); if (x(0,0) >= 0) alpha = -alpha; x(0,0) -= alpha; x *= 1/sqrt((x.Transposed() * x)(0,0)); Matrix M(x.Columns(), x.Columns()); M.LoadIdentity(); M = M - 2.0f*x*x.Transposed(); return M; } // Getters: // Adjoint Matrix Adjoint() const; // Determinant (by Gaussian elimination) T Det() const; // Determinant (by Laplace expansion - slow) T Det2() const; // Displays using C++ streams void Display() const; // Is it diagonally dominant? bool IsStrictlyDiagonallyDominant() const; // l1 norm T L1Norm() const; // l infinity norm T LInfNorm() const; // Submatrix without a row and column Matrix SubmatrixWithout(u32 row, u32 col) const; // Returns the transpose Matrix Transposed() const; // Setters // Divides row by scalar void DivideRow(u32, T); // Inverts using Gauss-Jordan reduction with // partial pivoting bool Invert(); // Inverts using 1/det(a)*Adj(A) bool Invert2(); // Makes matrix into I void LoadIdentity(); // Fills with 0's void LoadZeroes(); // Swaps 2 rows void SwapRows(u32, u32); // Row subtraction void SubRows(u32 source, T multiple, u32 dest, u32 start, u32 end); // Tranposes the matrix void Transpose(); ~Matrix(); }; class SizeMismatch{}; // Functions // Gauss-Seidel solution of STRICTLY // DIAGONALLY DOMINANT matrix, with error less // than err, in at most maxIterations iterations template<typename T> void GaussSeidel(const Matrix<T> &A, const Matrix<T> &b, Matrix<T> &x, T err, u32 maxIterations = 100); // Gives the in-place LU decomposition of A // with rowSwap having at least as many elements as // A's rows. Returns true on success, false otherwise. template<typename T> bool LUDecompose(Matrix<T> &A, u32 rowSwap[]); // Solves for the column vector b, in place, using // the LU decomposition of A given by LUDecompose() template<typename T> void LUSolve(const Matrix<T> &A, Matrix<T> &b, const u32 rowSwap[]); // Positive integer power of a matrix template<typename T> Matrix<T> pow(const Matrix<T> &A, u32 exponent); // Solution of Ax = b by Gaussian elimination. // Matrix A may be any size that matches b. // A and b are always changed. // If unique solution, Solve() return 1, // b is changed to x. // If no solution, Solve() return 0. // If infinite solutions, basis contains vectors // of the form form // [x1 x2 x3 x4 ... xn], so that // x = x1 + x2*s + x3*t + ... + xn*u, // x's column vectors; s,t...u parameters; // n is returned. template<typename T> u32 Solve(Matrix<T> &A, Matrix<T> &b, Matrix<T> *basis = 0); (end of Matrix.h)------------------------------------------------------- Matrix.cpp (cut here)-------------------------------------------------------------- #include <iostream> #include "Matrix.h" #define abs(x) ((x) >= 0 ? (x) : (-(x))) template <typename T> Matrix<T>::Matrix(u32 newRows, u32 newColumns) { if ((newRows == 0) || (newColumns == 0)) throw SizeMismatch(); rows = newRows; columns = newColumns; entry = new T*[rows]; for (u32 i = 0; i < rows; i++) { entry[i] = new T[columns]; for (u32 j = 0; j < columns; j++) entry[i][j] = 0; } } template <class T> Matrix<T>::Matrix(const Matrix<T>& m) { rows = m.rows; columns = m.columns; entry = new T*[rows]; for (u32 i = 0; i < rows; i++) { entry[i] = new T[columns]; for (u32 j = 0; j < columns; j++) entry[i][j] = m.entry[i][j]; } } template <class T> template <typename RHS> Matrix<T>& Matrix<T>::operator=(const RHS& rhs) { if (rhs.Rows() != rows) { for (u32 i = 0; i < rows; i++) delete[] entry[i]; delete[] entry; columns = rhs.Columns(); rows = rhs.Rows(); entry = new T*[rows]; for (u32 i = 0; i < rows; i++) entry[i] = new T[columns]; } if (rhs.Columns() != columns) { columns = rhs.Columns(); for (u32 i = 0; i < rows; i++) { delete[] entry[i]; entry[i] = new T[columns]; } } for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) entry[i][j] = rhs(i,j); return(*this); } template <class T> template <typename RHS> bool Matrix<T>::operator==(const RHS& rhs) const { if ((rhs.Rows() != rows) || (rhs.Columns() != columns)) return false; for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) if (entry[i][j] != rhs(i,j)) return false; return true; } template <class T> Matrix<T>& Matrix<T>::operator+=(const Matrix<T>& m) { if ((m.Rows != rows) || (m.columns != columns)) throw SizeMismatch(); for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) entry[i][j] += m.entry[i][j]; return(*this); } template <class T> Matrix<T>& Matrix<T>::operator-=(const Matrix<T>& m) { if ((m.rows != rows) || (m.columns != columns)) throw SizeMismatch(); for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) entry[i][j] -= m.entry[i][j]; return(*this); } template <class T> Matrix<T>& Matrix<T>::operator*=(const T scalar) { for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) entry[i][j] *= scalar; return(*this); } template<class T> Matrix<T> Matrix<T>::Adjoint() const { if (rows != columns) throw SizeMismatch(); Matrix<T> Adj(rows, columns); T sign; for (u32 i = 0; i < rows; i++) { if (i & 1) sign = -1; else sign = 1; for (u32 j = 0; j < columns; j++) { Adj(j,i) = sign*SubmatrixWithout(i,j).Det(); sign = -sign; } } return Adj; } template<class T> T Matrix<T>::Det() const { if (rows != columns) throw SizeMismatch(); Matrix<T> A = *this; u32 i,j; T pivot, determinant = 1; u32 pivotRow; // Gaussian elimination into upper triangular form for (j = 0; j < A.Columns(); j++) { pivot = abs(A(j,j)); pivotRow = j; // Find the biggest partial pivot for (i = j+1; i < A.Rows(); i++) if (abs(A(i,j)) > pivot) { pivot = abs(A(i,j)); pivotRow = i; } // Is the matrix singular? if (pivot == 0) return false; else { if (pivotRow != j) { A.SwapRows(j, pivotRow); determinant = -determinant; } // Reduce entries below pivot to 0 for (i = j+1; i < A.Rows(); i++) { pivot = A(i,j) / A(j,j); A.SubRows(i, pivot, j, j+1, A.Columns()); A(i,j) = pivot; } } } for (i = 0; i < A.Rows(); i++) determinant *= A(i,i); return determinant; } template <class T> T Matrix<T>::Det2() const { if (rows != columns) throw SizeMismatch(); if (rows == 1) return entry[0][0]; else if (rows == 2) return entry[0][0]*entry[1][1] - entry[0][1]*entry[1][0]; else { T sign = 1; T answer = 0; for (u32 i = 0; i < rows; i++) { answer += sign*entry[0][i]*SubmatrixWithout(0,i).Det(); sign = -sign; } return answer; } } template <class T> void Matrix<T>::Display() const { using namespace std; for (u32 i = 0; i < rows; i++) { for (u32 j = 0; j < columns; j++) cout << entry[i][j] << " "; cout << endl; } } template <class T> void Matrix<T>::DivideRow(u32 i, T value) { for (u32 k = 0; k < columns; k++) entry[i][k] /= value; } template<class T> void GaussSeidel(const Matrix<T> &A, const Matrix<T> &b, Matrix<T> *q, T err, u32 maxIterations) { Matrix<T> &x = *q; u32 i,j; T sum, error, maxError; u32 iteration = 0; do { maxError = 0; for (i = 0; i < A.Rows(); i++) { sum = 0; for (j = 0; j < i; j++) sum += x(j,0)*A(i,j); for (j++; j < A.Columns(); j++) sum += x(j,0)*A(i,j); error = x(i,0); x(i,0) = (b(i,0)-sum) / A(i,i); error = abs(x(i,0) - error); if (error > maxError) maxError = error; } iteration++; } while ((iteration < maxIterations) && (maxError > err)); } template <class T> bool Matrix<T>::Invert() { if (rows != columns) throw SizeMismatch(); s32 i, j, pivotRow; T pivot; // Keep track of row swaps u32 rowPosition[rows]; // Gaussian elimination into upper triangular form for (j = 0; j < columns; j++) { pivot = abs(entry[j][j]); pivotRow = j; // Find the biggest partial pivot for (i = j+1; i < rows; i++) if (abs(entry[i][j]) > pivot) { pivot = abs(entry[i][j]); pivotRow = i; } // Is the matrix singular? if (pivot == 0) return false; else { if (pivotRow != j) SwapRows(j, pivotRow); rowPosition[j] = pivotRow; // Divide row pivot = entry[j][j]; entry[j][j] = 1.0; DivideRow(j, pivot); // Reduce entries below pivot to 0 for (i = j+1; i < rows; i++) { pivot = entry[i][j]; entry[i][j] = 0.0; SubRows(i, pivot, j, 0, columns); } // Reduce entries above pivot to 0 for (i = j-1; i >= 0; i--) { pivot = entry[i][j]; entry[i][j] = 0.0; SubRows(i, pivot, j, 0, columns); } } } // Fix columns in the inverse for (j = columns-1; j >= 0; j--) for (i = rows-1; i >= 0; i--) if (j != rowPosition[j]) Swap(entry[i][j], entry[i][rowPosition[j]]); return true; }; template <class T> bool Matrix<T>::Invert2() { T determinant = Det(); if (determinant == 0) return false; else { Matrix<T> A = Adjoint(); *this = A*(1.0/determinant); return true; } } template <class T> bool Matrix<T>::IsStrictlyDiagonallyDominant() const { u32 i,j; T sum; for (i = 0; i < rows; i++) { sum = 0; for (j = 0; j < i; j++) sum += abs(entry[i][j]); for (j++; j < columns; j++) sum += abs(entry[i][j]); if (abs(entry[i][i]) <= sum) return false; } return true; } template <class T> T Matrix<T>::L1Norm() const { T *sums = new T[columns]; for (u32 i = 0; i < columns; i++) sums[i] = 0; for (u32 i = 0; i < columns; i++) for (u32 j = 0; j < rows; j++) sums[i] += abs(entry[j][i]); T max = sums[0]; for (u32 i = 0; i < columns; i++) if (sums[i] > max) max = Sums[i]; delete[] sums; return max; } template <class T> T Matrix<T>::LInfNorm() const { T *sums = new T[rows]; for (u32 i = 0; i < rows; i++) sums[i] = 0; for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) sums[i] += abs(entry[i][j]); T max = sums[0]; for (u32 i = 0; i < rows; i++) if (sums[i] > max) max = sums[i]; delete[] sums; return max; } template <class T> void Matrix<T>::LoadIdentity() { if (rows == columns) { for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) if (i == j) entry[i][j] = 1; else entry[i][j] = 0; } else throw SizeMismatch(); } template <class T> void Matrix<T>::LoadZeroes() { for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) entry[i][j] = 0; } template<class T> bool LUDecompose(Matrix<T> &A, u32 rowSwap[]) { if (A.Rows() != A.Columns()) throw SizeMismatch(); s32 i, j, pivotRow; T pivot; for (i = 0; i < A.Rows(); i++) rowSwap[i] = i; // Gaussian elimination into upper triangular form for (j = 0; j < A.Columns(); j++) { pivot = abs(A(j,j)); pivotRow = j; // Find the biggest partial pivot for (i = j+1; i < A.Rows(); i++) if (abs(A(i,j)) > pivot) { pivot = abs(A(i,j)); pivotRow = i; } // Is the matrix singular? if (pivot == 0) return false; else { if (pivotRow != j) A.SwapRows(j, pivotRow); rowSwap[j] = pivotRow; // Reduce entries below pivot to 0 for (i = j+1; i < A.Rows(); i++) { pivot = A(i,j) / A(j,j); A.SubRows(i, pivot, j, j+1, A.Columns()); A(i,j) = pivot; } } } return true; } template<class T> void LUSolve(const Matrix<T> &A, Matrix<T> &b, const u32 * const rowSwap) { if ((A.Columns() != b.Rows()) || (b.Columns() != 1)) throw SizeMismatch(); s32 i,j; T sum; // Swap rows for (i = 0; i < b.Rows(); i++) b.SwapRows(i, rowSwap[i]); // Solve Ly = b for (i = 0; i < b.Rows(); i++) { sum = 0; for (j = 0; j < i; j++) sum += b(j,0)*A(i,j); b(i,0) -= sum; } // Solve Ux = y for (i--; i >= 0; i--) { sum = 0; for (j = i+1; j < A.Columns(); j++) sum += A(i,j)*b(j,0); b(i,0) = (b(i,0)-sum) / A(i,i); } } template<class T> Matrix<T> pow(const Matrix<T> &A, u32 exponent) { if (A.Rows() != A.Columns()) throw SizeMismatch(); if (exponent == 0) { Matrix<T> I(A.Rows(), A.Rows()); I.LoadIdentity(); return I; } Matrix<T> B = A; while (exponent != 1) { B = B*B; if ((exponent & 1) != 0) B = B*A; exponent >>= 1; } return B; } template<class T> u32 Solve(Matrix<T> &A, Matrix<T> &b, Matrix<T> *basis) { if ((A.Rows() != b.Rows()) || (b.Columns() != 1)) throw SizeMismatch(); // Current row and column s32 i = 0, j = 0; // Counter s32 k; T pivot; s32 pivotRow; // Gaussian elimination into upper triangular form do { pivot = abs(A(i,j)); pivotRow = i; // Find the biggest partial pivot, if any for (k = i+1; k < A.Rows(); k++) if (abs(A(k,j)) > pivot) { pivot = abs(A(k,j)); pivotRow = k; } // Is there a pivot? if (pivot == 0) j++; else { A.SwapRows(i, pivotRow); b.SwapRows(i, pivotRow); // Reduce entries below pivot to 0 for (k = i+1; k < A.Rows(); k++) { b(k,0) -= A(k,j)/A(i,j)*b(i,0); A.SubRows(k, A(k,j)/A(i,j), i, j+1, A.Columns()); A(k,j) = 0; } i++; j++; } } while ((j < A.Columns()) && (i < A.Rows())); // Count the number of zero rows u32 zeroRows = A.Columns() - i; if (zeroRows == 0) // Unique solution: find it { // Back-substitute T sum; for (i--; i >= 0; i--) { sum = 0; for (j = i+1; j < A.Columns(); j++) sum += A(i,j)*b(j,0); b(i,0) = (b(i,0) - sum)/A(i,i); } // Fix the size of b, if necessary if (b.Rows() > A.Columns()) { Matrix<T> c(A.Columns(), 1); for (k = 0; k < A.Columns(); k++) c(k,0) = b(k,0); b = c; } return 1; } else // No solution or infinite solutions { // Check if there is no solutions for (k = i; k < A.Rows(); k++) if (b(k,0) != 0) return 0; // There is infinite solutions: find them if (basis == 0) return zeroRows + 1; else { // Get into reduced row echelon form bool Found; for (i--; i >= 0; i--) { Found = false; for (j = i; (j < A.Columns()) && (!Found); j++) if (A(i,j) != 0) { for (k = i-1; k >= 0; k--) A.SubRows(k, A(k,j)/A(i,j), i, j, A.Columns()); b(i,0) /= A(i,j); A.DivideRow(i, A(i,j)); Found = true; } } // Form the basis of solutions Matrix<T> solutions(A.Columns(), zeroRows + 1); solutions.LoadZeroes(); s32 lastColumn = 1; i = 0; j = 0; Found = false; do { if (Found || (A(i,j) == 0)) { for (k = 0; (k < A.Rows()) && (k < solutions.Rows()); k++) solutions(k,lastColumn) = -A(k,j); solutions(j,lastColumn) = 1; lastColumn++; j++; } else { solutions(j,0) = b(i,0); if (i < A.Rows()) i++; else Found = true; j++; } } while (j < A.Columns()); *basis = solutions; return zeroRows + 1; } } } template<class T> Matrix<T> Matrix<T>::SubmatrixWithout(u32 row, u32 col) const { Matrix<T> subMatrix(rows-1, columns-1); u32 x, y = 0; for (u32 i = 0; i < rows; i++) { x = 0; if (i != row) { for (u32 j = 0; j < columns; j++) { if (j != col) { subMatrix(x,y) = entry[i][j]; x++; } } y++; } } return subMatrix; } template <class T> void Matrix<T>::SubRows(u32 i, T pivot, u32 j, u32 start, u32 end) { for (u32 k = start; k < end; k++) entry[i][k] -= (pivot*entry[j][k]); } template <class T> void Matrix<T>::SwapRows(u32 i, u32 j) { Swap(entry[i], entry[j]); } template <class T> void Matrix<T>::Transpose() { T **newEntry = new T*[columns]; for (u32 i = 0; i < columns; i++) newEntry[i] = new T[rows]; for (u32 i = 0; i < columns; i++) for (u32 j = 0; j < rows; j++) newEntry[i][j] = entry[j][i]; for (u32 i = 0; i < rows; i++) delete[] entry[i]; delete[] entry; entry = newEntry; Swap(rows, columns); } template<class T> Matrix<T> Matrix<T>::Transposed() const { Matrix<T> result(columns, rows); for (u32 i = 0; i < rows; i++) for (u32 j = 0; j < columns; j++) result.entry[j][i] = entry[i][j]; return result; } template <class T> Matrix<T>::~Matrix() { for (unsigned long i = 0; i < rows; i++) delete[] entry[i]; delete[] entry; } int main() { using namespace std; Matrix<float> x(3,1); x(0,0) = 2; x(1,0) = 1; x(2,0) = 2; Matrix<float>::Householder(x).Display(); return 0; } (end of Matrix.cpp)-------------------------------------------
Same thing still with 3.3.3. Can you please attach the source + header and preprocessed source here, as it is easier to work with?
Created attachment 26453 [details] The header file used by Matrix.cpp
Created attachment 26454 [details] The file which causes the internal compiler error when compiled Try compiling with: gcc Matrix.cpp -o Matrix
This still fails for me with 3.3.5 and doesn't compile at all with 4.0.1. I think you should report it to gcc.gnu.org
