Question

c++ 使用Vector类,定义派生类matrix,实现矩阵基本运算

如题，求大神解救QAQ

Answer 1

#include <malloc.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>

// 矩阵的数据结构
struct CMatData
{
	size_t  nWidth;  // 矩阵宽度
	size_t  nHeigth; // 矩阵高度
	long    nRefs;   // 矩阵被引用次数
	BYTE *  GetData() {return (BYTE *)(this + 1);} // 内存地址
};

#define esp  1e-8 // 约等于零的数

template<class T>
class CMat  
{
private:
	size_t    mWidth;  // 矩阵宽度
	size_t    mHeigth; // 矩阵高度
	BYTE   *  mpData;  // 数据地址
	T      ** mpRows;  // 矩阵行指针
	// 申请分配内存
	void AllocBuffer(size_t w, size_t h)
	{
		if(w * h == 0) return;
		else
		{
			size_t rowSize = ((sizeof(T) * w + 3) & ~3); // 对齐到4BYTE
			size_t allocSize = sizeof(CMatData) + rowSize * h; 
			CMatData * p = (CMatData *)malloc(allocSize); 
			mpRows       = (T **)malloc(sizeof(T *) * h);
			if(p == NULL || mpRows == NULL) // 申请失败
			{
				Init();
				if(p) free(p);
				if(mpRows) free(mpRows);
			}
			else
			{
				p->nRefs   = 1;
				mWidth  = p->nWidth  = w;
				mHeigth = p->nHeigth = h;
				mpData  = p->GetData();
				for(size_t i=0; i<mHeigth; i++) mpRows[i] = (T *)(mpData + rowSize * i);
			}
		}
	}
	// 释放内存
	static void FreeData(void * p) {free(p);}
	// 引用释放
	void Release()
	{
		if(mpData == NULL) return;
		CMatData * p = GetData();
		assert(p->nRefs != 0);
		// 引用减少1，当引用次数等于0，释放内存
		if(--p->nRefs <= 0) FreeData(p);
		free(mpRows);
		Init();
	}
	// 引用释放
	static void Release(CMatData * p)
	{
		if(p)
		{
			assert(p->nRefs != 0);
			if(--p->nRefs <= 0) FreeData(p);
		}
	}
	// 释放旧数据并申请新的内存空间
	void AllocBeforeWrite(size_t w, size_t h)
	{
		if(mpData)
		{
			CMatData * p = GetData();
			if(p->nRefs > 1 || p->nHeigth != h || p->nWidth != w)
			{
				Release();
				Init();
			}
		}
		if(mpData == NULL) AllocBuffer(w, h);
		assert(GetData()->nRefs == 1);
	}
	// 更改数据前脱离引用
	void CopyBeforWrite()
	{
		if(mpData == NULL) return;
		if(GetData()->nRefs <= 1) return;
		CMat tmp(*this);
		Release();
		AssignCopy(tmp);
	}
                 
    template<class T1, class T2>
		T2 * Copy(T1 * f, T1 * t, T2 * o) const
	{
		while(f < t) *o++ = *f++;
		return o;
	}
    template<class T1, class V>
		T1 * Fill(T1 * f, T1 * t, V v) const
	{
		while(f < t) *f++ = v;
		return f;
	}
 
