matrix.h

复数矩阵的运算

////////////////////////////////
// Matrix TCL Lite v1.13
// Copyright (c) 1997-2002 Techsoft Pvt. Ltd. (See License.Txt file.)
//
// Matrix.h: Matrix C++ template class include file 
// Web: http://www.techsoftpl.com/matrix/
// Email: [e-mail]
//

//////////////////////////////
// Installation:
//
// Copy this "matrix.h" file into include directory of your compiler.
//

//////////////////////////////
// Note: This matrix template class defines majority of the matrix
// operations as overloaded operators or methods. It is assumed that
// users of this class is familiar with matrix algebra. We have not
// defined any specialization of this template here, so all the instances
// of matrix will be created implicitly by the compiler. The data types
// tested with this class are float, double, long double, complex<float>,
// complex<double> and complex<long double>. Note that this class is not 
// optimized for performance.
//
// Since implementation of exception, namespace and template are still
// not standardized among the various (mainly old) compilers, you may 
// encounter compilation error with some compilers. In that case remove 
// any of the above three features by defining the following macros:
//
//  _NO_NAMESPACE:  Define this macro to remove namespace support.
//
//  _NO_EXCEPTION:  Define this macro to remove exception handling
//                  and use old style of error handling using function.
//
//  _NO_TEMPLATE:   If this macro is defined matrix class of double
//                  type will be generated by default. You can also
//                  generate a different type of matrix like float.
//
//  _SGI_BROKEN_STL: For SGI C++ v.7.2.1 compiler.
//
//  Since all the definitions are also included in this header file as
//  inline function, some compiler may give warning "inline function
//  can't be expanded". You may ignore/disable this warning using compiler
//  switches. All the operators/methods defined in this class have their
//  natural meaning except the followings:
//
//  Operator/Method                          Description
//  ---------------                          -----------
//   operator ()   :   This function operator can be used as a
//                     two-dimensional subscript operator to get/set
//                     individual matrix elements.
//
//   operator !    :   This operator has been used to calculate inversion
//                     of matrix.
//
//   operator ~    :   This operator has been used to return transpose of
//                     a matrix.
//
//   operator ^    :   It is used calculate power (by a scalar) of a matrix.
//                     When using this operator in a matrix equation, care
//                     must be taken by parenthesizing it because it has
//                     lower precedence than addition, subtraction,
//                     multiplication and division operators.
//
//   operator >>   :   It is used to read matrix from input stream as per
//                     standard C++ stream operators.
//
//   operator <<   :   It is used to write matrix to output stream as per
//                     standard C++ stream operators.
//
// Note that professional version of this package, Matrix TCL Pro 2.11
// is optimized for performance and supports many more matrix operations.
// It is available from our web site at <http://www.techsoftpl.com/matrix/>.
//

#ifndef __cplusplus
#error Must use C++ for the type matrix.
#endif

#if !defined(__STD_MATRIX_H)
#define __STD_MATRIX_H

//////////////////////////////
// First deal with various shortcomings and incompatibilities of
// various (mainly old) versions of popular compilers available.
//

#if defined(__BORLANDC__)
#pragma option -w-inl -w-pch
#endif

#if ( defined(__BORLANDC__) || _MSC_VER <= 1000 ) && !defined( __GNUG__ )
#  include <stdio.h>
#  include <stdlib.h>
#  include <math.h>
#  include <iostream.h>
#  include <string.h>
#else
#  include <cmath>
#  include <cstdio>
#  include <cstdlib>
#  include <string>
#  include <iostream>
#endif

#if defined(_MSC_VER) && _MSC_VER <= 1000
#  define _NO_EXCEPTION        // stdexception is not fully supported in MSVC++ 4.0
typedef int bool;
#  if !defined(false)
#    define false  0
#  endif
#  if !defined(true)
#    define true   1
#  endif
#endif

