矩阵运算概要

最新推荐文章于 2023-11-27 00:22:53 发布
最新推荐文章于 2023-11-27 00:22:53 发布
阅读量730 收藏
点赞数 2
分类专栏： 数学 文章标签： 矩阵运算 矩阵乘 行列式 逆矩阵 特征值
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/linolzhang/article/details/55549018
版权

       提到矩阵，相信大家都不陌生，作为最基本的数学概念之一，我们 知道常用的矩阵运算包括：矩阵加减乘、行列式、转置矩阵、逆矩阵、特征向量（特征值）。

具体概念不做过多讲讲，直接上代码（相信代码是最好的讲述方式）。

/* 矩阵定义 - 4*4齐次矩阵 用于openGL
   linolzhang, 2010.2.10
*/

#ifndef MATRIX4_H
#define MATRIX4_H

#include "Tuple4.h"
#include "BaseOperation.h"


template <typename T>
class Matrix4
{
public:
	T  m[16]; // 一维数组，与openGL矩阵一致 - 按列优先

	// [ m0, m4, m8, m12 ]
	// [ m1, m5, m9, m13 ]
	// [ m2, m6, m10,m14 ]
	// [ m3, m4, m11,m15 ]

public:
	Matrix4() // 默认构造函数，设置单位矩阵
	{
		setIdentity();
	}

	Matrix4( T  m0, T  m1, T  m2, T  m3,
			 T  m4, T  m5, T  m6, T  m7,
			 T  m8, T  m9, T m10, T m11,
			 T m12, T m13, T m14, T m15 ) // 构造函数，指定值
	{
			m[ 0] =  m0; m[ 1] =  m1; m[ 2] =  m2; m[ 3] =  m3;
			m[ 4] =  m4; m[ 5] =  m5; m[ 6] =  m6; m[ 7] =  m7;
			m[ 8] =  m8; m[ 9] =  m9; m[10] = m10; m[11] = m11;
			m[12] = m12; m[13] = m13; m[14] = m14; m[15] = m15;
	}

	Matrix4(const T* elements) // 构造函数，指针指定
	{
		memcpy(m, elements, sizeof(T)*16);
	}

	Matrix4(const Matrix4 &mat)  // 拷贝构造函数
	{
		memcpy(m,mat.m, sizeof(T)*16);
	}
	const Matrix4 &operator  = (const Matrix4 &mat) // 赋值操作符
	{
		memcpy(m,mat.m, sizeof(T)*16);
		return *this;
	}

	// --------------------------------------------------------
	// 构造函数 - 由其它转换
	Matrix4(const Tuple3<T> &rAxis, float rAngle)  // 由 旋转轴向量，旋转角 得到
	{
		float sinAngle = fastSin(rAngle),
			  cosAngle = fastCos(rAngle),
			  oneMinusCosAngle = 1.0f - cosAngle;

		m[ 0] = (rAxis.x)*(rAxis.x) + cosAngle*(1-(rAxis.x)*(rAxis.x));
		m[ 4] = (rAxis.x)*(rAxis.y)*(oneMinusCosAngle) - sinAngle*rAxis.z;
		m[ 8] = (rAxis.x)*(rAxis.z)*(oneMinusCosAngle) + sinAngle*rAxis.y;

		m[ 1] = (rAxis.x)*(rAxis.y)*(oneMinusCosAngle) + sinAngle*rAxis.z;
		m[ 5] = (rAxis.y)*(rAxis.y) + cosAngle*(1-(rAxis.y)*(rAxis.y));
		m[ 9] = (rAxis.y)*(rAxis.z)*(oneMinusCosAngle) - sinAngle*rAxis.x;

		m[ 2] = (rAxis.x)*(rAxis.z)*(oneMinusCosAngle) - sinAngle*rAxis.y;
		m[ 6] = (rAxis.y)*(rAxis.z)*(oneMinusCosAngle) + sinAngle*rAxis.x;
		m[10] = (rAxis.z)*(rAxis.z) + cosAngle*(1-(rAxis.z)*(rAxis.z));

		m[3] = m[7] = m[11] = 0;
		m[12] = m[13] = m[14] = 0;
		m[15] = 1;
	}

