www.pudn.com > Image_segment.rar > Blas2.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 
{ 
 
 
	/// m*v for general matrix 
	Vector Gemv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y(A.Rows(), 0); 
		cblas_sgemv(CblasColMajor, CblasNoTrans, A.Rows(), A.Columns(), 1, A.Data(), A.Rows(), x.Data(), 1, 1, y.Data(), 1); 
		return y; 
	} 
 
	Vector Gemv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
 		Vector y(A.Rows(), 0); 
		cblas_dgemv(CblasColMajor, CblasNoTrans, A.Rows(), A.Columns(), 1, A.Data(), A.Rows(), x.Data(), 1, 1, y.Data(), 1); 
		return y; 
	} 
 
	Vector Gemv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		ComplexFloat temp(1,0); 
		Vector y(A.Rows(), ComplexFloat(0,0)); 
		cblas_cgemv(CblasColMajor, CblasNoTrans, A.Rows(), A.Columns(), &temp, A.Data(), A.Rows(), x.Data(), 1, &temp, y.Data(), 1); 
		return y; 
	} 
 
	Vector Gemv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		ComplexDouble temp(1,0); 
		Vector y(A.Rows(), ComplexDouble(0,0)); 
		cblas_zgemv(CblasColMajor, CblasNoTrans, A.Rows(), A.Columns(), &temp, A.Data(), A.Rows(), x.Data(), 1, &temp, y.Data(), 1); 
		return y; 
	} 
 
 
	/// Rank one update, general matrix 
	/// A += alpha*x*transp(y) 
	/// A += alpha*x*conjug_transp(y) if conjugate == true 
	void Ger(float alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
		cblas_sger(CblasColMajor, A.Rows(), A.Columns(), alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Ger(double alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
		cblas_dger(CblasColMajor, A.Rows(), A.Columns(), alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Ger(ComplexFloat alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		cblas_cgeru(CblasColMajor, A.Rows(), A.Columns(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Ger(ComplexDouble alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		cblas_zgeru(CblasColMajor, A.Rows(), A.Columns(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Ger(ComplexFloat alpha, Vector x, Vector y, Matrix& A, bool conjugated) 
	{ 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		if(conjugated) 
		{ 
			cblas_cgerc(CblasColMajor, A.Rows(), A.Columns(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
		} 
		else 
		{ 
			cblas_cgeru(CblasColMajor, A.Rows(), A.Columns(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
		} 
	} 
 
	void Ger(ComplexDouble alpha, Vector& x, Vector& y, Matrix& A, bool conjugated) 
	{ 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		if(conjugated) 
		{ 
			cblas_zgerc(CblasColMajor, A.Rows(), A.Columns(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
		} 
		else 
		{ 
			cblas_zgeru(CblasColMajor, A.Rows(), A.Columns(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
		} 
	} 
 
 
 
	/// Symv, Hemv: m*v for symmetric or hermitian matrix 
	Vector Symv(Matrix& A, Vector x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y(A.Rows(), 0); 
		cblas_ssymv(CblasColMajor, CblasUpper, A.Rows(), 1, A.Data(), A.Rows(), x.Data(), 1, 1, y.Data(), 1); 
		return y; 
	} 
 
	Vector Symv(Matrix& A, Vector x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y(A.Rows(), 0); 
		cblas_dsymv(CblasColMajor, CblasUpper, A.Rows(), 1, A.Data(), A.Rows(), x.Data(), 1, 1, y.Data(), 1); 
		return y; 
	} 
 
	Vector Hemv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		ComplexFloat temp(1,0); 
		Vector y(A.Rows(), ComplexFloat(0,0)); 
		cblas_chemv(CblasColMajor, CblasUpper, A.Rows(), &temp, A.Data(), A.Rows(), x.Data(), 1, &temp, y.Data(), 1); 
		return y; 
	} 
 
	Vector Hemv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		ComplexDouble temp(1,0); 
		Vector y(A.Rows(), ComplexDouble(0,0)); 
		cblas_zhemv(CblasColMajor, CblasUpper, A.Rows(), &temp, A.Data(), A.Rows(), x.Data(), 1, &temp, y.Data(), 1); 
		return y; 
	} 
 
 
 
 
 
	/// Rank-one update of a symmetric or hermitian matrix 
	/// A += alpha*x*transp(x) 
	/// A += alpha*x*conjug_transp(x) 
	void Syr(float alpha, Vector& x, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
		cblas_ssyr(CblasColMajor, CblasUpper, A.Rows(), alpha, x.Data(), 1, A.Data(), A.Rows()); 
	} 
 
 
	void Syr(double alpha, Vector& x, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
		cblas_dsyr(CblasColMajor, CblasUpper, A.Rows(), alpha, x.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Her(float alpha, Vector& x, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
		cblas_cher(CblasColMajor, CblasUpper, A.Rows(), alpha, x.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Her(double alpha, Vector& x, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
		cblas_zher(CblasColMajor, CblasUpper, A.Rows(), alpha, x.Data(), 1, A.Data(), A.Rows()); 
	} 
 
 
 
 
	/// Rank-two update of a symmetric or hermitian matrix 
	/// A += alpha*x*transp(y) + alpha*y*transp(x) 
	/// A += alpha*x*conjug_transp(y) + alpha*y*conjug_transp(x) 
	void Syr2(float alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		cblas_ssyr2(CblasColMajor, CblasUpper, A.Rows(), alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Syr2(double alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		cblas_dsyr2(CblasColMajor, CblasUpper, A.Rows(), alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Her2(ComplexFloat alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		cblas_cher2(CblasColMajor, CblasUpper, A.Rows(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
	void Her2(ComplexDouble alpha, Vector& x, Vector& y, Matrix& A) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(x.Length() != A.Rows() || y.Length() != A.Columns()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Some of the Vector lengths does not match matrix dimensions!"); 
		} 
 
		cblas_zher2(CblasColMajor, CblasUpper, A.Rows(), &alpha, x.Data(), 1, y.Data(), 1, A.Data(), A.Rows()); 
	} 
 
 
 
 
	/// Trmv: m*v for triangular matrix 
	Vector Trmv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y = x.Clone(); 
		cblas_strmv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), y.Data(), 1); 
		return y; 
	} 
 
	Vector Trmv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y = x.Clone(); 
		cblas_dtrmv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), y.Data(), 1); 
		return y; 
	} 
 
	Vector Trmv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y = x.Clone(); 
		cblas_ctrmv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), y.Data(), 1); 
		return y; 
	} 
 
	Vector Trmv(Matrix& A, Vector& x) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Columns() != x.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector y = x.Clone(); 
		cblas_ztrmv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), y.Data(), 1); 
		return y; 
	} 
 
 
	/// Trsv: Solver of a system of linear equations with a triangular matrix 
	Vector Trsv(Matrix& A, Vector& b) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Rows() != b.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector x = b.Clone(); 
		cblas_strsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), x.Data(), 1); 
		return x; 
	} 
 
	Vector Trsv(Matrix& A, Vector& b) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Rows() != b.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector x = b.Clone(); 
		cblas_dtrsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), x.Data(), 1); 
		return x; 
	} 
 
	Vector Trsv(Matrix& A, Vector& b) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Rows() != b.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector x = b.Clone(); 
		cblas_ctrsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), x.Data(), 1); 
		return x; 
	} 
 
	Vector Trsv(Matrix& A, Vector& b) 
	{ 
		if(A.Columns() != A.Rows()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix is not square!"); 
		} 
		if(A.Rows() != b.Length()) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Matrix columns and Vector length are not the same!"); 
		} 
		Vector x = b.Clone(); 
		cblas_ztrsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, A.Rows(), A.Data(), A.Rows(), x.Data(), 1); 
		return x; 
	} 
 
 
};