www.pudn.com > Image_segment.rar > 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  
 
 
 
/// \brief Blas Toolbox 
namespace Blas 
{ 
 
// BLAS LEVEL 1 
 
	/// |x(0)|+...+|x(n)| 
	 float Asum(Vector& x); 
	 double Asum(Vector& x); 
	 float Asum(Vector& x); 
	 double Asum(Vector& x); 
 
	/// y += alpha*x 
	 void Axpy(float alpha, Vector& x, Vector& y); 
	 void Axpy(double alpha, Vector& x, Vector& y); 
	 void Axpy(ComplexFloat alpha, Vector& x, Vector& y); 
	 void Axpy(ComplexDouble alpha, Vector& x, Vector& y); 
 
	/// y = x 
	 void Copy(Vector& x, Vector& y); 
	 void Copy(Vector& x, Vector& y); 
	 void Copy(Vector& x, Vector& y); 
	 void Copy(Vector& x, Vector& y); 
 
	/// Dot product of x and y 
	 float Dot(Vector& x, Vector& y); 
	 double Dot(Vector& x, Vector& y); 
	 ComplexFloat Dot(Vector& x, Vector& y, bool conjugated); 
	 ComplexDouble Dot(Vector& x, Vector& y, bool conjugated); 
	 ComplexFloat Dot(Vector& x, Vector& y); 
	 ComplexDouble Dot(Vector& x, Vector& 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& x); 
	 double Nrm2(Vector& x); 
	 float Nrm2(Vector& x); 
	 double Nrm2(Vector& x); 
 
	/// Apply givens plane rotation 
	/// x(i) = c*x(i) + s*y(i)     y(i) = -s*x(i) + c*y(i) 
	 void Rot(Vector& x, Vector& y, const float c, const float s); 
	 void Rot(Vector& x, Vector& y, const double c, const double s); 
	 void Rot(Vector& x, Vector& y, const float c, const float s); 
	 void Rot(Vector& x, Vector& 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& x, Vector& y, float* param); 
	 void Rotm(Vector& x, Vector& 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& x); 
	 void Scal(double alpha, Vector& x); 
	 void Scal(float alpha, Vector& x); 
	 void Scal(double alpha, Vector& x); 
	 void Scal(ComplexFloat alpha, Vector& x); 
	 void Scal(ComplexDouble alpha, Vector& x); 
 
	/// Swap elements of two vectors 
	 void Swap(Vector& x, Vector& y); 
	 void Swap(Vector& x, Vector& y); 
	 void Swap(Vector& x, Vector& y); 
	 void Swap(Vector& x, Vector& y); 
 
 
	/// index of the maximum element (First element is indexed as 1) 
	 int IAmax(Vector& x); 
	 int IAmax(Vector& x); 
	 int IAmax(Vector& x); 
	 int IAmax(Vector& x); 
 
 
	/// index of the minimum element (First element is indexed as 1) 
	 int IAmin(Vector& x); 
	 int IAmin(Vector& x); 
	 int IAmin(Vector& x); 
	 int IAmin(Vector& x); 
 
 
 
// BLAS LEVEL 2 
 
 
	/// m*v for general matrix 
	 Vector Gemv(Matrix& A, Vector& x); 
	 Vector Gemv(Matrix& A, Vector& x); 
	 Vector Gemv(Matrix& A, Vector& x); 
	 Vector Gemv(Matrix& A, Vector& 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& x, Vector& y, Matrix& A); 
	 void Ger(double alpha, Vector& x, Vector& y, Matrix& A); 
	 void Ger(ComplexFloat alpha, Vector& x, Vector& y, Matrix& A); 
	 void Ger(ComplexDouble alpha, Vector& x, Vector& y, Matrix& A); 
	 void Ger(ComplexFloat alpha, Vector x, Vector y, Matrix& A, bool conjugated); 
	 void Ger(ComplexDouble alpha, Vector& x, Vector& y, Matrix& A, bool conjugated); 
	 
	 
	/// Symv, Hemv: m*v for symmetric or hermitian matrix 
	 Vector Symv(Matrix& A, Vector x); 
	 Vector Symv(Matrix& A, Vector x); 
	  