	// --------------------------------------------------------
	// 重载 +=, -=, *=矩阵, *=常量
	void operator +=(const Matrix4 &mat)
	{
		for(int i = 0; i < 16; i++)
			m[i] += mat.m[i];
	}

	void operator -=(const Matrix4 &mat)
	{
		for(int i = 0; i < 16; i++)
			m[i] -= mat.m[i];
	}
	void operator *= (const Matrix4 &mat)
	{
		setElements( m[0]*mat.m[ 0]+m[4]*mat.m[ 1]+m[ 8]*mat.m[ 2]+ m[12]*mat.m[ 3],
					 m[1]*mat.m[ 0]+m[5]*mat.m[ 1]+m[ 9]*mat.m[ 2]+ m[13]*mat.m[ 3],
					 m[2]*mat.m[ 0]+m[6]*mat.m[ 1]+m[10]*mat.m[ 2]+ m[14]*mat.m[ 3],
					 m[3]*mat.m[ 0]+m[7]*mat.m[ 1]+m[11]*mat.m[ 2]+ m[15]*mat.m[ 3],
					 m[0]*mat.m[ 4]+m[4]*mat.m[ 5]+m[ 8]*mat.m[ 6]+ m[12]*mat.m[ 7],
					 m[1]*mat.m[ 4]+m[5]*mat.m[ 5]+m[ 9]*mat.m[ 6]+ m[13]*mat.m[ 7],
					 m[2]*mat.m[ 4]+m[6]*mat.m[ 5]+m[10]*mat.m[ 6]+ m[14]*mat.m[ 7],
					 m[3]*mat.m[ 4]+m[7]*mat.m[ 5]+m[11]*mat.m[ 6]+ m[15]*mat.m[ 7],
					 m[0]*mat.m[ 8]+m[4]*mat.m[ 9]+m[ 8]*mat.m[10]+ m[12]*mat.m[11],
					 m[1]*mat.m[ 8]+m[5]*mat.m[ 9]+m[ 9]*mat.m[10]+ m[13]*mat.m[11],
					 m[2]*mat.m[ 8]+m[6]*mat.m[ 9]+m[10]*mat.m[10]+ m[14]*mat.m[11],
					 m[3]*mat.m[ 8]+m[7]*mat.m[ 9]+m[11]*mat.m[10]+ m[15]*mat.m[11],
					 m[0]*mat.m[12]+m[4]*mat.m[13]+m[ 8]*mat.m[14]+ m[12]*mat.m[15],
					 m[1]*mat.m[12]+m[5]*mat.m[13]+m[ 9]*mat.m[14]+ m[13]*mat.m[15],
					 m[2]*mat.m[12]+m[6]*mat.m[13]+m[10]*mat.m[14]+ m[14]*mat.m[15],
					 m[3]*mat.m[12]+m[7]*mat.m[13]+m[11]*mat.m[14]+ m[15]*mat.m[15] );
	} 
	void operator *= (T scale) // 乘常量
	{
		for(int i=0; i<16; i++)
			m[i] *= scale;
	}

	// ----------------------------------
	// ==, !=
	bool operator == (const Matrix4 &mat)
	{
		return memcmp(m, mat.m, sizeof(T)*16) == 0;
	}

	bool operator != (const Matrix4 &mat)
	{
		return memcmp(m, mat.m, sizeof(T)*16) != 0;
	}

	
	// --------------------------------------------------------
	// 矩阵求解
	const T getDeterminant()  // 得到行列式的值（行列式 !=0 可逆）
	{	
		T det;
		det =  m[0] * m[5] * m[10];
		det += m[4] * m[9] * m[2];
		det += m[8] * m[1] * m[6];
		det -= m[8] * m[5] * m[2];
		det -= m[4] * m[1] * m[10];
		det -= m[0] * m[9] * m[6];
		return det;
	}

	
	// --------------------------------------------------------
	// 矩阵变换
	void setElements(const T* elements) // 重新设置矩阵元素
	{
		memcpy(m, elements, sizeof(T)*16);
	}

