lapackdriver.cpp

图像分割算法

			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<ComplexFloat> Lapack::Sysv(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();
	int bs = 64*lurows;
	Vector<ComplexFloat> work(bs);
	CSYSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex8*>(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 CSYSV 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<ComplexDouble> Lapack::Sysv(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();
	int bs = 64*lurows;
	Vector<ComplexDouble> work(bs);
	ZSYSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex16*>(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 ZSYSV 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;
}

//Sysvx




//Hesv
Matrix<ComplexFloat> Lapack::Hesv(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();
	int bs = 64*lurows;
	Vector<ComplexFloat> work(bs);
	CHESV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex8*>(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 CHESV 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<ComplexDouble> Lapack::Hesv(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();
	int bs = 64*lurows;
	Vector<ComplexDouble> work(bs);
	ZHESV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex16*>(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 ZHESV 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;
}


//Hesvx















//=========================================================
// Linear Least Squares problems
//=========================================================



// Gels
Matrix<float> Lapack::Gels(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;
	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;
	max = (max>nrhs)?max:nrhs;
	int bs = min+64*max;
	Vector<float> work(bs);
	SGELS("N", &m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, 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 SGELS 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 SGELS!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<double> Lapack::Gels(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;
	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;
	max = (max>nrhs)?max:nrhs;
	int bs = min+64*max;
	Vector<double> work(bs);
	DGELS("N", &m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, 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 DGELS 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 DGELS!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexFloat> Lapack::Gels(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;
	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;
	max = (max>nrhs)?max:nrhs;
	int bs = min+64*max;
	Vector<ComplexFloat> work(bs);
	CGELS("N", &m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, reinterpret_cast<MKL_Complex8*>(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 CGELS 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 CGELS!");
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexDouble> Lapack::Gels(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;
	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;
	max = (max>nrhs)?max:nrhs;
	int bs = min+64*max;
	Vector<ComplexDouble> work(bs);
	ZGELS("N", &m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, reinterpret_cast<MKL_Complex16*>(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 ZGELS 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 ZGELS!");
		}
	}

	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}

	return X;
}


// Gelsy
Matrix<float> Lapack::Gelsy(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 = min+2*n+64*(n+1);
	int temp = 2*min+64*nrhs;
	bs = (bs>temp)?bs:temp;
	Vector<float> work(bs);
	Vector<int> jpvt(n);
	float rcond = max*numeric_limits<float>::epsilon();

	SGELSY(&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)
	{