//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<float> Gemv(Matrix<float>& A, Vector<float>& x)
	{
		if(A.Columns() != x.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix columns and Vector length are not the same!");
		}
		Vector<float> 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<double> Gemv(Matrix<double>& A, Vector<double>& x)
	{
		if(A.Columns() != x.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix columns and Vector length are not the same!");
		}
 		Vector<double> 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<ComplexFloat> Gemv(Matrix<ComplexFloat>& A, Vector<ComplexFloat>& 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<ComplexFloat> 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<ComplexDouble> Gemv(Matrix<ComplexDouble>& A, Vector<ComplexDouble>& 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<ComplexDouble> 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<float>& x, Vector<float>& y, Matrix<float>& 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<double>& x, Vector<double>& y, Matrix<double>& 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<ComplexFloat>& x, Vector<ComplexFloat>& y, Matrix<ComplexFloat>& 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<ComplexDouble>& x, Vector<ComplexDouble>& y, Matrix<ComplexDouble>& 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<ComplexFloat> x, Vector<ComplexFloat> y, Matrix<ComplexFloat>& 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<ComplexDouble>& x, Vector<ComplexDouble>& y, Matrix<ComplexDouble>& 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<float> Symv(Matrix<float>& A, Vector<float> 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<float> 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<double> Symv(Matrix<double>& A, Vector<double> 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<double> 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<ComplexFloat> Hemv(Matrix<ComplexFloat>& A, Vector<ComplexFloat>& 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<ComplexFloat> 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<ComplexDouble> Hemv(Matrix<ComplexDouble>& A, Vector<ComplexDouble>& 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<ComplexDouble> 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<float>& x, Matrix<float>& A)
	{
		if(A.Columns() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");