	void setElements( T  m0, T  m1, T  m2, T  m3,
					  T  m4, T  m5, T  m6, T  m7,
					  T  m8, T  m9, T m10, T m11,
					  T m12, T m13, T m14, T m15 ) // 重新设置矩阵元素
	{
		m[ 0] =  m0; m[ 1] =  m1; m[ 2] =  m2; m[ 3] =  m3;
		m[ 4] =  m4; m[ 5] =  m5; m[ 6] =  m6; m[ 7] =  m7;
		m[ 8] =  m8; m[ 9] =  m9; m[10] = m10; m[11] = m11;
		m[12] = m12; m[13] = m13; m[14] = m14; m[15] = m15;
	}

	void setElements(const Matrix4  &matrix) // 重新设置矩阵元素
	{
		setElements(matrix.m);
	}

	// ----------------------------------
	void setZero() // 设置为0矩阵
	{
		for(int i = 0; i < 16; i++)
			m[i] = 0;
	}

	void setIdentity() // 设置为单位矩阵
	{
		m[ 0] = 1; m[ 1] = 0; m[ 2] = 0; m[ 3] = 0;
		m[ 4] = 0; m[ 5] = 1; m[ 6] = 0; m[ 7] = 0;
		m[ 8] = 0; m[ 9] = 0; m[10] = 1; m[11] = 0;
		m[12] = 0; m[13] = 0; m[14] = 0; m[15] = 1;
	}

	// ----------------------------------
	void setTranspose() // 矩阵转置
	{
		T temp = 0;
		temp  = m[4];  m[4]  = m[1];
		m[1]  = temp;  temp  = m[8];
		m[8]  = m[2];  m[2]  = temp;
		temp  = m[12]; m[12] = m[3];
		m[3]  = temp;  temp  = m[9];
		m[9]  = m[6];  m[6]  = temp;
		temp  = m[13]; m[13] = m[7];
		m[7]  = temp;  temp  = m[14];
		m[14] = m[11]; m[11] = temp;
	}

	void setNegate() // 矩阵取反
	{
		for(int i = 0; i < 16; i++)
			m[i] =-m[i];
	}

	// --------------------------------------------------------
	bool setInverse() // 求逆矩阵
	{
		int i, j, k, swap;
		T temp[16], t;

		memcpy(temp, m, 16*sizeof(T)); // 当前矩阵拷贝到暂存temp
		setIdentity();

		for (i = 0; i < 4; i++)
		{
			swap = i;
			for (j = i + 1; j < 4; j++)
			{
				if (fabs(temp[j*4 + i]) >  fabs(temp[i*4 + i]))
				{
					swap = j;
				}
			}

			if (swap != i)
			{
				for (k = 0; k < 4; k++) 
				{
					t                = temp[i*4 + k];
					temp[i*4 + k]    = temp[swap*4 + k];
					temp[swap*4 + k] = t;

					t             = m[i*4 + k];
					m[i*4 + k]    = m[swap*4 + k];
					m[swap*4 + k] = t;
				}
			}

			if(!temp[i*4 + i])
				return false;

			t = temp[i*4 + i];

			for (k = 0; k < 4; k++) 
			{
				temp[i*4 + k] /= t;
				m[i*4 + k]     = m[i*4 + k] / t;
			}

			for (j = 0; j < 4; j++) 
			{
				if (j != i) 
				{
					t = temp[j*4 + i];
					for (k = 0; k < 4; k++) 
					{
						temp[j*4 + k] -= temp[i*4 + k] * t;
						m[j*4 + k]     = m[j*4 + k] - m[i*4 + k] * t;
					}
				}
			}
		}
		return true;
	}

	
	// --------------------------------------------------------
	// 矩阵转化，由其它表达方式转换得到
	void setTranslate(const Tuple3<T> &t) // 设置为平移矩阵 - 由三元组指定平移量
	{
		setIdentity();
		m[12] =  t.x;
		m[13] =  t.y;
		m[14] =  t.z;
	}

	void setTranslate(T x,T y,T z) // 设置为平移矩阵 - 由x,y,z指定平移量
	{
		setIdentity();
		m[12] =  x;
		m[13] =  y;
		m[14] =  z;
	}

	void setRotationX(float xAngle) // 设置为旋转矩阵 - 指定绕 X 旋转角度 - 欧拉角
	{
		setIdentity();
		m[ 5] =  fastCos(xAngle);
		m[ 6] =  fastSin(xAngle);
		m[ 9] = -m[6];
		m[10] =  m[5];
	}

