blas.h

图像分割算法

//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.



#pragma once
#ifndef BLAS_H
#define BLAS_H

#include "cimpl.h"
using namespace CIMPL;


#include <cmath>



/// \brief Blas Toolbox
namespace Blas
{

// BLAS LEVEL 1

	/// |x(0)|+...+|x(n)|
	 float Asum(Vector<float>& x);
	 double Asum(Vector<double>& x);
	 float Asum(Vector<ComplexFloat>& x);
	 double Asum(Vector<ComplexDouble>& x);

	/// y += alpha*x
	 void Axpy(float alpha, Vector<float>& x, Vector<float>& y);
	 void Axpy(double alpha, Vector<double>& x, Vector<double>& y);
	 void Axpy(ComplexFloat alpha, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y);
	 void Axpy(ComplexDouble alpha, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y);

	/// y = x
	 void Copy(Vector<float>& x, Vector<float>& y);
	 void Copy(Vector<double>& x, Vector<double>& y);
	 void Copy(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y);
	 void Copy(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y);

	/// Dot product of x and y
	 float Dot(Vector<float>& x, Vector<float>& y);
	 double Dot(Vector<double>& x, Vector<double>& y);
	 ComplexFloat Dot(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y, bool conjugated);
	 ComplexDouble Dot(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, bool conjugated);
	 ComplexFloat Dot(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y);
	 ComplexDouble Dot(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y);

	/// Calculates the Euclidean norm of x.
	/// Euclidean norm is the square	root of	the conjugated dot product of a	vector with itself.
	 float Nrm2(Vector<float>& x);
	 double Nrm2(Vector<double>& x);
	 float Nrm2(Vector<ComplexFloat>& x);
	 double Nrm2(Vector<ComplexDouble>& x);

	/// Apply givens plane rotation
	/// x(i) = c*x(i) + s*y(i)     y(i) = -s*x(i) + c*y(i)
	 void Rot(Vector<float>& x, Vector<float>& y, const float c, const float s);
	 void Rot(Vector<double>& x, Vector<double>& y, const double c, const double s);
	 void Rot(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y, const float c, const float s);
	 void Rot(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, const double c, const double s);


	/// Generate elements for a givens plane rotation
	 void Rotg(float& a, float& b, float& c, float& s);
	 void Rotg(double& a, double& b, double& c, double& s);
	 void Rotg(ComplexFloat& a, ComplexFloat& b, float& c, ComplexFloat& s);
	 void Rotg(ComplexDouble& a, ComplexDouble& b, double& c, ComplexDouble& s);



	/// \brief Apply modified givens transformation
	///
	///| x(i) |       | x(i) |
	///|      | = H * |	   |
	///| y(i) |       | y(i) |
	///
	///Depending on the value of PARAM(1), the transformation matrix	is defined as
	///follows:
	///
	///PARAM(1)= -1.0
	///H(11)	H(12)
	///H(21)	H(22)
	///
	///PARAM(1)= 0.0
	///1.0	H(12)
	///H(21)	  1.0
	///
	///PARAM(1)= 1.0
	///H(11)	  1.0
	///-1.0	H(22)
	///
	///PARAM(1)= -2.0
	///1.0	  0.0
	///0.0	  1.0
	///
	///The array PARAM is generated by a call to the	routine	_ROTMG.
	 void Rotm(Vector<float>& x, Vector<float>& y, float* param);
	 void Rotm(Vector<double>& x, Vector<double>& y, double* param);


	/// Generate elements for a modified Givens transform
	 void Rotmg(float& d1, float& d2, float& b1, float& b2, float* param);
	 void Rotmg(double& d1, double& d2, double& b1, double& b2, double* param);


	/// x *= alpha
	 void Scal(float alpha, Vector<float>& x);
	 void Scal(double alpha, Vector<double>& x);
	 void Scal(float alpha, Vector<ComplexFloat>& x);
	 void Scal(double alpha, Vector<ComplexDouble>& x);
	 void Scal(ComplexFloat alpha, Vector<ComplexFloat>& x);
	 void Scal(ComplexDouble alpha, Vector<ComplexDouble>& x);

