转载矩阵运算代码

wxhwhmanshan

文章标题：矩阵运算代码
原 作 者：黄运新
原 出 处：不详
发 布 者：wxh
发布类型：转载

功能分析：
本程序能完成矩阵的输入、输出。具有相同行数和列数的矩阵间的加法、减法。符合矩阵乘法规则要求的矩阵间的乘法。方阵间的除法，方阵的求逆。矩阵的求转置矩阵等功能。
设计思想：
本程序要完成矩阵的一般运算。同时本程序的截面要友好，操作要简单，并具有较强的自我适应能力，即根据每人输入矩阵的方法不同，在生成矩阵时，程序能适当根据输入的各种格式来达到输入矩阵的功能；
本程序为每个矩阵变量都分配一个矩阵变量名，用户可以通过操作矩阵变量来达到运算目的。如果用户对同一个矩阵变量定义了多次，则程序能根据系统内是否存在这个矩阵变量名，来给同一个矩阵变量进行多次赋值。
本程序能根据矩阵表达式来进行表达式求解。也就是说程序能判定各个运算符的优先权，然后进行分析计算。
结构分析：
class CArrayMatrix : public CObject
该类是Cobject类的派生类。是用来存储数据和进行一般矩阵运算的类。在vc用户手册中强烈建议如果程序自己定义的类，最要是从Cobject派生出来。详细原因请参考MSDN。
数据成员：
CArray<float,float &> m_Array;
CString m_Name;
UINT m_wRow;
UINT m_wCon;
m_Array是CArray类的一个对象，用来保存矩阵中每个元素的值。
m_Name是一个字符串类的对象，用来保存矩阵对象的名称。
m_wRow是用来保存矩阵的行数。
m_wCon是用来保存矩阵的列数。
成员函数：
BOOL CArrayMatrix::CanAddSub(const CArrayMatrix &m)const
这个函数是用来判断矩阵是否可以相加。如果是返回TRUE；反之返回FALSE。
BOOL CArrayMatrix::CanContrary()const
这个函数是用来判断矩阵是否可逆。如果是返回TRUE；反之返回FALSE。
BOOL CArrayMatrix::CanMutiply(const CArrayMatrix &m)const
这个函数是用来判断矩阵是否可相乘。如果是返回TRUE；反之返回FALSE。
CArrayMatrix::CArrayMatrix()
无参构造函数。
CArrayMatrix::CArrayMatrix(const CArrayMatrix & m)
拷贝构造函数。进行深拷贝。
CArrayMatrix::~CArrayMatrix()
析构函数。
BOOL CArrayMatrix::ConvertToNum
(const CString & sRowString, const int & nAtCon, float & fResult)//nAtCon从0---m_wCon-1
此函数是把字符串中第nAtCon个子字符串转化成浮点型。
void CArrayMatrix:eleteHeadEnter(CString & RightString)
void CArrayMatrix:rocString(CString & RightString)
void CArrayMatrix::ProcStringBlank(CString & sHead)
void CArrayMatrix::ProcStringComma(CString & sHead)
void CArrayMatrix::ProcStringEnter(CString & sHead)
这些函数是用来过滤输入字符串的。
BOOL CArrayMatrix::DisPlay(CString & out)const
显示矩阵函数，把结果存到out里。
CArrayMatrix CArrayMatrix::GetAccompany()const
得到矩阵的伴随矩阵，并返回
CString CArrayMatrix::GetName()
得到矩阵的名称
float CArrayMatrix::GetRange() const
得到矩阵的行列式的值。
BOOL CArrayMatrix::InputCon(const CString & RS,CString * &amptRSBuffer)
根据过滤函数得出的字符串来给矩阵的类数赋值，并返回一个指向子字符串项的指针。
BOOL CArrayMatrix::InputRow(const CString & RS)
得到矩阵的行数。
CArrayMatrix CArrayMatrix:perator *(float m)
CArrayMatrix CArrayMatrix::operator *(const CArrayMatrix &m)
如果乘法运算符重载。
CArrayMatrix CArrayMatrix::operator +(const CArrayMatrix & m)
矩阵加法运算符重载。
CArrayMatrix CArrayMatrix::operator -(const CArrayMatrix & m)
矩阵减法运算符重载。
CArrayMatrix CArrayMatrix::operator /(float m)
CArrayMatrix CArrayMatrix::operator / (CArrayMatrix &m)
矩阵除法运算符重载。
CArrayMatrix & CArrayMatrix::operator =(float m)
CArrayMatrix & CArrayMatrix::operator =(const CArrayMatrix & m)
矩阵赋值运算符重载。
CArrayMatrix CArrayMatrix::operator ~()
求逆运算符重载。
BOOL CArrayMatrix::SetName(CString m)
给矩阵对象名称赋值。
BOOL CArrayMatrix::SetStringName(CString &DataString,CString & sName)
自动截断输入字符串等号，右边是矩阵对象的数据部分，左边是矩阵对象的名称。
CArrayMatrix CArrayMatrix::T()const
求转置矩阵的函数。
template<class T>
class CStack:public CList<T,T &>
标准栈类，是从链表类派生出来的。用来进行表达式求解时用到。
程序测试：
1、测试赋值和加法。
输入数据：
a=2 4 6 7
  6 9 3 4
  5 7 8 3
  3 2 5 7;