	void setRotationY(float yAngle) // 设置为旋转矩阵 - 指定绕 Y 旋转角度 - 欧拉角
	{
		setIdentity();
		m[ 0] =  fastCos(yAngle);
		m[ 2] =  fastSin(yAngle);
		m[ 8] = -m[2];
		m[10] =  m[0];
	}  

	void setRotationZ(float zAngle)  // 设置为旋转矩阵 - 指定绕 Z 旋转角度 - 欧拉角
	{
		setIdentity();
		m[0] =  fastCos(zAngle);
		m[1] =  fastSin(zAngle);
		m[4] = -m[1];
		m[5] =  m[0];
	}

	void setRotationXYZ(float x,float y,float z) // 设置为旋转矩阵 - 分别指定绕 X,Y,Z 旋转角度 - 欧拉角
	{
		float cosX = fastCos(x), sinX = fastSin(x),
			  cosY = fastCos(y), sinY = fastSin(y),
			  cosZ = fastCos(z), sinZ = fastSin(z);

		setElements( cosY * cosZ + sinX * sinY * sinZ,   -cosX * sinZ,    sinX * cosY * sinZ - sinY * cosZ,  0,
					 cosY * sinZ - sinX * sinY * cosZ,    cosX * cosZ,   -sinY * sinZ - sinX * cosY * cosZ,  0,
					 cosX * sinY,                         sinX,           cosX * cosY,                       0,
					 0,                                   0,              0,                                 1  );
	}
	void setRotateAxis(const Tuple3<T> &axis, float angle) // 设置为旋转矩阵 - 由旋转轴向量，旋转角指定
	{
		float sinAngle = fastSin(angle),
			  cosAngle = fastCos(angle),
			  oneMinusCosAngle = 1.0f - cosAngle;

		m[ 0] = (axis.x)*(axis.x) + cosAngle*(1-(axis.x)*(axis.x));
		m[ 4] = (axis.x)*(axis.y)*(oneMinusCosAngle) - sinAngle*axis.z;
		m[ 8] = (axis.x)*(axis.z)*(oneMinusCosAngle) + sinAngle*axis.y;

		m[ 1] = (axis.x)*(axis.y)*(oneMinusCosAngle) + sinAngle*axis.z;
		m[ 5] = (axis.y)*(axis.y) + cosAngle*(1-(axis.y)*(axis.y));
		m[ 9] = (axis.y)*(axis.z)*(oneMinusCosAngle) - sinAngle*axis.x;

		m[ 2] = (axis.x)*(axis.z)*(oneMinusCosAngle) - sinAngle*axis.y;
		m[ 6] = (axis.y)*(axis.z)*(oneMinusCosAngle) + sinAngle*axis.x;
		m[10] = (axis.z)*(axis.z) + cosAngle*(1-(axis.z)*(axis.z));

		m[3] = m[7] = m[11] = 0;
		m[12] = m[13] = m[14] = 0;
		m[15] = 1;
	}

	void setScale(T x,T y,T z)  // 设置为缩放矩阵，由 x,y,z 指定缩放量
	{
		setIdentity();
		m[ 0] = x;
		m[ 5] = y;
		m[10] = z;
	}


// --------------------------------------------------------
// 定义矩阵运算 - 左乘
	void multByMatrix(const Matrix4<T> &mat)
	{
		setElements(    mat.m[0]*m[ 0] + mat.m[4]*m[ 1] + mat.m[ 8]*m[ 2] + mat.m[12]*m[ 3],
						mat.m[1]*m[ 0] + mat.m[5]*m[ 1] + mat.m[ 9]*m[ 2] + mat.m[13]*m[ 3],
						mat.m[2]*m[ 0] + mat.m[6]*m[ 1] + mat.m[10]*m[ 2] + mat.m[14]*m[ 3],
						mat.m[3]*m[ 0] + mat.m[7]*m[ 1] + mat.m[11]*m[ 2] + mat.m[15]*m[ 3],
						mat.m[0]*m[ 4] + mat.m[4]*m[ 5] + mat.m[ 8]*m[ 6] + mat.m[12]*m[ 7],
						mat.m[1]*m[ 4] + mat.m[5]*m[ 5] + mat.m[ 9]*m[ 6] + mat.m[13]*m[ 7],
						mat.m[2]*m[ 4] + mat.m[6]*m[ 5] + mat.m[10]*m[ 6] + mat.m[14]*m[ 7],
						mat.m[3]*m[ 4] + mat.m[7]*m[ 5] + mat.m[11]*m[ 6] + mat.m[15]*m[ 7],
						mat.m[0]*m[ 8] + mat.m[4]*m[ 9] + mat.m[ 8]*m[10] + mat.m[12]*m[11],
						mat.m[1]*m[ 8] + mat.m[5]*m[ 9] + mat.m[ 9]*m[10] + mat.m[13]*m[11],
						mat.m[2]*m[ 8] + mat.m[6]*m[ 9] + mat.m[10]*m[10] + mat.m[14]*m[11],
						mat.m[3]*m[ 8] + mat.m[7]*m[ 9] + mat.m[11]*m[10] + mat.m[15]*m[11],
						mat.m[0]*m[12] + mat.m[4]*m[13] + mat.m[ 8]*m[14] + mat.m[12]*m[15],
						mat.m[1]*m[12] + mat.m[5]*m[13] + mat.m[ 9]*m[14] + mat.m[13]*m[15],
						mat.m[2]*m[12] + mat.m[6]*m[13] + mat.m[10]*m[14] + mat.m[14]*m[15],
						mat.m[3]*m[12] + mat.m[7]*m[13] + mat.m[11]*m[14] + mat.m[15]*m[15] );
	}
};


