cmatrix.cpp

GPS坐标转换软件与源程序 支持世界上大多数坐标框架下

i4 CMatrix::CalcuNorm()
{
	i4 nElementNum = m_nRowNumber * m_nColumnNumber;
	m_dNorm = 0;
	for(i4 i=0;i<nElementNum;i++)
	{
		m_dNorm += m_pData[i] * m_pData[i];
	}
	m_dNorm = sqrt(m_dNorm);
	return 1;
}


//功能：转换地址中的数值,矩阵转置核心函数
/*参数：m1:原数据首地址
		m2：转置后数据的首地址
		m：原数组行数
		n：原数组列数*/
void CMatrix::Transpose(const f8 *m1,f8 *m2,i4 m,i4 n)const
{
	i4 i,j;
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			m2[j*m+i]=m1[i*n+j];
	return;
}



//功能：两个矩阵相乘核心函数
/*参数：m1:		第一矩阵数据首地址
		m2：	第二矩阵数据首地址
		result：结果矩阵数据首地址
		i_1：	第一矩阵行数，即结果矩阵行数
		j_12：	第一矩阵列数，即第二矩阵的行数
		j_2：	第二矩阵的列数，即结果矩阵列数*/
void CMatrix::Mult(const f8 *m1,const f8 *m2, f8 *result,i4 i_1,i4 j_12,i4 j_2)const
{
	i4 i,j,k;
	for(i=0;i<i_1;i++)
		for(j=0;j<j_2;j++){
			result[i*j_2+j]=0.0;
			for(k=0;k<j_12;k++)
				result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];
							}
	return;
}



//功能：方矩阵求逆核心函数
/*参数：m1:		原矩阵数据首地址，也是逆阵的首地址
		n:		矩阵行数或列数
//返回值：成功返回1，失败返回0*/
i4 CMatrix::InversMatrix(f8 *m1,i4 n)const
{ 
	i4 *is,*js;
	i4 i,j,k,l,u,v;
	f8 temp,max_v;
	is = new i4[n];
	js = new i4[n];
	if(is==NULL||js==NULL)	return(0);

	for(k=0;k<n;k++)
	{
		max_v = 0.0;
		for(i=k;i<n;i++)
			for(j=k;j<n;j++)
			{
				temp = fabs(m1[i*n+j]);
				if( temp>max_v )
				{
			        max_v = temp; 
					is[k] = i; 
					js[k] = j;
				}
			}
		if(max_v==0.0)
		{
			delete []is; 
			delete []js;
			return(0);//printf("invers is not availble!\n");
		}
		if(is[k]!=k)
			for(j=0;j<n;j++)
			{
			   u = k*n+j;
			   v = is[k]*n+j;
			   temp = m1[u]; 
			   m1[u] = m1[v];
			   m1[v] = temp;
			}
		if(js[k]!=k)
			for(i=0;i<n;i++){
				u=i*n+k; v=i*n+js[k];
				temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
			}
		l=k*n+k;
		m1[l]=1.0/m1[l];
		for(j=0;j<n;j++)
			if(j!=k){
				u=k*n+j;
				m1[u]*=m1[l];
			}
		for(i=0;i<n;i++)
			if(i!=k)
				for(j=0;j<n;j++)
					if(j!=k){
					  u=i*n+j;
					  m1[u]-=m1[i*n+k]*m1[k*n+j];
					}
		for(i=0;i<n;i++)
			if(i!=k){
				u=i*n+k;
				m1[u]*=-m1[l];
			}
	}
	for(k=n-1;k>=0;k--){
		if(js[k]!=k)
			for(j=0;j<n;j++){
				u=k*n+j; v=js[k]*n+j;
				temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
			}
		if(is[k]!=k)
			for(i=0;i<n;i++){
				u=i*n+k; v=i*n+is[k];
				temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
			}
	}
	delete []is; delete []js;

  return 1;

}


//功能：将矩阵元素设零
//参数：无
//返回值：1
i4 CMatrix::SetZero()
{
	if(m_pData != NULL)
	{
		memset(m_pData,0,sizeof(f8)*m_nRowNumber * m_nColumnNumber);
	}
	return 1;
}


//功能：将矩阵元素设负
//参数：无
//返回值：1
i4 CMatrix::SetMinus()
{
	i4 nElementNumber = m_nRowNumber * m_nColumnNumber;
	for(i4 i=0; i<nElementNumber; i++)
	{
		*(m_pData + i) *= -1;
	}
	
	return 1;
}


//功能：将矩阵元素乘以倍数----注意：此函数和“*”操作符重载不同，这个函数视将自身元素改变，而后者则自身不变 
//参数：f8MulNum	倍数
//返回值：1
i4 CMatrix::Multiply(f8 f8MulNum)
{
	i4 nElementNumber = m_nRowNumber * m_nColumnNumber;
	for(i4 i=0;i<nElementNumber;i++)
	{
		m_pData[i] *= f8MulNum;
	}
	
	return 1;
}


