lapackdriver.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 "StdAfx.h"
#include "./Lapack.h"
#include "mkl.h"



//=========================================================
// Linear equations solver
//=========================================================

// Gesv
Matrix<float> Lapack::Gesv(Matrix<float>& A, Matrix<float>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Vector<int> ipiv(B.Rows());
	Matrix<float> LU = A.Clone();
	Matrix<float> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	SGESV(&xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SGESV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<double> Lapack::Gesv(Matrix<double>& A, Matrix<double>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Vector<int> ipiv(B.Rows());
	Matrix<double> LU = A.Clone();
	Matrix<double> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	DGESV(&xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DGESV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<ComplexFloat> Lapack::Gesv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Vector<int> ipiv(B.Rows());
	Matrix<ComplexFloat> LU = A.Clone();
	Matrix<ComplexFloat> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	CGESV(&xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of CGESV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<ComplexDouble> Lapack::Gesv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Vector<int> ipiv(B.Rows());
	Matrix<ComplexDouble> LU = A.Clone();
	Matrix<ComplexDouble> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	ZGESV(&xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of ZGESV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}


//Gesvx


//Posv
Matrix<float> Lapack::Posv(Matrix<float>& A, Matrix<float>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Matrix<float> LU = A.Clone();
	Matrix<float> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	SPOSV("U", &xrows, &xcols, LU.Data(), &lurows, X.Data(), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SPOSV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<double> Lapack::Posv(Matrix<double>& A, Matrix<double>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Matrix<double> LU = A.Clone();
	Matrix<double> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	DPOSV("U", &xrows, &xcols, LU.Data(), &lurows, X.Data(), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DPOSV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<ComplexFloat> Lapack::Posv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Matrix<ComplexFloat> LU = A.Clone();
	Matrix<ComplexFloat> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	CPOSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of CPOSV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<ComplexDouble> Lapack::Posv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Matrix<ComplexDouble> LU = A.Clone();
	Matrix<ComplexDouble> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	ZPOSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, &info);
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of ZPOSV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

//Posvx


//Sysv
Matrix<float> Lapack::Sysv(Matrix<float>& A, Matrix<float>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Vector<int> ipiv(B.Rows());
	Matrix<float> LU = A.Clone();
	Matrix<float> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	int bs = 64*lurows;
	Vector<float> work(bs);
	SSYSV("U", &xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, work.Data(), &bs, &info);

	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SSYSV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
			Utility::RunTimeError(s);
		}
	}
	return X;
}

Matrix<double> Lapack::Sysv(Matrix<double>& A, Matrix<double>& B)
{
	if(A.Columns() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	if(A.Columns() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	int info;
	Vector<int> ipiv(B.Rows());
	Matrix<double> LU = A.Clone();
	Matrix<double> X = B.Clone();
	int xrows = X.Rows();
	int xcols = X.Columns();
	int lurows = LU.Rows();
	int bs = 64*lurows;
	Vector<double> work(bs);
	DSYSV("U", &xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, work.Data(), &bs, &info);

	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DSYSV is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{