// --------------------------------------------------------
// 外部函数 - 矩阵运算 +, -, * ( Matrix4 vs Matrix4)
template<typename T>
inline const Matrix4<T> operator + (const Matrix4<T> &mat1,const Matrix4<T> &mat2)
{
	return Matrix4<T>( mat1.m[ 0] + mat2.m[ 0], mat1.m[ 1] + mat2.m[ 1], mat1.m[ 2] + mat2.m[ 2], mat1.m[ 3] + mat2.m[ 3],
					   mat1.m[ 4] + mat2.m[ 4], mat1.m[ 5] + mat2.m[ 5], mat1.m[ 6] + mat2.m[ 6], mat1.m[ 7] + mat2.m[ 7],
					   mat1.m[ 8] + mat2.m[ 8], mat1.m[ 9] + mat2.m[ 9], mat1.m[10] + mat2.m[10], mat1.m[11] + mat2.m[11],
					   mat1.m[12] + mat2.m[12], mat1.m[13] + mat2.m[13], mat1.m[14] + mat2.m[14], mat1.m[15] + mat2.m[15] );
}

template<typename T>
inline const Matrix4<T> operator - (const Matrix4<T> &mat1,const Matrix4<T> &mat2)
{
	return Matrix4<T>( mat1.m[ 0] - mat2.m[ 0], mat1.m[ 1] - mat2.m[ 1], mat1.m[ 2] - mat2.m[ 2], mat1.m[ 3] - mat2.m[ 3],
					   mat1.m[ 4] - mat2.m[ 4], mat1.m[ 5] - mat2.m[ 5], mat1.m[ 6] - mat2.m[ 6], mat1.m[ 7] - mat2.m[ 7],
					   mat1.m[ 8] - mat2.m[ 8], mat1.m[ 9] - mat2.m[ 9], mat1.m[10] - mat2.m[10], mat1.m[11] - mat2.m[11],
					   mat1.m[12] - mat2.m[12], mat1.m[13] - mat2.m[13], mat1.m[14] - mat2.m[14], mat1.m[15] - mat2.m[15]);
}