//功能：将矩阵设成单位阵
//参数：nRowNumber	单位阵的维数
//返回值：1
i4 CMatrix::SetIdentity(i4 nRowNumber)
{
	if(nRowNumber==0) return 0;

	this->InitEmptyMatrix(nRowNumber,nRowNumber);//开空间
	memset(m_pData,0,sizeof(f8)*nRowNumber * nRowNumber);//第一步全部置零
	for(i4 j=0;j<nRowNumber;j++)
	{
		*(m_pData + j*(nRowNumber+1) ) = 1.0;//第二步将对角线置1
	}

	return 1;
}


//功能：将矩阵设成简单相关阵
//参数：nRowNumber	相关阵的维数
//返回值：1
i4 CMatrix::SetCorrelativity(i4 nRowNumber)
{
	if(nRowNumber==0) return 0;

	if(nRowNumber * nRowNumber != m_nRowNumber * m_nColumnNumber)//空间不相同
	{
		this->InitEmptyMatrix(nRowNumber, nRowNumber);//重新开空间
	}
	
	for(i4 i=0;i<nRowNumber*nRowNumber;i++)
	{
		*(m_pData + i ) = 1.0;//第一步全部置1
	}
	for(i4 j=0;j<nRowNumber;j++)
	{
		*(m_pData + j*(nRowNumber+1) ) = 2.0;//第二步将对角线置2
	}

	return 1;
}




///////////////////////////////////////////////////
//公用操作符重载    operator override
///////////////////////////////////////////////////

//功能：赋值操作符重载
/*参数：SecMatrix:		源矩阵
//返回值：自身矩阵的引用，这样可以使用连等形式，如: A = B = C ，此时 A , B 都等于 C 的值;
//实现方法：调用InitMatrix()函数实现功能
使用示例：A=B=C 意为: B=C; A=B; */
CMatrix& CMatrix::operator= (const CMatrix & SecMatrix)
{
	if(&SecMatrix == this) return *this;
	this->InitMatrix(SecMatrix.GetDataAddress(),SecMatrix.GetRowNumber(),SecMatrix.GetColumnNumber());

	return *this;
}




//功能：矩阵加法操作符重载
/*参数：SecMatrix:		源矩阵
//返回值：相加后的结果矩阵
使用示例：A=B+C 意为: A中的元素为B、C中元素之和 */
CMatrix CMatrix::operator+ (const CMatrix & SecMatrix)const
{
	CMatrix TmpMatrix;

	//行、列不相等的矩阵不能相加
	if((this->m_nRowNumber != SecMatrix.GetRowNumber()) || (this->m_nColumnNumber != SecMatrix.GetColumnNumber())) return TmpMatrix;

	i4 nElementNumber = m_nRowNumber * m_nColumnNumber;
	f8* pTmpData2 = new f8[nElementNumber];
	SecMatrix.GetData(pTmpData2);
	for(i4 i=0; i<nElementNumber; i++)
	{
		*(pTmpData2 + i) += *(this->m_pData + i);
	}
	
	TmpMatrix.InitMatrix(pTmpData2,m_nRowNumber,m_nColumnNumber);
	delete [] pTmpData2;

	return TmpMatrix;
}




//功能：矩阵减法操作符重载
/*参数：SecMatrix:		源矩阵
//返回值：相加后的结果矩阵
使用示例：A=B-C 意为: A中的元素为B、C中元素之差*/
CMatrix CMatrix::operator- (const CMatrix & SecMatrix)const
{
	CMatrix TmpMatrix( SecMatrix );
	TmpMatrix.SetMinus();

	CMatrix SecTmpMatrix(* this + TmpMatrix);
	return SecTmpMatrix;
}




//功能：取负操作符重载
/*参数：无
//返回值：自身矩阵的负矩阵
//实现方法：调用SetMinus()函数实现功能
使用示例：A=-C 意为: A中的元素为B中元素的负值 */
CMatrix CMatrix::operator- ()
{
	CMatrix TmpMatrix( * this );
	TmpMatrix.SetMinus();
	return TmpMatrix;
}




//功能：矩阵相乘
/*参数：SecMatrix:		被乘矩阵
//返回值：矩阵相乘的结果矩阵
使用示例：C=A*B C为结果矩阵 */
CMatrix CMatrix::operator* (const CMatrix & SecMatrix)const
{
	i4 nElementNumber3 = m_nRowNumber * SecMatrix.GetColumnNumber();
	f8* pTmpData3 = new f8[nElementNumber3];
	this->Mult(this->m_pData,SecMatrix.GetDataAddress(),pTmpData3,m_nRowNumber,m_nColumnNumber,SecMatrix.GetColumnNumber());

	CMatrix TmpMatrix(pTmpData3,m_nRowNumber,SecMatrix.GetColumnNumber());
	delete [] pTmpData3;

	return TmpMatrix;
}




