blas3.cpp

图像分割算法

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



#include "./Blas.h"
#include "mkl.h"




namespace Blas
{


	/// Gemm: General m*m
	/// C = A*B
	/// A is an m by k matrix, B is a k by n matrix, and C is an m by n matrix.
	Matrix<float> Gemm(Matrix<float>& A, Matrix<float>& B)
	{
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		Matrix<float> C(A.Rows(), B.Columns(), 0);
		cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
		return C;
	}

	Matrix<double> Gemm(Matrix<double>& A, Matrix<double>& B)
	{
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		Matrix<double> C(A.Rows(), B.Columns(), 0);
		cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexFloat> Gemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
	{
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		ComplexFloat temp(1,0);
		Matrix<ComplexFloat> C(A.Rows(), B.Columns(), ComplexFloat(0,0));
		cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexDouble> Gemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
	{
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		ComplexDouble temp(1,0);
		Matrix<ComplexDouble> C(A.Rows(), B.Columns(), ComplexDouble(0,0));
		cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
		return C;
	}



	/// Hermitian or symmetric A
	///  C = A*B or C = B*A (if orderReversed == true)
	Matrix<float> Symm(Matrix<float>& A, Matrix<float>& B)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		Matrix<float> C(B.Rows(), B.Columns(), 0);
		cblas_ssymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
		return C;
	}

	Matrix<double> Symm(Matrix<double>& A, Matrix<double>& B)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		Matrix<double> C(B.Rows(), B.Columns(), 0);
		cblas_dsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexFloat> Symm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		ComplexFloat temp(1,0);
		Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
		cblas_csymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexDouble> Symm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		ComplexDouble temp(1,0);
		Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
		cblas_zsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexFloat> Hemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		ComplexFloat temp(1,0);
		Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
		cblas_chemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexDouble> Hemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(A.Columns() != B.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
		}
		ComplexDouble temp(1,0);
		Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
		cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
		return C;
	}


	Matrix<float> Symm(Matrix<float>& A, Matrix<float>& B, bool orderReversed)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(orderReversed)
		{
			if(A.Rows() != B.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			Matrix<float> C(B.Rows(), B.Columns(), 0);
			cblas_ssymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
			return C;
		}
		else
		{
			if(A.Columns() != B.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			Matrix<float> C(B.Rows(), B.Columns(), 0);
			cblas_ssymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
			return C;
		}
	}

	Matrix<double> Symm(Matrix<double>& A, Matrix<double>& B, bool orderReversed)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(orderReversed)
		{
			if(A.Rows() != B.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			Matrix<double> C(B.Rows(), B.Columns(), 0);
			cblas_dsymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
			return C;
		}
		else
		{
			if(A.Columns() != B.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			Matrix<double> C(B.Rows(), B.Columns(), 0);
			cblas_dsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
			return C;
		}
	}

	Matrix<ComplexFloat> Symm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, bool orderReversed)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(orderReversed)
		{
			if(A.Rows() != B.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexFloat temp(1,0);
			Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
			cblas_csymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
		else
		{
			if(A.Columns() != B.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexFloat temp(1,0);
			Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
			cblas_csymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
	}

	Matrix<ComplexDouble> Symm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, bool orderReversed)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(orderReversed)
		{
			if(A.Rows() != B.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexDouble temp(1,0);
			Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
			cblas_zsymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
		else
		{
			if(A.Columns() != B.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexDouble temp(1,0);
			Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
			cblas_zsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
	}

	Matrix<ComplexFloat> Hemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, bool orderReversed)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(orderReversed)
		{
			if(A.Rows() != B.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexFloat temp(1,0);
			Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
			cblas_chemm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
		else
		{
			if(A.Columns() != B.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexFloat temp(1,0);
			Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
			cblas_chemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
	}

	Matrix<ComplexDouble> Hemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, bool orderReversed)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(orderReversed)
		{
			if(A.Rows() != B.Columns())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;