template<typename T>
inline const Matrix4<T> operator * (const Matrix4<T> &mat1,const Matrix4<T> &mat2)
{
	return Matrix4<T>( mat1.m[0]*mat2.m[ 0] + mat1.m[4]*mat2.m[ 1] + mat1.m[ 8]*mat2.m[ 2] + mat1.m[12]*mat2.m[ 3],
					    mat1.m[1]*mat2.m[ 0] + mat1.m[5]*mat2.m[ 1] + mat1.m[ 9]*mat2.m[ 2] + mat1.m[13]*mat2.m[ 3],
					    mat1.m[2]*mat2.m[ 0] + mat1.m[6]*mat2.m[ 1] + mat1.m[10]*mat2.m[ 2] + mat1.m[14]*mat2.m[ 3],
					    mat1.m[3]*mat2.m[ 0] + mat1.m[7]*mat2.m[ 1] + mat1.m[11]*mat2.m[ 2] + mat1.m[15]*mat2.m[ 3],
					    mat1.m[0]*mat2.m[ 4] + mat1.m[4]*mat2.m[ 5] + mat1.m[ 8]*mat2.m[ 6] + mat1.m[12]*mat2.m[ 7],
					    mat1.m[1]*mat2.m[ 4] + mat1.m[5]*mat2.m[ 5] + mat1.m[ 9]*mat2.m[ 6] + mat1.m[13]*mat2.m[ 7],
					    mat1.m[2]*mat2.m[ 4] + mat1.m[6]*mat2.m[ 5] + mat1.m[10]*mat2.m[ 6] + mat1.m[14]*mat2.m[ 7],
					    mat1.m[3]*mat2.m[ 4] + mat1.m[7]*mat2.m[ 5] + mat1.m[11]*mat2.m[ 6] + mat1.m[15]*mat2.m[ 7],
					    mat1.m[0]*mat2.m[ 8] + mat1.m[4]*mat2.m[ 9] + mat1.m[ 8]*mat2.m[10] + mat1.m[12]*mat2.m[11],
					    mat1.m[1]*mat2.m[ 8] + mat1.m[5]*mat2.m[ 9] + mat1.m[ 9]*mat2.m[10] + mat1.m[13]*mat2.m[11],
					    mat1.m[2]*mat2.m[ 8] + mat1.m[6]*mat2.m[ 9] + mat1.m[10]*mat2.m[10] + mat1.m[14]*mat2.m[11],
					    mat1.m[3]*mat2.m[ 8] + mat1.m[7]*mat2.m[ 9] + mat1.m[11]*mat2.m[10] + mat1.m[15]*mat2.m[11],
					    mat1.m[0]*mat2.m[12] + mat1.m[4]*mat2.m[13] + mat1.m[ 8]*mat2.m[14] + mat1.m[12]*mat2.m[15],
					    mat1.m[1]*mat2.m[12] + mat1.m[5]*mat2.m[13] + mat1.m[ 9]*mat2.m[14] + mat1.m[13]*mat2.m[15],
					    mat1.m[2]*mat2.m[12] + mat1.m[6]*mat2.m[13] + mat1.m[10]*mat2.m[14] + mat1.m[14]*mat2.m[15],
					    mat1.m[3]*mat2.m[12] + mat1.m[7]*mat2.m[13] + mat1.m[11]*mat2.m[14] + mat1.m[15]*mat2.m[15] );
}

// --------------------------------------------------------
// 外部函数 - 矩阵运算( Matrix4 vs Others)
template<typename T>
inline const Tuple3<T> operator * (const Tuple3<T> &t,const Matrix4<T> &mat)
{
	return Tuple3<T>( mat.m[ 0]*t.x + mat.m[ 4]*t.y + mat.m[ 8]*t.z + mat.m[12],
					  mat.m[ 1]*t.x + mat.m[ 5]*t.y + mat.m[ 9]*t.z + mat.m[13],
					  mat.m[ 2]*t.x + mat.m[ 6]*t.y + mat.m[10]*t.z + mat.m[14] );
}

template<typename T>
inline const Tuple3<T> operator * (const Matrix4<T> &mat,const Tuple3<T> &t)
{
	return Tuple3<T>( mat.m[ 0]*t.x + mat.m[ 4]*t.y + mat.m[ 8]*t.z + mat.m[12],
					  mat.m[ 1]*t.x + mat.m[ 5]*t.y + mat.m[ 9]*t.z + mat.m[13],
					  mat.m[ 2]*t.x + mat.m[ 6]*t.y + mat.m[10]*t.z + mat.m[14] );
}

template<typename T>
inline const Tuple4<T> operator * (const Matrix4<T> &mat, const Tuple4<T> &t)
{
	return Tuple4<T>( mat.m[ 0]*t.x + mat.m[ 4]*t.y + mat.m[ 8]*t.z + mat.m[12]*t.w,
					   mat.m[ 1]*t.x + mat.m[ 5]*t.y + mat.m[ 9]*t.z + mat.m[13]*t.w,
					   mat.m[ 2]*t.x + mat.m[ 6]*t.y + mat.m[10]*t.z + mat.m[14]*t.w,
					   mat.m[ 3]*t.x + mat.m[ 7]*t.y + mat.m[11]*t.z + mat.m[15]*t.w  );
}