#if defined(__BORLANDC__) && !defined(__WIN32__)
#  define _NO_EXCEPTION        // std exception and namespace are not fully
#  define _NO_NAMESPACE        // supported in 16-bit compiler
#endif

#if defined(_MSC_VER) && !defined(_WIN32)
#  define _NO_EXCEPTION
#endif

#if defined(_NO_EXCEPTION)
#  define _NO_THROW
#  define _THROW_MATRIX_ERROR
#else
#  if defined(_MSC_VER)
#    if _MSC_VER >= 1020
#      include <stdexcept>
#    else
#      include <stdexcpt.h>
#    endif
#  elif defined(__MWERKS__)
#      include <stdexcept>
#  elif (__GNUC__ >= 2 || (__GNUC__ == 2 && __GNUC_MINOR__ >= 8))
#     include <stdexcept>
#  else
#     include <stdexcep>
#  endif
#  define _NO_THROW               throw ()
#  define _THROW_MATRIX_ERROR     throw (matrix_error)
#endif

#ifndef __MINMAX_DEFINED
#  define max(a,b)    (((a) > (b)) ? (a) : (b))
#  define min(a,b)    (((a) < (b)) ? (a) : (b))
#endif

#if defined(_MSC_VER)
#undef _MSC_EXTENSIONS      // To include overloaded abs function definitions!
#endif

//#if ( defined(__BORLANDC__) || _MSC_VER ) && !defined( __GNUG__ ) 
//inline float abs (float v) { return (float)fabs( v); } 
//inline double abs (double v) { return fabs( v); } 
//inline long double abs (long double v) { return fabsl( v); }
//#endif

#if defined(__GNUG__) || defined(__MWERKS__) || (defined(__BORLANDC__) && (__BORLANDC__ >= 0x540))
#define FRIEND_FUN_TEMPLATE <>
#else
#define FRIEND_FUN_TEMPLATE
#endif

#if defined(_MSC_VER) && _MSC_VER <= 1020   // MSVC++ 4.0/4.2 does not
#  define _NO_NAMESPACE                     // support "std" namespace
#endif

#if !defined(_NO_NAMESPACE)
#if defined( _SGI_BROKEN_STL )              // For SGI C++ v.7.2.1 compiler
namespace std { }
#endif
using namespace std;
#endif

#ifndef _NO_NAMESPACE
namespace math {
#endif

#if !defined(_NO_EXCEPTION)
    class matrix_error : public logic_error
    {
    public:
        matrix_error (const string& what_arg) : logic_error( what_arg) {}
    };
#define REPORT_ERROR(ErrormMsg)  throw matrix_error( ErrormMsg);
#else
    inline void _matrix_error (const char* pErrMsg)
    {
        cout << pErrMsg << endl;
        exit(1);
    }
#define REPORT_ERROR(ErrormMsg)  _matrix_error( ErrormMsg);
#endif

#if !defined(_NO_TEMPLATE)
#  define MAT_TEMPLATE  template <class T>
#  define matrixT  matrix<T>
#else
#  define MAT_TEMPLATE
#  define matrixT  matrix
#  ifdef MATRIX_TYPE
    typedef MATRIX_TYPE T;
#  else
    typedef double T;
#  endif
#endif


    MAT_TEMPLATE
    class matrix
    {
    public:
        // Constructors
        // 构造函数
        matrix (const matrixT& m);
        matrix (size_t row = 0, size_t col = 0);

        // Destructor
        ~matrix ();

        // Assignment operators
        matrixT& operator = (const matrixT& m) _NO_THROW;

        // Value extraction method
        size_t get_rowsize () const { return _m->Row; }
        size_t get_colsize () const { return _m->Col; }

        // Subscript operator
        // 下标操作符
        T& operator [][] (size_t row, size_t col) _THROW_MATRIX_ERROR;
        T  operator [][] (size_t row, size_t col) const _THROW_MATRIX_ERROR;

        // Unary operators
        matrixT operator + () _NO_THROW { return *this; }
        matrixT operator - () _NO_THROW;