	// 检测数据有效性
	inline void Youxiao(size_t w, size_t h) const
	{
		if(mpData == NULL || w >= mWidth || h >= mHeigth) assert(0);
	}
	inline void Youxiao(size_t w) const
	{
		if(mpData == NULL || w >= mWidth) assert(0);
	}
    // 复制矩阵
	void AssignCopy(const CMat & r)
	{
		assert(mpData == NULL);
		if(r.IsEmpty()) return; // 直接返回
		AllocBuffer(r.mWidth, r.mHeigth); // 分配内存
		for(size_t i=0; i<mHeigth; i++) // 复制数据
		{
			const T * p = r.GetRows(i);
			Copy(p, p + mWidth, GetRows(i));
		}
	}
	// 引用矩阵
	void RefCopy(const CMat & r)
	{
		assert(mpData == NULL);
		if(r.IsEmpty()) return; // 直接返回
		mpRows = (T **)malloc(sizeof(T *) * r.mHeigth);
		if(mpRows == NULL) return; // 内存分配错误
		
		mpData   = r.mpData; // 直接复制内存地址
		mWidth   = r.mWidth;
		mHeigth  = r.mHeigth;
		T ** p = r.mpRows;
		Copy(p, p + r.mHeigth, mpRows);
		++GetData()->nRefs; // 引用计数增加1
	}
 
protected:
	// 初始化
	void Init()
	{
		mpData   = NULL;
		mpRows   = NULL;
		mWidth   = 0;
		mHeigth  = 0;
	}
	// 取得数据实际地址
	CMatData * GetData() const
	{
		if(mpData) return (CMatData *)mpData - 1;
		return NULL;
	}
	// 交换a b两行
	inline void SwapRows(size_t a, size_t b)
	{
		Youxiao(a);
		Youxiao(b);
		CopyBeforWrite();
		T * tmp = mpRows[a];
		mpRows[a] = mpRows[b];
		mpRows[b] = tmp;
	}
	// row[a] += row[b]; b行加a行，结果保存到a行
	inline void AddRows(size_t a, size_t b)
	{
		Youxiao(a);
		Youxiao(b);
		CopyBeforWrite();
		T * pal = mpRows[a];
		T * pah = pal + mWidth;
		T * pbl = mpRows[b];
		while(pal < pah) *pal++ += *pbl++;
    }
	// row[a] -= row[b]; a行减b行，结果保存到a行
	inline void SubRows(size_t a, size_t b)
	{
		Youxiao(a);
		Youxiao(b);
		CopyBeforWrite();
		T * pal = mpRows[a];
		T * pah = pal + mWidth;
		T * pbl = mpRows[b];
		while(pal < pah) *pal++ -= *pbl++;
	}
	// row[a] -= row[b] * k; a行减b行乘以常数k，结果保存到a行
	inline void SubRows(size_t a, size_t b, T k)
	{
		Youxiao(a);
		Youxiao(b);
		CopyBeforWrite();
		T * pal = mpRows[a];
		T * pah = pal + mWidth;
		T * pbl = mpRows[b];
		while(pal < pah) *pal++ -= *pbl++ * k;
	}
	// row[a] *= k;  a行乘以常数k
	inline void MulRows(size_t a, T k)
	{
		Youxiao(a);
		CopyBeforWrite();
		T * pal = mpRows[a];
		T * pah = pal + mWidth;
		while(pal < pah) *pal++ *= k;
    }
	// row[a] /= k; a行除以常数k
	inline void DivRows(size_t a, T k)
	{
		Youxiao(a);
		CopyBeforWrite();
		k = T(1) / k;
		MulRows(a, k);
	} 
public:
	// 构造、析构函数

	CMat() {Init();}
	CMat(size_t w, size_t h)
	{
		Init();
		AllocBuffer(w, h);
		SetData();
	}
	CMat(size_t w, size_t h, const T * data)
	{
		Init();
		AllocBuffer(w, h);
		SetData(data);
	}
	CMat(const CMat & r) // 复制构造函数
	{
		Init();
		if(!r.IsEmpty()) 
		{
			if(r.GetData()->nRefs < 1) AssignCopy(r); // 源数据被锁定，需要复制数据
			else RefCopy(r);
		}
	}
	virtual ~CMat() // 
	{
		CMatData * p = GetData();
		if(p) if(--p->nRefs <= 0) FreeData(p);
		if(mpRows) free(mpRows);
	}

	// 复制函数
	const CMat & operator=(const CMat & r)
	{
		CMatData * p1 = GetData();
		CMatData * p2 = r.GetData();
		if(p1 == p2)
		{
			mWidth  = r.mWidth;
			mHeigth = r.mHeigth;
			free(mpRows);
			mpRows = new T*[mHeigth];
			for(int i = 0; i < mHeigth; ++i) mpRows[i] = r.mpRows[i];
			return *this;
		}
		Release();
		if(!r.IsEmpty()) 
		{
			if(r.GetData()->nRefs < 1) AssignCopy(r); // 源数据被锁定，需要复制数据
			else RefCopy(r);
		}
		return *this;
	}

	// 矩阵运算
	