	/// Swap elements of two vectors
	 void Swap(Vector<float>& x, Vector<float>& y);
	 void Swap(Vector<double>& x, Vector<double>& y);
	 void Swap(Vector<ComplexFloat>& x, Vector<ComplexFloat>& y);
	 void Swap(Vector<ComplexDouble>& x, Vector<ComplexDouble>& y);


	/// index of the maximum element (First element is indexed as 1)
	 int IAmax(Vector<float>& x);
	 int IAmax(Vector<double>& x);
	 int IAmax(Vector<ComplexFloat>& x);
	 int IAmax(Vector<ComplexDouble>& x);


	/// index of the minimum element (First element is indexed as 1)
	 int IAmin(Vector<float>& x);
	 int IAmin(Vector<double>& x);
	 int IAmin(Vector<ComplexFloat>& x);
	 int IAmin(Vector<ComplexDouble>& x);



// BLAS LEVEL 2


	/// m*v for general matrix
	 Vector<float> Gemv(Matrix<float>& A, Vector<float>& x);
	 Vector<double> Gemv(Matrix<double>& A, Vector<double>& x);
	 Vector<ComplexFloat> Gemv(Matrix<ComplexFloat>& A, Vector<ComplexFloat>& x);
	 Vector<ComplexDouble> Gemv(Matrix<ComplexDouble>& A, Vector<ComplexDouble>& x);

	/// 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);
	 void Ger(double alpha, Vector<double>& x, Vector<double>& y, Matrix<double>& A);
	 void Ger(ComplexFloat alpha, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y, Matrix<ComplexFloat>& A);
	 void Ger(ComplexDouble alpha, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, Matrix<ComplexDouble>& A);
	 void Ger(ComplexFloat alpha, Vector<ComplexFloat> x, Vector<ComplexFloat> y, Matrix<ComplexFloat>& A, bool conjugated);
	 void Ger(ComplexDouble alpha, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, Matrix<ComplexDouble>& A, bool conjugated);
	
	
	/// Symv, Hemv: m*v for symmetric or hermitian matrix
	 Vector<float> Symv(Matrix<float>& A, Vector<float> x);
	 Vector<double> Symv(Matrix<double>& A, Vector<double> x);
	 
	 Vector<ComplexFloat> Hemv(Matrix<ComplexFloat>& A, Vector<ComplexFloat>& x);
	 Vector<ComplexDouble> Hemv(Matrix<ComplexDouble>& A, Vector<ComplexDouble>& x);

	/// Syr, Her: Rank one update, symmetric or hermitian matrix
	/// 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);
	 void Syr(double alpha, Vector<double>& x, Matrix<double>& A);
	 
	 void Her(float alpha, Vector<ComplexFloat>& x, Matrix<ComplexFloat>& A);
	 void Her(double alpha, Vector<ComplexDouble>& x, Matrix<ComplexDouble>& A);


	/// Syr2, Her2: Rank two update, symmetric or hermitian matrix
	/// Rank-two update of a symmetric or hermitian matrix
	/// A += alpha*x*transp(y) + alpha*y*transp(x)
	/// A += alpha*x*conjug_transp(x) + alpha*x*conjug_transp(y)
	 void Syr2(float alpha, Vector<float>& x, Vector<float>& y, Matrix<float>& A);
	 void Syr2(double alpha, Vector<double>& x, Vector<double>& y, Matrix<double>& A);
	 
	 void Her2(ComplexFloat alpha, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y, Matrix<ComplexFloat>& A);
	 void Her2(ComplexDouble alpha, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y, Matrix<ComplexDouble>& A);



	/// Trmv: m*v for triangular matrix
	 Vector<float> Trmv(Matrix<float>& A, Vector<float>& x);
	 Vector<double> Trmv(Matrix<double>& A, Vector<double>& x);
	 Vector<ComplexFloat> Trmv(Matrix<ComplexFloat>& A, Vector<ComplexFloat>& x);
	 Vector<ComplexDouble> Trmv(Matrix<ComplexDouble>& A, Vector<ComplexDouble>& x);
	