        // Combined assignment - calculation operators
        matrixT& operator += (const matrixT& m) _THROW_MATRIX_ERROR;
        matrixT& operator -= (const matrixT& m) _THROW_MATRIX_ERROR;
        matrixT& operator *= (const matrixT& m) _THROW_MATRIX_ERROR;
        matrixT& operator *= (const T& c) _NO_THROW;
        matrixT& operator /= (const T& c) _NO_THROW;
        matrixT& operator ^= (const size_t& pow) _THROW_MATRIX_ERROR;
        
        // Miscellaneous -methods
        // 清零
        void Null (const size_t& row, const size_t& col) _NO_THROW;
        void Null () _NO_THROW;
        // 置成单位矩阵
        void Unit (const size_t& row) _NO_THROW;
        void Unit () _NO_THROW;
        void SetSize (size_t row, size_t col) _NO_THROW;

        // Utility methods
        matrixT Solve (const matrixT& v) const _THROW_MATRIX_ERROR;
        matrixT Adj () _THROW_MATRIX_ERROR;
        matrixT Inv () _THROW_MATRIX_ERROR;
        T Det () const _THROW_MATRIX_ERROR;
        T Norm () _NO_THROW;
        T Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR;
        T Cond () _NO_THROW;

        // Type of matrices
        bool IsSquare () _NO_THROW { return (_m->Row == _m->Col); } 
        bool IsSingular () _NO_THROW;
        bool IsDiagonal () _NO_THROW;
        bool IsScalar () _NO_THROW;
        bool IsUnit () _NO_THROW;
        bool IsNull () _NO_THROW;
        bool IsSymmetric () _NO_THROW;
        bool IsSkewSymmetric () _NO_THROW;
        bool IsUpperTriangular () _NO_THROW;
        bool IsLowerTriangular () _NO_THROW;

    private:
        struct base_mat
        {
            T **Val;
            size_t Row, Col, RowSiz, ColSiz;
            int Refcnt;

            base_mat (size_t row, size_t col, T** v)
            {
                Row = row; RowSiz = row;
                Col = col; ColSiz = col;
                Refcnt = 1;

                Val = new T* [row];
                size_t rowlen = col * sizeof(T);

                for (size_t i=0; i < row; i++)
                {
                    Val[i] = new T [col];
                    if (v) memcpy( Val[i], v[i], rowlen);
                }
            }
            ~base_mat ()
            {
                for (size_t i=0; i < RowSiz; i++)
                    delete [] Val[i];
                delete [] Val;
            }
        };
        base_mat *_m;

        void clone ();
        void realloc (size_t row, size_t col);
        int pivot (size_t row);
    };

#if defined(_MSC_VER) && _MSC_VER <= 1020
#  undef  _NO_THROW               // MSVC++ 4.0/4.2 does not support 
#  undef  _THROW_MATRIX_ERROR     // exception specification in definition
#  define _NO_THROW
#  define _THROW_MATRIX_ERROR
#endif

    // constructor
    MAT_TEMPLATE inline
        matrixT::matrix (size_t row, size_t col)
    {
        _m = new base_mat( row, col, 0);
    }

    // copy constructor
    MAT_TEMPLATE inline
        matrixT::matrix (const matrixT& m)
    {
        _m = m._m;
        _m->Refcnt++;
    }

    // Internal copy constructor
    MAT_TEMPLATE inline void
        matrixT::clone ()
    {
        _m->Refcnt--;
        _m = new base_mat( _m->Row, _m->Col, _m->Val);
    }

    // destructor
    MAT_TEMPLATE inline
        matrixT::~matrix ()
    {
        if (--_m->Refcnt == 0) delete _m;
    }

    // assignment operator
    MAT_TEMPLATE inline matrixT&
        matrixT::operator = (const matrixT& m) _NO_THROW
    {
        m._m->Refcnt++;
        if (--_m->Refcnt == 0) delete _m;