blas1.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
{


	// ?ASUM

	float Asum(Vector<float>& x)
	{
		return cblas_sasum(x.Length(), x.DataPtr(), 1);
	}

	double Asum(Vector<double>& x)
	{
		return cblas_dasum(x.Length(), x.DataPtr(), 1);
	}

	float Asum(Vector<ComplexFloat>& x)
	{
		return cblas_scasum(x.Length(), x.DataPtr(), 1);
	}

	double Asum(Vector<ComplexDouble>& x)
	{
		return cblas_dzasum(x.Length(), x.DataPtr(), 1);
	}



	// ?AXPY

	void Axpy(float alpha, Vector<float>& x, Vector<float>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_saxpy(x.Length(), alpha, x.DataPtr(), 1, y.Data(), 1);
	}

	void Axpy(double alpha, Vector<double>& x, Vector<double>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_daxpy(x.Length(), alpha, x.DataPtr(), 1, y.Data(), 1);
	}

	void Axpy(ComplexFloat alpha, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_caxpy(x.Length(), &alpha, x.DataPtr(), 1, y.Data(), 1);
	}

	void Axpy(ComplexDouble alpha, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_zaxpy(x.Length(), &alpha, x.DataPtr(), 1, y.Data(), 1);
	}





	// ?COPY

	void Copy(Vector<float>& x, Vector<float>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_scopy(x.Length(), x.DataPtr(), 1, y.Data(), 1);
	}

	void Copy(Vector<double>& x, Vector<double>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_dcopy(x.Length(), x.DataPtr(), 1, y.Data(), 1);
	}

	void Copy(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_ccopy(x.Length(), x.DataPtr(), 1, y.Data(), 1);
	}

	void Copy(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_zcopy(x.Length(), x.DataPtr(), 1, y.Data(), 1);
	}


	// ?DOT


	float Dot(Vector<float>& x, Vector<float>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		return cblas_sdot(x.Length(), x.DataPtr(), 1, y.Data(), 1);
	}

	double Dot(Vector<double>& x, Vector<double>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		return cblas_ddot(x.Length(), x.DataPtr(), 1, y.Data(), 1);
	}

	ComplexFloat Dot(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y, bool conjugated)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		
		ComplexFloat result;
		if(conjugated)
		{
			cblas_cdotc_sub(x.Length(), x.DataPtr(), 1, y.Data(), 1, &result);
		}
		else
		{
			cblas_cdotu_sub(x.Length(), x.DataPtr(), 1, y.Data(), 1, &result);
		}
		return result;
	}

	ComplexDouble Dot(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, bool conjugated)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		
		ComplexDouble result;
		if(conjugated)
		{
			cblas_zdotc_sub(x.Length(), x.DataPtr(), 1, y.Data(), 1, &result);
		}
		else
		{
			cblas_zdotu_sub(x.Length(), x.DataPtr(), 1, y.Data(), 1, &result);
		}
		return result;
	}

	ComplexFloat Dot(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
	{
		return Dot(x, y, false);
	}

	ComplexDouble Dot(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
	{
		return Dot(x, y, false);
	}



	// ?NRM2

	float Nrm2(Vector<float>& x)
	{
		return cblas_snrm2(x.Length(), x.DataPtr(), 1);
	}

	double Nrm2(Vector<double>& x)
	{
		return cblas_dnrm2(x.Length(), x.DataPtr(), 1);
	}

	float Nrm2(Vector<ComplexFloat>& x)
	{
		return cblas_scnrm2(x.Length(), x.DataPtr(), 1);
	}

	double Nrm2(Vector<ComplexDouble>& x)
	{
		return cblas_dznrm2(x.Length(), x.DataPtr(), 1);
	}






	// ?ROT

	void Rot(Vector<float>& x, Vector<float>& y, const float c, const float s)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_srot(x.Length(), x.Data(), 1, y.Data(), 1, c, s);
	}

	void Rot(Vector<double>& x, Vector<double>& y, const double c, const double s)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_drot(x.Length(), x.Data(), 1, y.Data(), 1, c, s);
	}

	void Rot(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y, const float c, const float s)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_csrot(x.Length(), x.Data(), 1, y.Data(), 1, c, s);
	}

	void Rot(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, const double c, const double s)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_zdrot(x.Length(), x.Data(), 1, y.Data(), 1, c, s);
	}





	// ?ROTG

	void Rotg(float& a, float& b, float& c, float& s)
	{
		cblas_srotg(&a, &b, &c, &s);
	}

	void Rotg(double& a, double& b, double& c, double& s)
	{
		cblas_drotg(&a, &b, &c, &s);
	}

	void Rotg(ComplexFloat& a, ComplexFloat& b, float& c, ComplexFloat& s)
	{
		cblas_crotg(&a, &b, &c, &s);
	}

	void Rotg(ComplexDouble& a, ComplexDouble& b, double& c, ComplexDouble& s)
	{
		cblas_zrotg(&a, &b, &c, &s);
	}




	// ?ROTM

	void Rotm(Vector<float>& x, Vector<float>& y, float* param)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_srotm(x.Length(), x.Data(), 1, y.Data(), 1, param);
	}

	void Rotm(Vector<double>& x, Vector<double>& y, double* param)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_drotm(x.Length(), x.Data(), 1, y.Data(), 1, param);
	}




	// ?ROTMG
	void Rotmg(float& d1, float& d2, float& b1, float& b2, float* param)
	{
		cblas_srotmg(&d1, &d2, &b1, &b2, param);
	}

	void Rotmg(double& d1, double& d2, double& b1, double& b2, double* param)
	{
		cblas_drotmg(&d1, &d2, &b1, &b2, param);
	}






	void Scal(float alpha, Vector<float>& x)
	{
		cblas_sscal(x.Length(), alpha, x.Data(), 1);
	}

	void Scal(double alpha, Vector<double>& x)
	{
		cblas_dscal(x.Length(), alpha, x.Data(), 1);
	}

	void Scal(float alpha, Vector<ComplexFloat>& x)
	{
		cblas_csscal(x.Length(), alpha, x.Data(), 1);
	}

	void Scal(double alpha, Vector<ComplexDouble>& x)
	{
		cblas_zdscal(x.Length(), alpha, x.Data(), 1);
	}

	void Scal(ComplexFloat alpha, Vector<ComplexFloat>& x)
	{
		cblas_cscal(x.Length(), &alpha, x.Data(), 1);
	}

	void Scal(ComplexDouble alpha, Vector<ComplexDouble>& x)
	{
		cblas_zscal(x.Length(), &alpha, x.Data(), 1);
	}







	// ?SWAP

	void Swap(Vector<float>& x, Vector<float>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_sswap(x.Length(), x.Data(), 1, y.Data(), 1);
	}


	void Swap(Vector<double>& x, Vector<double>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_dswap(x.Length(), x.Data(), 1, y.Data(), 1);
	}

	void Swap(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_cswap(x.Length(), x.Data(), 1, y.Data(), 1);
	}

	void Swap(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
	{
		if(x.Length() != y.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Vector lengths are not the same!");
		}
		cblas_zswap(x.Length(), x.Data(), 1, y.Data(), 1);
	}





	// I?AMAX

	int IAmax(Vector<float>& x)
	{
		return (int)cblas_isamax(x.Length(), x.DataPtr(), 1);
	}

	int IAmax(Vector<double>& x)
	{
		return (int)cblas_idamax(x.Length(), x.DataPtr(), 1);
	}

	int IAmax(Vector<ComplexFloat>& x)
	{
		return (int)cblas_icamax(x.Length(), x.DataPtr(), 1);
	}

	int IAmax(Vector<ComplexDouble>& x)
	{
		return (int)cblas_izamax(x.Length(), x.DataPtr(), 1);
	}




	// I?AMIN

	int IAmin(Vector<float>& x)
	{
		return (int)cblas_isamin(x.Length(), x.DataPtr(), 1);
	}

	int IAmin(Vector<double>& x)
	{
		return (int)cblas_idamin(x.Length(), x.DataPtr(), 1);
	}

	int IAmin(Vector<ComplexFloat>& x)
	{
		return (int)cblas_icamin(x.Length(), x.DataPtr(), 1);
	}

	int IAmin(Vector<ComplexDouble>& x)
	{
		return (int)cblas_izamin(x.Length(), x.DataPtr(), 1);
	}

};