template<typename T>
inline const Tuple4<T> operator * (const Tuple4<T> &t,const Matrix4<T> &mat)
{
	return Tuple4<T>( mat.m[ 0]*t.x + mat.m[ 4]*t.y + mat.m[ 8]*t.z + mat.m[12]*t.w,
					   mat.m[ 1]*t.x + mat.m[ 5]*t.y + mat.m[ 9]*t.z + mat.m[13]*t.w,
					   mat.m[ 2]*t.x + mat.m[ 6]*t.y + mat.m[10]*t.z + mat.m[14]*t.w,
					   mat.m[ 3]*t.x + mat.m[ 7]*t.y + mat.m[11]*t.z + mat.m[15]*t.w  );
}

// -------------------------------------
typedef Matrix4<int>    Matrix4i;
typedef Matrix4<float>  Matrix4f;
typedef Matrix4<double> Matrix4d;

#endif


• 求解线性方程组

        线性方程组 是我们经常会遇到的一个问题，常用的 解线性方程组常用的方法有两个：

1. 高斯消元法；

2. 利用行列式，即Crammer法则 求解。

        这里我们只讲 克莱姆法则，该方法的前提条件是：

1. 线性方程组 AX=b 方程的个数与未知量的个数相同，即系数矩阵A是一个方阵；

2. 系数矩阵A的行列式 |A| ≠ 0。则方程组有唯一解： xi = Di/D

    D=|A|

    Di 是 D 中第 i 列换成 b 得到的行列式。

        下面给出Crammer法则 的代码：