//功能：矩阵乘以实数
/*参数： dMultiNum:		被乘数
//返回值：矩阵乘以实数的结果矩阵
使用示例：C= A * dMultiNum C为结果矩阵 */
CMatrix CMatrix::operator *(const f8 dMultiNum) const
{
	i4 nElementNumber = m_nRowNumber * m_nColumnNumber;
	f8* pTmpData = new f8[nElementNumber];
	this->GetData(pTmpData);
	for(i4 i=0;i<nElementNumber;i++)
	{
		pTmpData[i] *= dMultiNum;
	}
	
	CMatrix TmpMatrix(pTmpData,m_nRowNumber,m_nColumnNumber);
	delete [] pTmpData;

	return TmpMatrix;

}




//功能：矩阵除以实数
/*参数： dMultiNum:		被除数，不能为0
//返回值：矩阵乘以实数的结果矩阵
使用示例：C= A * dMultiNum C为结果矩阵 */
CMatrix CMatrix::operator /(const f8 dDevideNum) const
{
	CMatrix TmpMatrix;
	if(dDevideNum == 0) return TmpMatrix;
	
	TmpMatrix = this->operator *(1/dDevideNum);

	return TmpMatrix;
}





//功能：矩阵正向除
/*参数：SecMatrix:		被除矩阵
//返回值：矩阵正向除的结果矩阵 */
//使用方法：C=A/B C为A乘以B的逆阵的结果阵
CMatrix CMatrix::operator/ (const CMatrix & SecMatrix)const
{
	//先求B的逆阵
	CMatrix TmpMatrix1; SecMatrix.GetInverseMatrix(TmpMatrix1);
	if(!TmpMatrix1.GetRowNumber()) return TmpMatrix1;//如果求逆失败，则返回空阵

	//再矩阵相乘
	CMatrix TmpMatrix2 = (*this) * TmpMatrix1;

	return TmpMatrix2;
}



//功能：矩阵逆向除，A%B:表示A的逆阵乘以B
/*参数：SecMatrix:		被操作矩阵
//返回值：矩阵逆向除的结果矩阵 */
//使用方法：C=A%B C为A的逆阵乘以B的结果阵
CMatrix CMatrix::operator% (const CMatrix & SecMatrix)const
{
	//先求A的逆阵
	CMatrix TmpMatrix1; this->GetInverseMatrix(TmpMatrix1);
	if(!TmpMatrix1.GetRowNumber()) return TmpMatrix1;//如果求逆失败，则返回空阵

	//再矩阵相乘
	CMatrix TmpMatrix2 = TmpMatrix1 * SecMatrix;

	return TmpMatrix2;
}

i4	CMatrix::PrintFile(FILE& fp)
{
	for (i4 row=0;row<m_nRowNumber;row++)
	{
		for (i4 column=0;column<m_nColumnNumber;column++)
		{
			if(row >= m_nRowNumber || column>= m_nColumnNumber) return 0;
			f8 data=*(m_pData+ row * m_nColumnNumber + column);
			if ((fabs(data)>0.001 && data <1e5 && data >-1e5) || fabs(data) <1e-10)
				fprintf(&fp,"%9.3f,",data);
			else if(data<-1e5)
			{
				//modify the index digit to 2;
				int     decimal,   sign;	//decimal point position and sign.
				char    *buffer;			//store string
				char	csign;
				int     precision = 3;		//Number of digits stored
				buffer = _ecvt( data, precision, &decimal, &sign );
				sign == 0 ?csign=' ':csign='-';
				fprintf(&fp,"%c%c.%se%d,",csign,buffer[0],buffer+1,decimal-1);

			}
 			else if(data>1e5)
 			{
				//modify the index digit to 2;
				int     decimal,   sign;	//decimal point position and sign.
				char    *buffer;			//store string
				char	csign;
				int     precision = 4;		//Number of digits stored
				buffer = _ecvt( data, precision, &decimal, &sign );
				sign == 0 ?csign=' ':csign='-';
				fprintf(&fp," %c%c.%se%d,",csign,buffer[0],buffer+1,decimal-1);
 
 			}
			else
			{
				//modify the index digit to 2;
				int     decimal,   sign;	//decimal point position and sign.
				char    *buffer;			//store string
				char	csign;
				int     precision = 4;		//Number of digits stored
				buffer = _ecvt( data, precision, &decimal, &sign );
				sign == 0 ?csign=' ':csign='-';
				decimal<0 ? fprintf(&fp,"%c%c.%se%d,",csign,buffer[0],buffer+1,decimal-1):
							fprintf(&fp," %c%c.%se%d,",csign,buffer[0],buffer+1,decimal-1);
			}

		}
		fprintf(&fp,"\n");
	}
	return 1;
}