	// 矩阵转置
	inline CMat Transpose() const
	{
		CMat t(mHeigth, mWidth); // 建立一个新矩阵
		size_t i, j;
		for(j=0; j<mHeigth; j++)
		{
			const T * p = (const T *)mpRows[j];
			for(i=0; i<mWidth; i++) t.At(j, i) = *p++;
		}
		return t;
	}
	// 矩阵加法
	inline CMat Add(const CMat & l) const
	{
		if(l.mHeigth != mHeigth || l.mWidth != mWidth)
		{
			printf("CMat::Add error 矩阵参数不一致\n");
			return CMat(); // 参数不一致返回空矩阵
		}
		CMat t(mWidth, mHeigth);
		size_t i;
		for(i=0; i<mHeigth; i++)
		{
			const T * prl = (const T *)mpRows[i];
			const T * prh = prl + mWidth;
			const T * pl  = (const T *)l.mpRows[i];
			T * pt        = t.mpRows[i];
			while(prl < prh) *pt++ = *prl++ + *pl++;
		}
		return t;
	}
	// 矩阵减法
	inline CMat Sub(const CMat & l) const
	{
		if(l.mHeigth != mHeigth || l.mWidth != mWidth)
		{
			printf("CMat::Add error 矩阵参数不一致\n");
			return CMat(); // 参数不一致返回空矩阵
		}
		CMat t(mWidth, mHeigth);
		size_t i;
		for(i=0; i<mHeigth; i++)
		{
			const T * prl = (const T *)mpRows[i];
			const T * prh = prl + mWidth;
			const T * pl  = (const T *)l.mpRows[i];
			T * pt        = t.mpRows[i];
			while(prl < prh) *pt++ = *prl++ - *pl++;
		}
		return t;
	}
	// 矩阵乘以常数
	inline CMat Mul(const T & k) const
	{
		CMat t(mWidth, mHeigth);
		size_t i;
		for(i=0; i<mHeigth; i++)
		{
			const T * prl = (const T *)mpRows[i];
			const T * prh = prl + mWidth;
			T       * pt  = t.mpRows[i];
			while(prl < prh) *pt++ = *prl++ * k;
		}
		return t;
	}
	// 矩阵除以常数
	inline CMat Div(const T & k) const
	{
		T d = T(1) / k;
		return Mul(d); // 除以k就相当于乘以 1/k
	}
	// 两行的数量积
	inline T DotRows(size_t a, size_t b) const
	{
		T r = 0;
		T * pal = mpRows[a];
		T * pah = pal + mWidth;
		T * pbl = mpRows[b];
		while(pal < pah) r += *pal++ * *pbl++;
		return r;
	}
	// 两个矩阵乘法
	inline CMat Mul(const CMat & l) const
	{
		if(mWidth != l.mHeigth)
		{
			printf("CMat::Mul error 矩阵参数不一致\n");
			return CMat(); // 参数不一致返回空矩阵
		}
		size_t i, j, k;
		CMat t(l.mWidth, mHeigth);
		const T** plRows = (const T**)l.mpRows;
		for(i=0; i<mHeigth; i++)
		{
			T       * ptl = t.mpRows[i];
			const T * pr  = (const T *)mpRows[i];
			for(j=0; j<l.mWidth; j++)
			{
				const T *prl = pr;
				for(k=0; k<mWidth; k++) *ptl += *prl++ * l.At(j, k);
				ptl++;
			}
		}
		return t;
	}
	// 矩阵乘以单行数据
	inline T * Mul(T * p) const
	{
		CMat l(1, mWidth, p);
		CMat r = Mul(l);
		return r.GetList(p, mWidth);
	}
	// 高斯消元，用于求解线性方程和求矩阵的逆矩阵
	inline CMat GaussianElimination() const
	{
		if(mWidth <= mHeigth) return CMat();
		bool loop = true;
		CMat r(*this);
		size_t i, j, k;
		for(i=0; i<r.mHeigth && loop; i++)
		{
			size_t pos = i;
			T max = r.At(i, i);
			for(k = i + 1; k < r.mHeigth; k++) // 找出最大的数值，减小浮点数计算误差
			{
				T y = r.At(k, i);
				if(y > max) {pos = k; max = y;}
			}
			if(i != pos) r.SwapRows(i, pos); // 交换i行与找到的最大数据行
			if(r.Iszero(i, i)) {loop = false; break;} // 如果最大的数值为零，无解
 
			// 使i行第i个元素等于1
			r.DivRows(i, r.At(i, i));
			// 主计算循环 j行减去i行乘以j行的第i个元素, 使j行的第i个元素为零
			for(j = 0; j < r.mHeigth; j++) if(i != j) r.SubRows(j, i, r.At(i, j));
		}
		if(loop == false) return CMat();
		r.Resize(r.mWidth - r.mHeigth, r.mHeigth, r.mHeigth, 0);
		return r;
	}
	// 转换为阶梯形矩阵
	inline CMat Echelon() const
	{
		if(mpData == NULL) return CMat();
		CMat r(*this);
		bool loop = true;
		size_t i, j, k;
		for(i=0; i<r.mHeigth && loop; i++)
		{
			size_t pos = i;
			T max = r.At(i, i);
			for(k = i + 1; k < r.mHeigth; k++) // 找出最大的数值，减小浮点数计算误差
			{
				T y = r.At(k, i);
				if(y > max) {pos = k; max = y;}
			}
			if(i != pos) r.SwapRows(i, pos); // 交换i行与找到的最大数据行
			// 主计算循环
			for(j = i + 1; j < r.mHeigth; j++) r.SubRows(j, i, r.At(i, j) / r.At(i, i));
        }
        if(!loop) r.Resize(r.mWidth, i, 0, 0);
        return r;
    }
    // 求行列式的值，为减少计算，1，2，3系列的直接计算
    inline T Det() const
	{
		if(mWidth != mHeigth) return 0;
		switch(mWidth)
		{
		case 1:
			return At(0,0);
		case 2:
			return At(0,0) * At(1,1) - At(0,1) * At(1, 0);
		case 3:
			return At(0,0) * At(1,1) * At(2,2) + At(0,1) * At(1,2) * At(2,0) + At(0,2) * At(1,0) * At(2,1)
				- (At(2,0) * At(1,1) * At(0,2) + At(2,1) * At(1,2) * At(0,0) + At(2,2) * At(1,0) * At(0,1));
        default:
			{
				CMat r = Echelon();
				if(r.mpData == NULL)  return 0;
				if(r.mWidth != r.mHeigth) return 0;
				size_t i;
				T v = r.At(0,0);
				for(i=1; i<mWidth; i++) v *= r.At(i, i);
				return v;
			}
		}
	}
	// 求矩阵的逆矩阵
	inline CMat Inv() const
	{
		if(mWidth != mHeigth) return CMat();
		CMat r(2 * mWidth, mWidth); // 扩展矩阵
		size_t i;
		for(i=0; i<mWidth; i++) // 复制数据到扩展矩阵
		{
			const T * psl = (const T *)mpRows[i];
			const T * psh = psl + mWidth;
			T       * pd  = r.mpRows[i];
			Copy(psl, psh, pd);
			r.At(mWidth + i, i) = 1;
		}
		return r.GaussianElimination(); // 返回扩展矩阵的解
	}
    // 矩阵的一个应用：最小二乘法
    static T * LeastSquares(T * r, const T * y, const T * x, size_t num, size_t n = 1);
 