b=1 0 -1 2
-1 1 3 0
0 5 7 -6
6 8 0 4;
c=a+b;
c=;
输出结果：
c=
3   4   5   9   
5   10   6   4   
5   12   15   -3   
9   10   5   11   ;
2、减法测试
输入数据：
d=c-b;
d=;
输出结果：
d=
2   4   6   7   
6   9   3   4   
5   7   8   3   
3   2   5   7   ;
3、乘法测试
输入数据：
a=1 0 -1 2
-1 1 3 0
0 5 7 -6;
b=0 -4
1 2
-3 -2
-1 1;
c=a*b;
c=;
输入结果：
c=
1    0   
-8   0   
-10   -10   ;
4、求转置测试：
输入：
d=c';
d=;
输出：
d=
1    -8   -10   
0    0    -10   ;
5、求逆
输入：
a=3 2 1
  1 1 1
  1 0 1;
b=~a;
b=;
输出：
b=
0.5   -1   0.5   
0    1    -1   
-0.5  1    0.5   ;
6、除法
输入：
a=1 2 -3
  3 2 -4
  2 -1 0;
b=1 -3 0
  10 2 7
  10 7 8;
c=b/a;
c=;
输出：
c=
20    -15   13   
-105   77    -58   
-152   112   -87   ;

Pointer3D

代码哩?

SevenCat

// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2004.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

#ifndef WMLMATRIX_H
#define WMLMATRIX_H

// Matrix operations are applied on the left.  For example, given a matrix M
// and a vector V, matrix-times-vector is M*V.  That is, V is treated as a
// column vector.  Some graphics APIs use V*M where V is treated as a row
// vector.  In this context the "M" matrix is really a transpose of the M as
// represented in Wild Magic.  Similarly, to apply two matrix operations M0
// and M1, in that order, you compute M1*M0 so that the transform of a vector
// is (M1*M0)*V = M1*(M0*V).  Some graphics APIs use M0*M1, but again these
// matrices are the transpose of those as represented in Wild Magic.  You
// must therefore be careful about how you interface the transformation code
// with graphics APIS.
//
// For memory organization it might seem natural to chose Real[N][N] for the
// matrix storage, but this can be a problem on a platform/console that
// chooses to store the data in column-major rather than row-major format.
// To avoid potential portability problems, the matrix is stored as Real[N*N]
// and organized in row-major order.  That is, the entry of the matrix in row
// r (0 <= r < N) and column c (0 <= c < N) is stored at index i = c+N*r
// (0 <= i < N*N).

#include "WmlVector.h"

namespace Wml
{

template <int N, class Real>
class WML_ITEM Matrix
{
public:
    // Construction.  In the second constructor, if bZero is 'true', the
    // matrix is set to the zero matrix.  If bZero is 'false', the matrix
    // is set to the identity matrix.
    Matrix ();
    Matrix (bool bZero);
    Matrix (const Matrix& rkM);

    void MakeZero ();
    void MakeIdentity ();
    void MakeDiagonal (const Real* afDiag);

    // member access
    operator const Real* () const;
    operator Real* ();
    const Real* operator[] (int iRow) const;
    Real* operator[] (int iRow);
    Real operator() (int iRow, int iCol) const;
    Real& operator() (int iRow, int iCol);
    void SetRow (int iRow, const Vector<N,Real>& rkV);
    Vector<N,Real> GetRow (int iRow) const;
    void SetColumn (int iCol, const Vector<N,Real>& rkV);
    Vector<N,Real> GetColumn (int iCol) const;
    void GetColumnMajor (Real* afCMajor) const;

    // assignment
    Matrix& operator= (const Matrix& rkM);

    // comparison
    bool operator== (const Matrix& rkM) const;
    bool operator!= (const Matrix& rkM) const;
    bool operator<  (const Matrix& rkM) const;
    bool operator<= (const Matrix& rkM) const;
    bool operator>  (const Matrix& rkM) const;
    bool operator>= (const Matrix& rkM) const;

    // arithmetic operations
    Matrix operator+ (const Matrix& rkM) const;
    Matrix operator- (const Matrix& rkM) const;
    Matrix operator* (const Matrix& rkM) const;
    Matrix operator* (Real fScalar) const;
    Matrix operator/ (Real fScalar) const;
    Matrix operator- () const;

    // arithmetic updates
    Matrix& operator+= (const Matrix& rkM);
    Matrix& operator-= (const Matrix& rkM);
    Matrix& operator*= (Real fScalar);
    Matrix& operator/= (Real fScalar);

    // matrix products
    Matrix Transpose () const;  // M^T
    Matrix TransposeTimes (const Matrix& rkM) const;  // this^T * M
    Matrix TimesTranspose (const Matrix& rkM) const;  // this * M^T

    // matrix-vector operations
    Vector<N,Real> operator* (const Vector<N,Real>& rkV) const;  // M * v
    Real QForm (const Vector<N,Real>& rkU, const Vector<N,Real>& rkV)
        const;  // u^T*M*v

protected:
    // for indexing into the 1D array of the matrix, iCol+N*iRow
    static int I (int iRow, int iCol);

    // support for comparisons
    int CompareArrays (const Matrix& rkM) const;

