Matrix

矩阵乘法

bool matrix_product(const double* A,const double * B,double* C,
	int ARowNum,int AColNum,int BColNum)
{
	// memory error
	if(A == NULL||B == NULL || C == NULL)
		return false;
	// column or row number error
	if(!(ARowNum * AColNum *BColNum) )
		return false;
	// calculate matrix product	
	for(int rowNum = 0;rowNum < ARowNum;rowNum++)
	{
		for(int colNum = 0;colNum < BColNum; colNum++)
		{
			double sum = 0;
			for(int index = 0;index < AColNum;index++)
			{
				sum+=A[rowNum*AColNum+index]*B[index*BColNum+colNum];
			}
			C[rowNum*BColNum+colNum]= sum;
		}
	}
	return true;
}

矩阵转秩

bool matrix_transpose(const double* A,double* B,int ARowNum,int AColNum)
{
	// memory error
	if(A == NULL||B == NULL || C == NULL)
		return false;
	// column or row number error
	if(!(ARowNum * AColNum *BColNum) )
		return false;
	// swap data a[i][j] = a[j][i]
	for(int rowNum = 0;rowNum < AColNum;rowNum++)
	{
		for(int colNum = 0;colNum < ARowNum; colNum++)
		{
			B[rowNum*ARowNum+colNum] = A[colNum*AColNum+rowNum];
		}
	}

	return true;
}

矩阵求逆

bool matrix_inverse(const double* A,double* A_inv,int rowNum)
{
	// memory error
	if(A == NULL||B == NULL || C == NULL)
		return false;
	// column or row number error
	if(!(ARowNum * AColNum *BColNum) )
		return false;
	// calculate determinant of A	
	double d = calculateDeterminant(rowNum,A);
	// A is degradation.
	if(d == 0)   
		return false;
		
	for(int i= 0;i<rowNum;i++)
	{
		for(int j=0;j<rowNum;j++)
		{
			A_inv[i*rowNum+j] = CalculateAlgebraicComplement(i,j,A,rowNum)/d;
		}
	}
	return true;
}

double calculateDeterminant(int rowNum,const double* A)
{
	double* B = new double[rowNum*rowNum];
	memcpy(B,A,rowNum*rowNum*sizeof(double));
	int sign = 0;
	double result = 1.0;
	for(int i= 0;i<rowNum;i++)
	{
		if(B[i*rowNum+i] == 0)
		{
			if(!AdjustDeterminant(B,i,rowNum))
			{
				delete B;
				return 0;
			}
				sign++;
		}

		Elimination(B,i,rowNum);	
		result*= B[i*rowNum+i];
	}	
	delete B;
	return result*(sign%2?-1:1);
}

bool AdjustDeterminant(double* A,int i,int rowNum)
{
	for(int k= i+1;k<rowNum;k++)
	{
		if(A[k*rowNum+i] != 0)
		{
			SwapRow(A,i,k,rowNum);
			break;
		}
	}	

	if(A[i*rowNum+i] == 0)
		return false;
	return true;
}

void Elimination(double*A,int i,int rowNum)
{
	for(int j= i+1;j<rowNum;j++)
	{
		double scale = -1*(A[j*rowNum+i]/A[i*rowNum+i]);
		for(int k= i;k<rowNum;k++)
		{
			A[j*rowNum+k] += A[i*rowNum+k]*scale;
		}
	}
}

void SwapRow(double* A,int i,int j,int rowNum )
{
	double mid;
	for(int k=0;k<rowNum;k++)
	{
		mid = A[i*rowNum+k];
		A[i*rowNum+k] = A[j*rowNum+k];
		A[j*rowNum+k] = mid;
	}	
 }

double CalculateAlgebraicComplement(int rowIndex,int colIndex,const double*A,int rowNum)
{
	double* adjointMatrix = GetAdjointMatrix(rowIndex,colIndex,A,rowNum);	
	double result =  calculateDeterminant(rowNum-1,adjointMatrix);
	delete adjointMatrix;
	result *= (rowIndex+colIndex)%2?-1:1;	
	return result;
}

double* GetAdjointMatrix(int rowIndex,int colIndex,const double*A,int rowNum)
{
	double* adjointMatrix = new double[(rowNum-1)*(rowNum-1)];
	int count = 0;
	for(int i= 0;i<rowNum;i++)
	{
		if(i==rowIndex)
			continue;
		for(int j=0;j<rowNum;j++)
		{
			if(j==colIndex)
				continue;
			adjointMatrix[count] = A[i*rowNum+j];
			count++;
		}
	}
	return adjointMatrix;
}

矩阵数乘

bool matrix_scale(double *A,double scale,int rowNum,int colNum)
{
	if(A == NULL)
		return false;
	if(!(rowNum*colNum))
		return false;
	for(int i = 0;i<rowNum;i++)
		for(int j=0;j<colNum;j++)
			A[i*colNum+j]  *=scale;
	return true;		
}

矩阵加法

bool matrix_Sum(const double* A,const double*B,double* C,int ARowNum,int AColNum,int BRowNum,int BColNum)
{
	if(A == NULL||B == NULL || C == NULL)
		return false;	
	if(ARowNum != BRowNum||AColNum != BColNum)
		return false;
	if(!(ARowNum*AColNum))
		return false;
	for(int i= 0;i<ARowNum;i++)
		for(int j=0;j<AColNum;j++)
			C[i*AColNum+j] = A[i*AColNum+j]+B[i*AColNum+j];


	return true;
}