	// 属性

	// 获取分配的矩阵宽度
	size_t GetAllocWidth() const;
	// 获取分配的矩阵高度
	size_t GetAllocHeigth() const;
	// 获取矩阵宽度
	size_t GetWidth() const;
	// 获取矩阵高度
	size_t GetHeigth() const;
	// 判断矩阵是否为空
	bool IsEmpty() const;
	// 取得h行的地址
    inline T * GetRows(int h);
	// 取得矩阵的h行，w列的引用
    inline T & At(int w, int h);
	// 取得h行的地址
	inline const T * GetRows(int h) const;
	// 取得矩阵的h行，w列的引用
    inline const T & At(int w, int h) const;
    // 取得矩阵数据
	inline T * GetList(T * p, size_t max) const;
	// 设置数据
	int SetData(const T * d = NULL)
	{
		assert(mpData);
		assert(mpRows);
		if(d) for(size_t i = 0; i < mHeigth; i++)
		{
			Copy(d, d + mWidth, mpRows[i]);
			d += mWidth;
		}
		else for(size_t i = 0; i < mHeigth; i++) Fill(mpRows[i], mpRows[i] + mWidth, T(0));
		return mHeigth * mWidth;
	}
	// 改变矩阵大小
	inline void Resize(size_t w, size_t h, size_t xOff, size_t yOff)
	{
		Youxiao(w + xOff - 1, h + yOff - 1);
		size_t i;
		mWidth  = w;
		mHeigth = h;
		for(i=0; i<mHeigth; i++) mpRows[i] = mpRows[i + yOff] + xOff;
	}

    inline static bool Eqzero(const T & k) 
    {
        if((k < T(0) ? -k : k) < esp) return true;
        return false;
    }
    inline bool Iszero(size_t w, size_t h) const
    {
        return Eqzero(At(w, h));
    }
};

Answer 2

