C++: 矩阵操作

在原作的基础上添加对矩阵扩充的支持：

void insert(unsigned row, unsigned col, TYPE const& value)；

matrix.h   ：

/*****************************************************************************/
/* Name: matrix.h                                                            */
/* Uses: Class for matrix math functions.                                    */
/* Date: 4/19/2011                                                           */
/* Author: Andrew Que <http://www.DrQue.net/>                                */
/* Revisions:                                                                */
/*   0.1 - 2011/04/19 - QUE - Creation.                                      */
/*   0.5 - 2011/04/24 - QUE - Most functions are complete.                   */
/*   0.8 - 2011/05/01 - QUE -                                                */
/*     = Bug fixes.                                                          */
/*     + Dot product.                                                        */
/*   1.0 - 2011/11/26 - QUE - Release.                                       */
/*                                                                           */
/* Notes:                                                                    */
/*   This unit implements some very basic matrix functions, which include:   */
/*    + Addition/subtraction                                                 */
/*    + Transpose                                                            */
/*    + Row echelon reduction                                                */
/*    + Determinant                                                          */
/*    + Dot product                                                          */
/*    + Matrix product                                                       */
/*    + Scalar product                                                       */
/*    + Inversion                                                            */
/*    + LU factorization/decomposition                                       */
/*     There isn't much for optimization in this unit as it was designed as  */
/*   more of a learning experience.                                          */
/*                                                                           */
/* License:                                                                  */
/*   This program is free software: you can redistribute it and/or modify    */
/*   it under the terms of the GNU General Public License as published by    */
/*   the Free Software Foundation, either version 3 of the License, or       */
/*   (at your option) any later version.                                     */
/*                                                                           */
/*   This program is distributed in the hope that it will be useful,         */
/*   but WITHOUT ANY WARRANTY; without even the implied warranty of          */
/*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           */
/*   GNU General Public License for more details.                            */
/*                                                                           */
/*   You should have received a copy of the GNU General Public License       */
/*   along with this program.  If not, see <http://www.gnu.org/licenses/>.   */
/*                                                                           */
/*                     (C) Copyright 2011 by Andrew Que                      */
/*                           http://www.DrQue.net/                           */
/*****************************************************************************/
#ifndef _MATRIX_H_
#define _MATRIX_H_

#include <iostream>
#include <cassert>
#include <climits>
#include <vector>

// Class forward for identity matrix.
template< class TYPE > class IdentityMatrix;

//=============================================================================
// Matrix template class
//   Contains a set of matrix manipulation functions.  The template is designed
// so that the values of the matrix can be of any type that allows basic
// arithmetic.
//=============================================================================
template< class TYPE = int >
  class Matrix
  {
    protected:
      // Matrix data.
      unsigned rows;
      unsigned columns;

      // Storage for matrix data.
      std::vector< std::vector< TYPE > > matrix;

      // Order sub-index for rows.
      //   Use: matrix[ order[ row ] ][ column ].
      unsigned * order;

      //-------------------------------------------------------------
      // Return the number of leading zeros in the given row.
      //-------------------------------------------------------------
      unsigned getLeadingZeros
      (
        // Row to count
        unsigned row
      ) const
      {
        TYPE const ZERO = static_cast< TYPE >( 0 );
        unsigned column = 0;
        while ( ZERO == matrix[ row ][ column ] )
          ++column;

        return column;
      }

      //-------------------------------------------------------------
      // Reorder the matrix so the rows with the most zeros are at
      // the end, and those with the least at the beginning.
      //
      // NOTE: The matrix data itself is not manipulated, just the
      // 'order' sub-indexes.
      //-------------------------------------------------------------
      void reorder()
      {
        unsigned * zeros = new unsigned[ rows ];

        for ( unsigned row = 0; row < rows; ++row )
        {
          order[ row ] = row;
          zeros[ row ] = getLeadingZeros( row );
        }

        for ( unsigned row = 0; row < (rows-1); ++row )
        {
          unsigned swapRow = row;
          for ( unsigned subRow = row + 1; subRow < rows; ++subRow )
          {
            if ( zeros[ order[ subRow ] ] < zeros[ order[ swapRow ] ] )
              swapRow = subRow;
          }

          unsigned hold    = order[ row ];
          order[ row ]     = order[ swapRow ];
          order[ swapRow ] = hold;
        }

        delete zeros;
      }

      //-------------------------------------------------------------
      // Divide a row by given value.  An elementary row operation.
      //-------------------------------------------------------------
      void divideRow
      (
        // Row to divide.
        unsigned row,

        // Divisor.
        TYPE const & divisor
      )
      {
        for ( unsigned column = 0; column < columns; ++column )
          matrix[ row ][ column ] /= divisor;
      }

      //-------------------------------------------------------------
      // Modify a row by adding