	/// Trsv: Solver of a system of linear equations with a triangular matrix
	 Vector<float> Trsv(Matrix<float>& A, Vector<float>& b);
	 Vector<double> Trsv(Matrix<double>& A, Vector<double>& b);
	 Vector<ComplexFloat> Trsv(Matrix<ComplexFloat>& A, Vector<ComplexFloat>& b);
	 Vector<ComplexDouble> Trsv(Matrix<ComplexDouble>& A, Vector<ComplexDouble>& b);
	







// BLAS LEVEL 3

	/// 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);
	 Matrix<double> Gemm(Matrix<double>& A, Matrix<double>& B);
	 Matrix<ComplexFloat> Gemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B);
	 Matrix<ComplexDouble> Gemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B);

	
	// Hemm, Symm: Hermitian, symmetric (A is sym. or hermitian)  m*m
	///  C = A*B or C = B*A
	 Matrix<float> Symm(Matrix<float>& A, Matrix<float>& B);
	 Matrix<double> Symm(Matrix<double>& A, Matrix<double>& B);
	 Matrix<ComplexFloat> Symm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B);
	 Matrix<ComplexDouble> Symm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B);
	 

	 Matrix<float> Symm(Matrix<float>& A, Matrix<float>& B, bool orderReversed);
	 Matrix<double> Symm(Matrix<double>& A, Matrix<double>& B, bool orderReversed);
	 Matrix<ComplexFloat> Symm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, bool orderReversed);
	 Matrix<ComplexDouble> Symm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, bool orderReversed);
	 
	 Matrix<ComplexFloat> Hemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B);
	 Matrix<ComplexDouble> Hemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B);
	 Matrix<ComplexFloat> Hemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, bool orderReversed);
	 Matrix<ComplexDouble> Hemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, bool orderReversed);


	/// Rank k update of a symmetric matrix
	/// C += alpha*A*transp(A)
	/// A is of size n by k
	
	 void Syrk(float alpha, Matrix<float>& A, Matrix<float>& C);
	 void Syrk(double alpha, Matrix<double>& A, Matrix<double>& C);
	 void Syrk(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& C);
	 void Syrk(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& C);


	/// rank 2k update of a symmetric matrix
	/// C += alpha*A*transp(B) + alpha*B*transp(A)
	/// A and B are of size n by k

	 void Syr2k(float alpha, Matrix<float>& A, Matrix<float>& B, Matrix<float>& C);
	 void Syr2k(double alpha, Matrix<double>& A, Matrix<double>& B, Matrix<double>& C);
	 void Syr2k(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& C);
	 void Syr2k(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& C);


	// Herk: rank k update of an hermitian matrix
	/// C += alpha*A*conjug_transp(A)
	/// A is of size n by k

	 void Herk(float alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& C);
	 void Herk(double alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& C);

	// Her2k: rank 2k update of an hermitian matrix
	/// C += alpha*A*conjug_transp(B) + alpha*B*conjug_transp(A)
	/// A and B are of size n by k

	 void Her2k(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& C);
	 void Her2k(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& C);


	// Trmm: Triangular (A is triangular) m*m
	///  C = A*B
	 Matrix<float> Trmm(Matrix<float>& A, Matrix<float>& B);
	 Matrix<double> Trmm(Matrix<double>& A, Matrix<double>& B);
	 Matrix<ComplexFloat> Trmm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B);
	 Matrix<ComplexDouble> Trmm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B);

	// Trsm: Solve a triangular system of equations with a triangular coefficient matrix
	/// Solve AX=B for X
	 Matrix<float> Trsm(Matrix<float>& A, Matrix<float>& B);
	 Matrix<double> Trsm(Matrix<double>& A, Matrix<double>& B);
	 Matrix<ComplexFloat> Trsm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B);
	 Matrix<ComplexDouble> Trsm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B);


};




#endif