microsoft 的D3DX9库什么矩阵都有的随便用。
我有一些代码发给你。
D3DXMATRIX* WINAPI D3DXMatrixInit(D3DXMATRIX* pOut,
float m11, float m12, float m13, float m14,
float m21, float m22, float m23, float m24,
float m31, float m32, float m33, float m34,
float m41, float m42, float m43, float m44)
{
pOut->_11 = m11; pOut->_12 = m12; pOut->_13 = m13; pOut->_14 = m14;
pOut->_21 = m21; pOut->_22 = m22; pOut->_23 = m23; pOut->_24 = m24;
pOut->_31 = m31; pOut->_32 = m32; pOut->_33 = m33; pOut->_34 = m34;
pOut->_41 = m41; pOut->_42 = m42; pOut->_43 = m43; pOut->_44 = m44;
return pOut;
}
这个就是矩阵看的懂英文的就不用我说了。
D3DXMATRIX* WINAPI D3DXMatrixMultiply
( D3DXMATRIX *pOut, CONST D3DXMATRIX *pM1, CONST D3DXMATRIX *pM2 )
{
ATLASSERT(pOut!=NULL&&pM1!=NULL&&pM2!=NULL);
D3DXMATRIX matTemp;
float fSum = 0;
for (int nRow=0; nRow<4; nRow++)
{
for (int nCol=0; nCol<4; nCol++)
{
fSum=0;
for (int nIndex=0;nIndex<4; nIndex++)
{
fSum+=pM1->m[nRow][nIndex]*pM2->m[nIndex][nCol];
}
matTemp.m[nRow][nCol] = fSum;
}
}
pOut[0]=matTemp;
return pOut;
}
D3DXMATRIX* WINAPI D3DXMatrixScaling
( D3DXMATRIX *pOut, FLOAT sx, FLOAT sy, FLOAT sz )
{
D3DXMatrixIdentity(pOut);
pOut->_11=sx;
pOut->_22=sy;
pOut->_33=sz;
return pOut;
}
D3DXMATRIX* WINAPI D3DXMatrixTranslation
( D3DXMATRIX *pOut, FLOAT x, FLOAT y, FLOAT z )
{
D3DXMatrixIdentity(pOut);
pOut->_41=x;
pOut->_42=y;
pOut->_43=z;
return pOut;
}
D3DXMATRIX* WINAPI D3DXMatrixRotationX
( D3DXMATRIX *pOut, FLOAT Angle )
{
FLOAT fCos=cos(Angle);
FLOAT fSin=sin(Angle);
D3DXMatrixInit(pOut,1,0,0,0,
0,fCos,fSin,0,
0,-fSin,fCos,0,
0,0,0,1);
return pOut;
}
D3DXMATRIX* WINAPI D3DXMatrixRotationY
( D3DXMATRIX *pOut, FLOAT Angle )
{
FLOAT fCos=cos(Angle);
FLOAT fSin=sin(Angle);
D3DXMatrixInit(pOut,
fCos,0,-fSin,0,
0,1,0,0,
fSin,0,fCos,0,
0,0,0,1);
return pOut;
}
D3DXMATRIX* WINAPI D3DXMatrixRotationZ
( D3DXMATRIX *pOut, FLOAT Angle )
{
FLOAT fCos=cos(Angle);
FLOAT fSin=sin(Angle);
D3DXMatrixInit(pOut,
fCos,0,fSin,0,
-fSin,fCos,0,0,
0,0,1,0,
0,0,0,1);
return pOut;
}
D3DXMATRIX* WINAPI D3DXMatrixRotationYawPitchRoll
( D3DXMATRIX *pOut, FLOAT Yaw, FLOAT Pitch, FLOAT Roll )
{
D3DXMATRIX matTemp;
D3DXMatrixRotationY(pOut,Yaw);
D3DXMatrixRotationX(&matTemp,Pitch);
D3DXMatrixMultiply(pOut,pOut,&matTemp);
D3DXMatrixRotationZ(&matTemp,Roll);
D3DXMatrixMultiply(pOut,pOut,&matTemp);
return pOut;
}
D3DXVECTOR4* WINAPI D3DXVec4Transform
( D3DXVECTOR4 *pOut, CONST D3DXVECTOR4 *pV, CONST D3DXMATRIX *pM )
{
FLOAT fSum = 0;
ATLASSERT(pOut!=pV);
for (int nCol=0; nCol < 4; nCol++)
{
fSum=0;
for (int nRow=0; nRow<4; nRow++)
{
fSum+=((FLOAT*)pV)[nRow]*pM->m[nRow][nCol];
}
((FLOAT*)pOut)[nCol] = fSum;
}
return pOut;
}
D3DXVECTOR4* WINAPI D3DXVec4TransformArray
( D3DXVECTOR4 *pOut, UINT OutStride,
CONST D3DXVECTOR4 *pV, UINT VStride,
CONST D3DXMATRIX *pM, UINT n )
{
ATLASSERT(pOut!=pV);
我没有完成的；
return pOut;
}