// 利用克莱姆法则求方程解
// a11 a12 a13   x1   b1
// a21 a22 a23 * x2 = b2
// a31 a32 a33   x3   b3
// 或者写成矩阵形式为Ax=b，其中A为n*n方阵，x为n个变量构成列向量，b为n个常数项构成列向量。
bool CramerSolver(double a11,double a12,double a13,double a21,double a22,double a23,double a31,double a32,double a33,double b1,double b2,double b3,double& x1,double& x2,double& x3)
{
	// 首先计算行列式的值
	double D = Determinant( a11, a12, a13, 
				a21, a22, a23, 
				a31, a32, a33 );
	if(0==D)
	    return false;

	// 当它的系数矩阵可逆，或者说对应的行列式|A|不等于0的时候，它有唯一解xi=|Ai|/|A|，其中Ai〔i = 1，2，……，n〕是矩阵A中第i列的a 1i，a 2i，……a ni （即第i列）依次换成b1，b2，……bn所得的矩阵。
	//D1=|b1 a12 a13| 
	//	 |b2 a22 a23| 
	//	 |b3 a32 a33|
	double D1 = Determinant( b1, a12, a13, 
				b2, a22, a23, 
				b3, a32, a33 );
	double D2 = Determinant( a11, b1, a13, 
				a21, b2, a23, 
				a31, b3, a33 );
	double D3 = Determinant( a11, a12, b1, 
				a21, a22, b2, 
	                    	a31, a32, b3 );

	x1 = D1/D;
	x2 = D2/D;
	x3 = D3/D;

	return true;
}
linolzhang
关注
专栏目录
矩阵的运算
也许会有_hui的博客
06-14 1033
目录 1.矩阵与数相乘 每一项都要乘 2.矩阵的加减运算 每一项都要乘​ 3.矩阵相乘 4.矩阵对应元素相乘 ---同型矩阵 5.矩阵的转置 T 6.矩阵的共轭转置 H 7.矩阵的逆 I 8.矩阵的试图 A 1.矩阵与数相乘 每一项都要乘 import numpy as np m1 = np.mat([[1, 2, 3], [2, 3, 4]]) pr...
c语言稀疏矩阵运算器,数据结构：稀疏矩阵运算器
weixin_35574537的博客
05-19 3179
题目：设计一个程序，实现一个能进行稀疏矩阵基本运算的计算器。按照教科书《数据结构(C语言版)》(严蔚敏等)5.3.2节中描述的方法，以十字链表表示稀疏矩阵。一、需求分析稀疏矩阵是指那些多数元素为零的矩阵。利用“稀疏”特点进修学校储存和计算可以大大节省储存空间，提高计算效率。实现一个能进行稀疏矩阵基本运算的计算器。以“带行逻辑链接信息”的三元组顺序表表示稀疏矩阵，实现两个矩阵相加、相减和相乘的运算。...
矩阵运算
aoniebu42564的博客
10-06 236
矩阵加、减法 矩阵加法非常简单，对应位置直接加减即可，但是前提是两个矩阵大小相同（即一个矩阵是N*M的，另一个与之相加的矩阵的大小也要是N*M）。就像这样： 矩阵乘法 矩阵乘法就相对比较复杂了。他需要满足的前提是第一个矩阵的列数要等于第二个矩阵的行数，这样的两个矩阵才可以相乘。下面我用一个图来解释怎么样进行矩阵乘法： 矩阵乘法性质： ①矩阵乘法满足乘...
原创CODESYS操作Matrix33矩阵运算的功能块库文件
05-22
内容概要：原创的CODESYS操作Matrix3阶方阵矩阵运算的功能块的编译库。调用库内功能块可便捷实现对Matrix3阶方矩矩阵运算。 1，3阶方矩加法，减法，点乘，叉乘，左乘常数，右乘1列。 2，3阶方矩的求行列式，转置，...
用Matlab学习线性代数-线性方程组与矩阵代数概要.pdf
11-15
在学习线性代数时，矩阵代数和线性方程组是核心概念。Matlab 是一个强大的工具，可以帮助我们理解和解决这些问题。...通过实际操作，我们不仅掌握了基本的矩阵运算，还学会了如何在 Matlab 中高效地解决线性代数问题。
稀疏矩阵三元组表的矩阵运算与Python实现
最新发布
10-18
内容概要：本文详细介绍了使用三元组表表示稀疏矩阵及其常见的运算方法，包括矩阵加法、乘法和转置运算。文中提供了基于Python语言的具体实现代码，涵盖了三元组类定义、矩阵创建以及各类运算的方法。 适合人群：...
数据结构课程设计特殊矩阵运算知识.pdf
02-05
数据结构课程设计中，特殊矩阵运算是一个重要的实践课题，主要目标是实现矩阵的基本运算，如加法、减法、转置以及数乘等。在处理这类问题时，通常会运用到两种数据结构：二维数组和类似图的数据结构，以适应不同类型...
矩陣的運算
I am Rocky
02-14 3562
http://blog.csdn.net/as3_flash/article/details/2378453 把一組數字記錄成矩陣形式是沒有意義的，學習矩陣的關鍵在於掌握矩陣之間的運算。 1．矩陣加法運算 矩陣之間也可以相加。把兩個矩陣對應位置的單個元素相加，得到的新矩陣就是矩陣加法的結果。由其運算法則可知，只有行數和列數完全相同的矩陣才能進行加法運算。 矩陣之間相加沒有順序，假設A、B都
线性代数--矩阵基本计算（加减乘法）
qq_51216031的博客
11-27 6万+
对线性代数的矩阵的运算有详细的总结，来源于B站宋浩的视频所作的学习笔记。
数学基础 - 矩阵的基本运算（Matrix Operations）
DarkRabbit的专栏
04-21 4万+
矩阵的基本运算(Matrix Operations) 三个初等行（列）变换 加法(Plus) 乘法(Multiply) 与数的乘法 与矩阵的乘法 转置(Transpose) 方阵的行列式(Determinant) 克莱默法则(Cramer) 雅可比行列式(Jacobi) 逆(Inverse) 秩(Rank) 迹(Trace) 矩阵的基本运算(Matrix Operati...
矩阵的基本运算
热门推荐
palawind
05-26 10万+
一、矩阵的加法与减法 　　1、运算规则 　　设矩阵，， 　　则 　　　　　 　　简言之，两个矩阵相加减，即它们相同位置的元素相加减！ 　　注意：只有对于两个行数、列数分别相等的矩阵（即同型矩阵），加减法运算才有意义，即加减运算是可行的．　　2、 运算性质 （假设运算都是可行的） 　　满足交换律和结合律　　交换律　 ； 　　结合律　 ．二、矩阵与数的乘法　　1、 运算规则 　　数乘矩阵A，就是...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值