blas3.cpp

图像分割算法

				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexDouble temp(1,0);
			Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
			cblas_zhemm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
		else
		{
			if(A.Columns() != B.Rows())
			{
				cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
				Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
			}
			ComplexDouble temp(1,0);
			Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
			cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
			return C;
		}
	}



	/// 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)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_ssyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
	}

	void Syrk(double alpha, Matrix<double>& A, Matrix<double>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
	}

	void Syrk(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		ComplexFloat temp(1,0);
		cblas_csyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), &temp, C.Data(), C.Rows());
	}

	void Syrk(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		ComplexDouble temp(1,0);
		cblas_zsyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), &temp, C.Data(), C.Rows());
	}



	/// 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)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_ssyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
	}

	void Syr2k(double alpha, Matrix<double>& A, Matrix<double>& B, Matrix<double>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_dsyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
	}

	void Syr2k(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		ComplexFloat temp(1,0);
		cblas_csyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
	}

	void Syr2k(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		ComplexDouble temp(1,0);
		cblas_zsyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
	}



	// 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)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_cherk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
	}

	void Herk(double alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
	}


	// 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)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_cher2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
	}

	void Her2k(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& C)
	{
		if(C.Columns() != C.Rows())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix is not square!");
		}
		if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match!");
		}
		cblas_zher2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
	}



	/// Triangular A
	///  C = A*B
	Matrix<float> Trmm(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 for matrix multiplication!");
		}
		Matrix<float> C = B.Clone();
		cblas_strmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), 1, A.Data(), A.Rows(), C.Data(), C.Rows());
		return C;
	}

	Matrix<double> Trmm(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 for matrix multiplication!");
		}
		Matrix<double> C = B.Clone();
		cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), 1, A.Data(), A.Rows(), C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexFloat> Trmm(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 for matrix multiplication!");
		}
		ComplexFloat temp(1,0);
		Matrix<ComplexFloat> C = B.Clone();
		cblas_ctrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), &temp, A.Data(), A.Rows(), C.Data(), C.Rows());
		return C;
	}

	Matrix<ComplexDouble> Trmm(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 for matrix multiplication!");
		}
		ComplexDouble temp(1,0);
		Matrix<ComplexDouble> C = B.Clone();
		cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), &temp, A.Data(), A.Rows(), C.Data(), C.Rows());
		return C;
	}





	// 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)
	{
		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 for matrix multiplication!");
		}
		Matrix<float> X = B.Clone();
		cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), 1, A.Data(), A.Rows(), X.Data(), X.Rows());
		return X;
	}

	Matrix<double> Trsm(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 for matrix multiplication!");
		}
		Matrix<double> X = B.Clone();
		cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), 1, A.Data(), A.Rows(), X.Data(), X.Rows());
		return X;
	}

	Matrix<ComplexFloat> Trsm(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 for matrix multiplication!");
		}
		ComplexFloat temp(1,0);
		Matrix<ComplexFloat> X = B.Clone();
		cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), &temp, A.Data(), A.Rows(), X.Data(), X.Rows());
		return X;
	}

	Matrix<ComplexDouble> Trsm(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 for matrix multiplication!");
		}
		ComplexDouble temp(1,0);
		Matrix<ComplexDouble> X = B.Clone();
		cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), &temp, A.Data(), A.Rows(), X.Data(), X.Rows());
		return X;
	}




};