	 Vector Hemv(Matrix& A, Vector& x); 
	 Vector Hemv(Matrix& A, Vector& 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& x, Matrix& A); 
	 void Syr(double alpha, Vector& x, Matrix& A); 
	  
	 void Her(float alpha, Vector& x, Matrix& A); 
	 void Her(double alpha, Vector& x, Matrix& 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& x, Vector& y, Matrix& A); 
	 void Syr2(double alpha, Vector& x, Vector& y, Matrix& A); 
	  
	 void Her2(ComplexFloat alpha, Vector& x, Vector& y, Matrix& A); 
	 void Her2(ComplexDouble alpha, Vector& x, Vector& y, Matrix& A); 
 
 
 
	/// Trmv: m*v for triangular matrix 
	 Vector Trmv(Matrix& A, Vector& x); 
	 Vector Trmv(Matrix& A, Vector& x); 
	 Vector Trmv(Matrix& A, Vector& x); 
	 Vector Trmv(Matrix& A, Vector& x); 
	 
	/// Trsv: Solver of a system of linear equations with a triangular matrix 
	 Vector Trsv(Matrix& A, Vector& b); 
	 Vector Trsv(Matrix& A, Vector& b); 
	 Vector Trsv(Matrix& A, Vector& b); 
	 Vector Trsv(Matrix& A, Vector& 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 Gemm(Matrix& A, Matrix& B); 
	 Matrix Gemm(Matrix& A, Matrix& B); 
	 Matrix Gemm(Matrix& A, Matrix& B); 
	 Matrix Gemm(Matrix& A, Matrix& B); 
 
	 
	// Hemm, Symm: Hermitian, symmetric (A is sym. or hermitian)  m*m 
	///  C = A*B or C = B*A 
	 Matrix Symm(Matrix& A, Matrix& B); 
	 Matrix Symm(Matrix& A, Matrix& B); 
	 Matrix Symm(Matrix& A, Matrix& B); 
	 Matrix Symm(Matrix& A, Matrix& B); 
	  
 
	 Matrix Symm(Matrix& A, Matrix& B, bool orderReversed); 
	 Matrix Symm(Matrix& A, Matrix& B, bool orderReversed); 
	 Matrix Symm(Matrix& A, Matrix& B, bool orderReversed); 
	 Matrix Symm(Matrix& A, Matrix& B, bool orderReversed); 
	  
	 Matrix Hemm(Matrix& A, Matrix& B); 
	 Matrix Hemm(Matrix& A, Matrix& B); 
	 Matrix Hemm(Matrix& A, Matrix& B, bool orderReversed); 
	 Matrix Hemm(Matrix& A, Matrix& 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& A, Matrix& C); 
	 void Syrk(double alpha, Matrix& A, Matrix& C); 
	 void Syrk(ComplexFloat alpha, Matrix& A, Matrix& C); 
	 void Syrk(ComplexDouble alpha, Matrix& A, Matrix& 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& A, Matrix& B, Matrix& C); 
	 void Syr2k(double alpha, Matrix& A, Matrix& B, Matrix& C); 
	 void Syr2k(ComplexFloat alpha, Matrix& A, Matrix& B, Matrix& C); 
	 void Syr2k(ComplexDouble alpha, Matrix& A, Matrix& B, Matrix& 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& A, Matrix& C); 
	 void Herk(double alpha, Matrix& A, Matrix& 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& A, Matrix& B, Matrix& C); 
	 void Her2k(ComplexDouble alpha, Matrix& A, Matrix& B, Matrix& C); 
 
 
	// Trmm: Triangular (A is triangular) m*m 
	///  C = A*B 
	 Matrix Trmm(Matrix& A, Matrix& B); 
	 Matrix Trmm(Matrix& A, Matrix& B); 
	 Matrix Trmm(Matrix& A, Matrix& B); 
	 Matrix Trmm(Matrix& A, Matrix& B); 
 
	// Trsm: Solve a triangular system of equations with a triangular coefficient matrix 
	/// Solve AX=B for X 
	 Matrix Trsm(Matrix& A, Matrix& B); 
	 Matrix Trsm(Matrix& A, Matrix& B); 
	 Matrix Trsm(Matrix& A, Matrix& B); 
	 Matrix Trsm(Matrix& A, Matrix& B); 
 
 
}; 
 
 
 
 
#endif