lapackdriver.cpp

图像分割算法

		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SGELSY is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Something went wrong when running SGELSY!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}



Matrix<double> Lapack::Gelsy(Matrix<double>& A, Matrix<double>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<double> F = A.Clone();
	Matrix<double> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<double>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = min+2*n+64*(n+1);
	int temp = 2*min+64*nrhs;
	bs = (bs>temp)?bs:temp;
	Vector<double> work(bs);
	Vector<int> jpvt(n);
	double rcond = max*numeric_limits<double>::epsilon();

	DGELSY(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, jpvt.Data(), &rcond, &rank, work.Data(), &bs, &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DGELSY is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Something went wrong when running DGELSY!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexFloat> Lapack::Gelsy(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<ComplexFloat> F = A.Clone();
	Matrix<ComplexFloat> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = 2*min;
	int temp = 64*(n+1);
	bs = (bs>temp)?bs:temp;
	temp = min+64*min;
	bs = (bs>temp)?bs:temp;
	temp = min+64*nrhs;
	bs = (bs>temp)?bs:temp;
	bs += min;
	Vector<ComplexFloat> work(bs);
	Vector<float> rwork(2*n);
	Vector<int> jpvt(n);
	float rcond = max*numeric_limits<float>::epsilon();

	CGELSY(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, jpvt.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of CGELSY is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Something went wrong when running CGELSY!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexDouble> Lapack::Gelsy(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<ComplexDouble> F = A.Clone();
	Matrix<ComplexDouble> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = 2*min;
	int temp = 64*(n+1);
	bs = (bs>temp)?bs:temp;
	temp = min+64*min;
	bs = (bs>temp)?bs:temp;
	temp = min+64*nrhs;
	bs = (bs>temp)?bs:temp;
	bs += min;
	Vector<ComplexDouble> work(bs);
	Vector<double> rwork(2*n);
	Vector<int> jpvt(n);
	double rcond = max*numeric_limits<double>::epsilon();

	ZGELSY(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, jpvt.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of ZGELSY is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Something went wrong when running ZGELSY!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}


// Gelss
Matrix<float> Lapack::Gelss(Matrix<float>& A, Matrix<float>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<float> F = A.Clone();
	Matrix<float> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<float>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = max;
	int temp = 2*min;
	bs = (bs>temp)?bs:temp;
	bs = (bs>nrhs)?bs:nrhs;
	bs += 3*min;
	Vector<float> work(bs);
	Vector<float> s(min);
	float rcond = max*numeric_limits<float>::epsilon();

	SGELSS(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SGELSS is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<double> Lapack::Gelss(Matrix<double>& A, Matrix<double>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<double> F = A.Clone();
	Matrix<double> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<double>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = max;
	int temp = 2*min;
	bs = (bs>temp)?bs:temp;
	bs = (bs>nrhs)?bs:nrhs;
	bs += 3*min;
	Vector<double> work(bs);
	Vector<double> s(min);
	double rcond = max*numeric_limits<double>::epsilon();

	DGELSS(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DGELSS is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexFloat> Lapack::Gelss(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<ComplexFloat> F = A.Clone();
	Matrix<ComplexFloat> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = (max>nrhs)?max:nrhs;
	bs += 2*min;
	Vector<ComplexFloat> work(bs);
	Vector<float> s(min);
	Vector<float> rwork(5*min);
	float rcond = max*numeric_limits<float>::epsilon();

	CGELSS(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of CGELSS is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexDouble> Lapack::Gelss(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<ComplexDouble> F = A.Clone();
	Matrix<ComplexDouble> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int bs = (max>nrhs)?max:nrhs;
	bs += 2*min;
	Vector<ComplexDouble> work(bs);
	Vector<double> s(min);
	Vector<double> rwork(5*min);
	double rcond = max*numeric_limits<double>::epsilon();

	ZGELSS(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);