    Real m_afEntry[N*N];
};

// c * M
template <int N, class Real>
WML_ITEM Matrix<N,Real> operator* (Real fScalar, const Matrix<N,Real>& rkM);

// v^T * M
template <int N, class Real>
WML_ITEM Vector<N,Real> operator* (const Vector<N,Real>& rkV,
    const Matrix<N,Real>& rkM);

}

#endif
// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2004.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>::Matrix ()
{
    // the matrix is uninitialized
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>::Matrix (bool bZero)
{
    memset(m_afEntry,0,N*N*sizeof(Real));

    if ( !bZero )
    {
        for (int i = 0; i < N; i++)
            m_afEntry[I(i,i)] = (Real)1.0;
    }
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>::Matrix (const Matrix& rkM)
{
    memcpy(m_afEntry,rkM.m_afEntry,N*N*sizeof(Real));
}
//----------------------------------------------------------------------------
template <int N, class Real>
void Matrix<N,Real>::MakeZero ()
{
    memset(m_afEntry,0,N*N*sizeof(Real));
}
//----------------------------------------------------------------------------
template <int N, class Real>
void Matrix<N,Real>::MakeIdentity ()
{
    memset(m_afEntry,0,N*N*sizeof(Real));
    for (int i = 0; i < N; i++)
        m_afEntry[I(i,i)] = (Real)1.0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
void Matrix<N,Real>::MakeDiagonal (const Real* afDiag)
{
    memset(m_afEntry,0,N*N*sizeof(Real));
    for (int i = 0; i < N; i++)
        m_afEntry[I(i,i)] = afDiag;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>:perator const Real* () const
{
    return m_afEntry;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>::operator Real* ()
{
    return m_afEntry;
}
//----------------------------------------------------------------------------
template <int N, class Real>
const Real* Matrix<N,Real>::operator[] (int iRow) const
{
    assert( 0 <= iRow && iRow < N );
    return &m_afEntry[N*iRow];
}
//----------------------------------------------------------------------------
template <int N, class Real>
Real* Matrix<N,Real>::operator[] (int iRow)
{
    assert( 0 <= iRow && iRow < N );
    return &m_afEntry[N*iRow];
}
//----------------------------------------------------------------------------
template <int N, class Real>
Real Matrix<N,Real>::operator() (int iRow, int iCol) const
{
    return m_afEntry[I(iRow,iCol)];
}
//----------------------------------------------------------------------------
template <int N, class Real>
Real& Matrix<N,Real>::operator() (int iRow, int iCol)
{
    return m_afEntry[I(iRow,iCol)];
}
//----------------------------------------------------------------------------
template <int N, class Real>
void Matrix<N,Real>::SetRow (int iRow, const Vector<N,Real>& rkV)
{
    assert( 0 <= iRow && iRow < N );
    for (int iCol = 0, i = N*iRow; iCol < N; iCol++, i++)
        m_afEntry = rkV[iCol];
}
//----------------------------------------------------------------------------
template <int N, class Real>
Vector<N,Real> Matrix<N,Real>::GetRow (int iRow) const
{
    assert( 0 <= iRow && iRow < N );
    Vector<N,Real> kV;
    for (int iCol = 0, i = N*iRow; iCol < N; iCol++, i++)
        kV[iCol] = m_afEntry;
    return kV;
}
//----------------------------------------------------------------------------
template <int N, class Real>
void Matrix<N,Real>::SetColumn (int iCol, const Vector<N,Real>& rkV)
{
    assert( 0 <= iCol && iCol < N );
    for (int iRow = 0, i = iCol; iRow < N; iRow++, i += N)
        m_afEntry = rkV[iRow];
}
//----------------------------------------------------------------------------
template <int N, class Real>
Vector<N,Real> Matrix<N,Real>::GetColumn (int iCol) const
{
    assert( 0 <= iCol && iCol < N );
    Vector<N,Real> kV;
    for (int iRow = 0, i = iCol; iRow < N; iRow++, i += N)
        kV[iRow] = m_afEntry;
    return kV;
}
//----------------------------------------------------------------------------
template <int N, class Real>
void Matrix<N,Real>::GetColumnMajor (Real* afCMajor) const
{
    for (int iRow = 0, i = 0; iRow < N; iRow++)
    {
        for (int iCol = 0; iCol < N; iCol++)
            afCMajor[i++] = m_afEntry[I(iCol,iRow)];
    }
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>& Matrix<N,Real>::operator= (const Matrix& rkM)
{
    memcpy(m_afEntry,rkM.m_afEntry,N*N*sizeof(Real));
    return *this;
}
//----------------------------------------------------------------------------
template <int N, class Real>
bool Matrix<N,Real>::operator== (const Matrix& rkM) const
{
    return memcmp(m_afEntry,rkM.m_afEntry,N*N*sizeof(Real)) == 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
bool Matrix<N,Real>::operator!= (const Matrix& rkM) const
{
    return memcmp(m_afEntry,rkM.m_afEntry,N*N*sizeof(Real)) != 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
int Matrix<N,Real>::CompareArrays (const Matrix& rkM) const
{
    for (int i = 0; i < N*N; i++)
    {
        unsigned int uiTest0 = *(unsigned int*)&m_afEntry;
        unsigned int uiTest1 = *(unsigned int*)&rkM.m_afEntry;
        if ( uiTest0 < uiTest1 )
            return -1;
        if ( uiTest0 > uiTest1 )
            return +1;
    }
    return 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
bool Matrix<N,Real>::operator<  (const Matrix& rkM) const
{
    return CompareArrays(rkM) < 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
bool Matrix<N,Real>::operator<= (const Matrix& rkM) const
{
    return CompareArrays(rkM) <= 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
bool Matrix<N,Real>::operator>  (const Matrix& rkM) const
{
    return CompareArrays(rkM) > 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
bool Matrix<N,Real>::operator>= (const Matrix& rkM) const
{
    return CompareArrays(rkM) >= 0;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::operator+ (const Matrix& rkM) const
{
    Matrix<N,Real> kSum;
    for (int i = 0; i < N*N; i++)
        kSum.m_afEntry = m_afEntry + rkM.m_afEntry;
    return kSum;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::operator- (const Matrix& rkM) const
{
    Matrix<N,Real> kDiff;
    for (int i = 0; i < N*N; i++)
        kDiff.m_afEntry = m_afEntry - rkM.m_afEntry;
    return kDiff;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::operator* (const Matrix& rkM) const
{
    Matrix<N,Real> kProd;
    for (int iRow = 0; iRow < N; iRow++)
    {
        for (int iCol = 0; iCol < N; iCol++)
        {
            int i = I(iRow,iCol);
            kProd.m_afEntry = (Real)0.0f;
            for (int iMid = 0; iMid < N; iMid++)
            {
                kProd.m_afEntry +=
                    m_afEntry[I(iRow,iMid)]*rkM.m_afEntry[I(iMid,iCol)];
            }
        }
    }
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::operator* (Real fScalar) const
{
    Matrix<N,Real> kProd;
    for (int i = 0; i < N*N; i++)
        kProd.m_afEntry = fScalar*m_afEntry;
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::operator/ (Real fScalar) const
{
    Matrix<N,Real> kQuot;
    int i;

    if ( fScalar != (Real)0.0 )
    {
        Real fInvScalar = ((Real)1.0)/fScalar;
        for (i = 0; i < N*N; i++)
            kQuot.m_afEntry = fInvScalar*m_afEntry;
    }
    else
    {
        for (i = 0; i < N*N; i++)
            kQuot.m_afEntry = Math<Real>::MAX_REAL;
    }

    return kQuot;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::operator- () const
{
    Matrix<N,Real> kNeg;
    for (int i = 0; i < N*N; i++)
        kNeg.m_afEntry = -m_afEntry;
    return kNeg;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> operator* (Real fScalar, const Matrix<N,Real>& rkM)
{
    Matrix<N,Real> kProd;
    const Real* afMEntry = rkM;
    Real* afPEntry = kProd;
    for (int i = 0; i < N*N; i++)
        afPEntry = fScalar*afMEntry;
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>& Matrix<N,Real>::operator+= (const Matrix& rkM)
{
    for (int i = 0; i < N*N; i++)
        m_afEntry += rkM.m_afEntry;
    return *this;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>& Matrix<N,Real>::operator-= (const Matrix& rkM)
{
    for (int i = 0; i < N*N; i++)
        m_afEntry -= rkM.m_afEntry;
    return *this;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>& Matrix<N,Real>::operator*= (Real fScalar)
{
    for (int i = 0; i < N*N; i++)
        m_afEntry *= fScalar;
    return *this;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real>& Matrix<N,Real>::operator/= (Real fScalar)
{
    int i;

    if ( fScalar != (Real)0.0 )
    {
        Real fInvScalar = ((Real)1.0)/fScalar;
        for (i = 0; i < N*N; i++)
            m_afEntry *= fInvScalar;
    }
    else
    {
        for (i = 0; i < N*N; i++)
            m_afEntry = Math<Real>::MAX_REAL;
    }

    return *this;
}

SevenCat

//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::Transpose () const
{
    Matrix<N,Real> kTranspose;
    for (int iRow = 0; iRow < N; iRow++)
    {
        for (int iCol = 0; iCol < N; iCol++)
            kTranspose.m_afEntry[I(iRow,iCol)] = m_afEntry[I(iCol,iRow)];
    }
    return kTranspose;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::TransposeTimes (const Matrix& rkM) const
{
    // P = A^T*B, P[r][c] = sum_m A[m][r]*B[m][c]
    Matrix<N,Real> kProd;
    for (int iRow = 0; iRow < N; iRow++)
    {
        for (int iCol = 0; iCol < N; iCol++)
        {
            int i = I(iRow,iCol);
            kProd.m_afEntry = (Real)0.0;
            for (int iMid = 0; iMid < N; iMid++)
            {
                kProd.m_afEntry +=
                    m_afEntry[I(iMid,iRow)]*rkM.m_afEntry[I(iMid,iCol)];
            }
        }
    }
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Matrix<N,Real> Matrix<N,Real>::TimesTranspose (const Matrix& rkM) const
{
    // P = A*B^T, P[r][c] = sum_m A[r][m]*B[c][m]
    Matrix<N,Real> kProd;
    for (int iRow = 0; iRow < N; iRow++)
    {
        for (int iCol = 0; iCol < N; iCol++)
        {
            int i = I(iRow,iCol);
            kProd.m_afEntry = (Real)0.0;
            for (int iMid = 0; iMid < N; iMid++)
            {
                kProd.m_afEntry +=
                    m_afEntry[I(iRow,iMid)]*rkM.m_afEntry[I(iCol,iRow)];
            }
        }
    }
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Vector<N,Real> Matrix<N,Real>:perator* (const Vector<N,Real>& rkV) const
{
    Vector<N,Real> kProd;
    for (int iRow = 0; iRow < N; iRow++)
    {
        kProd[iRow] = (Real)0.0;
        for (int iCol = 0; iCol < N; iCol++)
            kProd[iRow] += m_afEntry[I(iRow,iCol)]*rkV[iCol];
            
    }
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Vector<N,Real> operator* (const Vector<N,Real>& rkV,
    const Matrix<N,Real>& rkM)
{
    Vector<N,Real> kProd;
    const Real* afMEntry = rkM;
    Real* afPEntry = kProd;
    for (int iCol = 0; iCol < N; iCol++)
    {
        afPEntry[iCol] = (Real)0.0;
        for (int iRow = 0; iRow < N; iRow++)
            afPEntry[iCol] += rkV[iRow]*afMEntry[iCol + N*iRow];
            
    }
    return kProd;
}
//----------------------------------------------------------------------------
template <int N, class Real>
Real Matrix<N,Real>:Form (const Vector<N,Real>& rkU,
    const Vector<N,Real>& rkV) const
{
    return rkU.Dot((*this)*rkV);
}
//----------------------------------------------------------------------------
template <int N, class Real>
int Matrix<N,Real>::I (int iRow, int iCol)
{
    assert( 0 <= iRow && iRow < N && 0 <= iCol && iCol < N );
    return iCol + N*iRow;
}
//----------------------------------------------------------------------------
// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2004.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

#ifndef WMLMATRIX4_H
#define WMLMATRIX4_H

#include "WmlMatrix.h"
#include "WmlVector4.h"

namespace Wml
{

template <class Real>
class WML_ITEM Matrix4 : public Matrix<4,Real>
{
public:
    // construction
    Matrix4 ();
    Matrix4 (const Matrix4& rkM);
    Matrix4 (const Matrix<4,Real>& rkM);

    // input Mrc is in row r, column c.
    Matrix4 (Real fM00, Real fM01, Real fM02, Real fM03,
             Real fM10, Real fM11, Real fM12, Real fM13,
             Real fM20, Real fM21, Real fM22, Real fM23,
             Real fM30, Real fM31, Real fM32, Real fM33);

    // Create a matrix from an array of numbers.  The input array is
    // interpreted based on the Boolean input as
    //   true:  entry[0..15]={m00,m01,m02,m03,m10,m11,m12,m13,m20,m21,m22,
    //                        m23,m30,m31,m32,m33} [row major]
    //   false: entry[0..15]={m00,m10,m20,m30,m01,m11,m21,m31,m02,m12,m22,
    //                        m32,m03,m13,m23,m33} [col major]
    Matrix4 (const Real afEntry[16], bool bRowMajor);

    // assignment
    Matrix4& operator= (const Matrix4& rkM);
    Matrix4& operator= (const Matrix<4,Real>& rkM);

    // matrix operations
    Matrix4 Inverse () const;
    Matrix4 Adjoint () const;
    Real Determinant () const;

    // special matrices
    static const Matrix4 ZERO;
    static const Matrix4 IDENTITY;
};

typedef Matrix4<float> Matrix4f;
typedef Matrix4<double> Matrix4d;

}

#endif
// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2004.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

#include "WmlMatrix4.h"
namespace Wml
{
#include "WmlMatrix.inl"
}
using namespace Wml;

//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4 ()
{
    // the matrix is uninitialized
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4 (const Matrix4& rkM)
{
    memcpy(m_afEntry,rkM.m_afEntry,16*sizeof(Real));
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4 (const Matrix<4,Real>& rkM)
{
    memcpy(m_afEntry,(const Real*)rkM,16*sizeof(Real));
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>& Matrix4<Real>::operator= (const Matrix4& rkM)
{
    memcpy(m_afEntry,rkM.m_afEntry,16*sizeof(Real));
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>& Matrix4<Real>::operator= (const Matrix<4,Real>& rkM)
{
    memcpy(m_afEntry,(const Real*)rkM,16*sizeof(Real));
    return *this;
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4 (Real fM00, Real fM01, Real fM02, Real fM03,
    Real fM10, Real fM11, Real fM12, Real fM13, Real fM20, Real fM21,
    Real fM22, Real fM23, Real fM30, Real fM31, Real fM32, Real fM33)
{
    m_afEntry[ 0] = fM00;
    m_afEntry[ 1] = fM01;
    m_afEntry[ 2] = fM02;
    m_afEntry[ 3] = fM03;
    m_afEntry[ 4] = fM10;
    m_afEntry[ 5] = fM11;
    m_afEntry[ 6] = fM12;
    m_afEntry[ 7] = fM13;
    m_afEntry[ 8] = fM20;
    m_afEntry[ 9] = fM21;
    m_afEntry[10] = fM22;
    m_afEntry[11] = fM23;
    m_afEntry[12] = fM30;
    m_afEntry[13] = fM31;
    m_afEntry[14] = fM32;
    m_afEntry[15] = fM33;
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4 (const Real afEntry[16], bool bRowMajor)
{
    if ( bRowMajor )
    {
        memcpy(m_afEntry,afEntry,16*sizeof(Real));
    }
    else
    {
        m_afEntry[ 0] = afEntry[ 0];
        m_afEntry[ 1] = afEntry[ 4];
        m_afEntry[ 2] = afEntry[ 8];
        m_afEntry[ 3] = afEntry[12];
        m_afEntry[ 4] = afEntry[ 1];
        m_afEntry[ 5] = afEntry[ 5];
        m_afEntry[ 6] = afEntry[ 9];
        m_afEntry[ 7] = afEntry[13];
        m_afEntry[ 8] = afEntry[ 2];
        m_afEntry[ 9] = afEntry[ 6];
        m_afEntry[10] = afEntry[10];
        m_afEntry[11] = afEntry[14];
        m_afEntry[12] = afEntry[ 3];
        m_afEntry[13] = afEntry[ 7];
        m_afEntry[14] = afEntry[11];
        m_afEntry[15] = afEntry[15];
    }
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real> Matrix4<Real>::Inverse () const
{
    Real fA0 = m_afEntry[ 0]*m_afEntry[ 5] - m_afEntry[ 1]*m_afEntry[ 4];
    Real fA1 = m_afEntry[ 0]*m_afEntry[ 6] - m_afEntry[ 2]*m_afEntry[ 4];
    Real fA2 = m_afEntry[ 0]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 4];
    Real fA3 = m_afEntry[ 1]*m_afEntry[ 6] - m_afEntry[ 2]*m_afEntry[ 5];
    Real fA4 = m_afEntry[ 1]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 5];
    Real fA5 = m_afEntry[ 2]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 6];
    Real fB0 = m_afEntry[ 8]*m_afEntry[13] - m_afEntry[ 9]*m_afEntry[12];
    Real fB1 = m_afEntry[ 8]*m_afEntry[14] - m_afEntry[10]*m_afEntry[12];
    Real fB2 = m_afEntry[ 8]*m_afEntry[15] - m_afEntry[11]*m_afEntry[12];
    Real fB3 = m_afEntry[ 9]*m_afEntry[14] - m_afEntry[10]*m_afEntry[13];
    Real fB4 = m_afEntry[ 9]*m_afEntry[15] - m_afEntry[11]*m_afEntry[13];
    Real fB5 = m_afEntry[10]*m_afEntry[15] - m_afEntry[11]*m_afEntry[14];

    Real fDet = fA0*fB5-fA1*fB4+fA2*fB3+fA3*fB2-fA4*fB1+fA5*fB0;
    if ( Math<Real>::FAbs(fDet) <= Math<Real>::EPSILON )
        return Matrix4::ZERO;

    Matrix4 kInv;
    kInv[0][0] = + m_afEntry[ 5]*fB5 - m_afEntry[ 6]*fB4 + m_afEntry[ 7]*fB3;
    kInv[1][0] = - m_afEntry[ 4]*fB5 + m_afEntry[ 6]*fB2 - m_afEntry[ 7]*fB1;
    kInv[2][0] = + m_afEntry[ 4]*fB4 - m_afEntry[ 5]*fB2 + m_afEntry[ 7]*fB0;
    kInv[3][0] = - m_afEntry[ 4]*fB3 + m_afEntry[ 5]*fB1 - m_afEntry[ 6]*fB0;
    kInv[0][1] = - m_afEntry[ 1]*fB5 + m_afEntry[ 2]*fB4 - m_afEntry[ 3]*fB3;
    kInv[1][1] = + m_afEntry[ 0]*fB5 - m_afEntry[ 2]*fB2 + m_afEntry[ 3]*fB1;
    kInv[2][1] = - m_afEntry[ 0]*fB4 + m_afEntry[ 1]*fB2 - m_afEntry[ 3]*fB0;
    kInv[3][1] = + m_afEntry[ 0]*fB3 - m_afEntry[ 1]*fB1 + m_afEntry[ 2]*fB0;
    kInv[0][2] = + m_afEntry[13]*fA5 - m_afEntry[14]*fA4 + m_afEntry[15]*fA3;
    kInv[1][2] = - m_afEntry[12]*fA5 + m_afEntry[14]*fA2 - m_afEntry[15]*fA1;
    kInv[2][2] = + m_afEntry[12]*fA4 - m_afEntry[13]*fA2 + m_afEntry[15]*fA0;
    kInv[3][2] = - m_afEntry[12]*fA3 + m_afEntry[13]*fA1 - m_afEntry[14]*fA0;
    kInv[0][3] = - m_afEntry[ 9]*fA5 + m_afEntry[10]*fA4 - m_afEntry[11]*fA3;
    kInv[1][3] = + m_afEntry[ 8]*fA5 - m_afEntry[10]*fA2 + m_afEntry[11]*fA1;
    kInv[2][3] = - m_afEntry[ 8]*fA4 + m_afEntry[ 9]*fA2 - m_afEntry[11]*fA0;
    kInv[3][3] = + m_afEntry[ 8]*fA3 - m_afEntry[ 9]*fA1 + m_afEntry[10]*fA0;

    Real fInvDet = ((Real)1.0)/fDet;
    for (int iRow = 0; iRow < 4; iRow++)
    {
        for (int iCol = 0; iCol < 4; iCol++)
            kInv[iRow][iCol] *= fInvDet;
    }

    return kInv;
}
//----------------------------------------------------------------------------
template <class Real>
Matrix4<Real> Matrix4<Real>::Adjoint () const
{
    Real fA0 = m_afEntry[ 0]*m_afEntry[ 5] - m_afEntry[ 1]*m_afEntry[ 4];
    Real fA1 = m_afEntry[ 0]*m_afEntry[ 6] - m_afEntry[ 2]*m_afEntry[ 4];
    Real fA2 = m_afEntry[ 0]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 4];
    Real fA3 = m_afEntry[ 1]*m_afEntry[ 6] - m_afEntry[ 2]*m_afEntry[ 5];
    Real fA4 = m_afEntry[ 1]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 5];
    Real fA5 = m_afEntry[ 2]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 6];
    Real fB0 = m_afEntry[ 8]*m_afEntry[13] - m_afEntry[ 9]*m_afEntry[12];
    Real fB1 = m_afEntry[ 8]*m_afEntry[14] - m_afEntry[10]*m_afEntry[12];
    Real fB2 = m_afEntry[ 8]*m_afEntry[15] - m_afEntry[11]*m_afEntry[12];
    Real fB3 = m_afEntry[ 9]*m_afEntry[14] - m_afEntry[10]*m_afEntry[13];
    Real fB4 = m_afEntry[ 9]*m_afEntry[15] - m_afEntry[11]*m_afEntry[13];
    Real fB5 = m_afEntry[10]*m_afEntry[15] - m_afEntry[11]*m_afEntry[14];

    Matrix4 kAdj;
    kAdj[0][0] = + m_afEntry[ 5]*fB5 - m_afEntry[ 6]*fB4 + m_afEntry[ 7]*fB3;
    kAdj[1][0] = - m_afEntry[ 4]*fB5 + m_afEntry[ 6]*fB2 - m_afEntry[ 7]*fB1;
    kAdj[2][0] = + m_afEntry[ 4]*fB4 - m_afEntry[ 5]*fB2 + m_afEntry[ 7]*fB0;
    kAdj[3][0] = - m_afEntry[ 4]*fB3 + m_afEntry[ 5]*fB1 - m_afEntry[ 6]*fB0;
    kAdj[0][1] = - m_afEntry[ 1]*fB5 + m_afEntry[ 2]*fB4 - m_afEntry[ 3]*fB3;
    kAdj[1][1] = + m_afEntry[ 0]*fB5 - m_afEntry[ 2]*fB2 + m_afEntry[ 3]*fB1;
    kAdj[2][1] = - m_afEntry[ 0]*fB4 + m_afEntry[ 1]*fB2 - m_afEntry[ 3]*fB0;
    kAdj[3][1] = + m_afEntry[ 0]*fB3 - m_afEntry[ 1]*fB1 + m_afEntry[ 2]*fB0;
    kAdj[0][2] = + m_afEntry[13]*fA5 - m_afEntry[14]*fA4 + m_afEntry[15]*fA3;
    kAdj[1][2] = - m_afEntry[12]*fA5 + m_afEntry[14]*fA2 - m_afEntry[15]*fA1;
    kAdj[2][2] = + m_afEntry[12]*fA4 - m_afEntry[13]*fA2 + m_afEntry[15]*fA0;
    kAdj[3][2] = - m_afEntry[12]*fA3 + m_afEntry[13]*fA1 - m_afEntry[14]*fA0;
    kAdj[0][3] = - m_afEntry[ 9]*fA5 + m_afEntry[10]*fA4 - m_afEntry[11]*fA3;
    kAdj[1][3] = + m_afEntry[ 8]*fA5 - m_afEntry[10]*fA2 + m_afEntry[11]*fA1;
    kAdj[2][3] = - m_afEntry[ 8]*fA4 + m_afEntry[ 9]*fA2 - m_afEntry[11]*fA0;
    kAdj[3][3] = + m_afEntry[ 8]*fA3 - m_afEntry[ 9]*fA1 + m_afEntry[10]*fA0;

    return kAdj;
}
//----------------------------------------------------------------------------
template <class Real>
Real Matrix4<Real>:eterminant () const
{
    Real fA0 = m_afEntry[ 0]*m_afEntry[ 5] - m_afEntry[ 1]*m_afEntry[ 4];
    Real fA1 = m_afEntry[ 0]*m_afEntry[ 6] - m_afEntry[ 2]*m_afEntry[ 4];
    Real fA2 = m_afEntry[ 0]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 4];
    Real fA3 = m_afEntry[ 1]*m_afEntry[ 6] - m_afEntry[ 2]*m_afEntry[ 5];
    Real fA4 = m_afEntry[ 1]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 5];
    Real fA5 = m_afEntry[ 2]*m_afEntry[ 7] - m_afEntry[ 3]*m_afEntry[ 6];
    Real fB0 = m_afEntry[ 8]*m_afEntry[13] - m_afEntry[ 9]*m_afEntry[12];
    Real fB1 = m_afEntry[ 8]*m_afEntry[14] - m_afEntry[10]*m_afEntry[12];
    Real fB2 = m_afEntry[ 8]*m_afEntry[15] - m_afEntry[11]*m_afEntry[12];
    Real fB3 = m_afEntry[ 9]*m_afEntry[14] - m_afEntry[10]*m_afEntry[13];
    Real fB4 = m_afEntry[ 9]*m_afEntry[15] - m_afEntry[11]*m_afEntry[13];
    Real fB5 = m_afEntry[10]*m_afEntry[15] - m_afEntry[11]*m_afEntry[14];
    Real fDet = fA0*fB5-fA1*fB4+fA2*fB3+fA3*fB2-fA4*fB1+fA5*fB0;
    return fDet;
}
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// explicit instantiation
//----------------------------------------------------------------------------
namespace Wml
{
#ifdef WML_INSTANTIATE_BEFORE
template class WML_ITEM Matrix<4,float>;

#ifdef WML_USING_VC6
template WML_ITEM Matrix<4,float> operator* (float,
    const Matrix<4,float>&);
template WML_ITEM Vector<4,float> operator*
    (const Vector<4,float>&, const Matrix<4,float>&);
#else
template WML_ITEM Matrix<4,float> operator*<4,float> (float,
    const Matrix<4,float>&);
template WML_ITEM Vector<4,float> operator*<4,float>
    (const Vector<4,float>&, const Matrix<4,float>&);
#endif

template<> const Matrix4<float> Matrix4<float>::ZERO(
    0.0f,0.0f,0.0f,0.0f,
    0.0f,0.0f,0.0f,0.0f,
    0.0f,0.0f,0.0f,0.0f,
    0.0f,0.0f,0.0f,0.0f);
template<> const Matrix4<float> Matrix4<float>::IDENTITY(
    1.0f,0.0f,0.0f,0.0f,
    0.0f,1.0f,0.0f,0.0f,
    0.0f,0.0f,1.0f,0.0f,
    0.0f,0.0f,0.0f,1.0f);
template class WML_ITEM Matrix4<float>;

template class WML_ITEM Matrix<4,double>;

#ifdef WML_USING_VC6
template WML_ITEM Matrix<4,double> operator* (double,
    const Matrix<4,double>&);
template WML_ITEM Vector<4,double> operator*
    (const Vector<4,double>&, const Matrix<4,double>&);
#else
template WML_ITEM Matrix<4,double> operator*<4,double> (double,
    const Matrix<4,double>&);
template WML_ITEM Vector<4,double> operator*<4,double>
    (const Vector<4,double>&, const Matrix<4,double>&);
#endif

template<> const Matrix4<double> Matrix4<double>::ZERO(
    0.0,0.0,0.0,0.0,
    0.0,0.0,0.0,0.0,
    0.0,0.0,0.0,0.0,
    0.0,0.0,0.0,0.0);
template<> const Matrix4<double> Matrix4<double>::IDENTITY(
    1.0,0.0,0.0,0.0,
    0.0,1.0,0.0,0.0,
    0.0,0.0,1.0,0.0,
    0.0,0.0,0.0,1.0);
template class WML_ITEM Matrix4<double>;

#else
template class WML_ITEM Matrix<4,float>;

#ifdef WML_USING_VC6
template WML_ITEM Matrix<4,float> operator* (float,
    const Matrix<4,float>&);
template WML_ITEM Vector<4,float> operator*
    (const Vector<4,float>&, const Matrix<4,float>&);
#else
template WML_ITEM Matrix<4,float> operator*<4,float> (float,
    const Matrix<4,float>&);
template WML_ITEM Vector<4,float> operator*<4,float>
    (const Vector<4,float>&, const Matrix<4,float>&);
#endif

template class WML_ITEM Matrix4<float>;
template<> const Matrix4<float> Matrix4<float>::ZERO(
    0.0f,0.0f,0.0f,0.0f,
    0.0f,0.0f,0.0f,0.0f,
    0.0f,0.0f,0.0f,0.0f,
    0.0f,0.0f,0.0f,0.0f);
template<> const Matrix4<float> Matrix4<float>::IDENTITY(
    1.0f,0.0f,0.0f,0.0f,
    0.0f,1.0f,0.0f,0.0f,
    0.0f,0.0f,1.0f,0.0f,
    0.0f,0.0f,0.0f,1.0f);

template class WML_ITEM Matrix<4,double>;

#ifdef WML_USING_VC6
template WML_ITEM Matrix<4,double> operator* (double,
    const Matrix<4,double>&);
template WML_ITEM Vector<4,double> operator*
    (const Vector<4,double>&, const Matrix<4,double>&);
#else
template WML_ITEM Matrix<4,double> operator*<4,double> (double,
    const Matrix<4,double>&);
template WML_ITEM Vector<4,double> operator*<4,double>
    (const Vector<4,double>&, const Matrix<4,double>&);
#endif

template class WML_ITEM Matrix4<double>;
template<> const Matrix4<double> Matrix4<double>::ZERO(
    0.0,0.0,0.0,0.0,
    0.0,0.0,0.0,0.0,
    0.0,0.0,0.0,0.0,
    0.0,0.0,0.0,0.0);
template<> const Matrix4<double> Matrix4<double>::IDENTITY(
    1.0,0.0,0.0,0.0,
    0.0,1.0,0.0,0.0,
    0.0,0.0,1.0,0.0,
    0.0,0.0,0.0,1.0);
#endif
}
//----------------------------------------------------------------------------