www.pudn.com > Image_segment.rar > Blas3.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 "./Blas.h"
#include "mkl.h"
namespace Blas
{
/// 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)
{
if(A.Columns() != B.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
Matrix C(A.Rows(), B.Columns(), 0);
cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
return C;
}
Matrix Gemm(Matrix& A, Matrix& B)
{
if(A.Columns() != B.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
Matrix C(A.Rows(), B.Columns(), 0);
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
return C;
}
Matrix Gemm(Matrix& A, Matrix& B)
{
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 C(A.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
Matrix Gemm(Matrix& A, Matrix& B)
{
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 C(A.Rows(), B.Columns(), ComplexDouble(0,0));
cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
/// Hermitian or symmetric A
/// C = A*B or C = B*A (if orderReversed == true)
Matrix Symm(Matrix& A, Matrix& 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 C(B.Rows(), B.Columns(), 0);
cblas_ssymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
return C;
}
Matrix Symm(Matrix& A, Matrix& 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 C(B.Rows(), B.Columns(), 0);
cblas_dsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
return C;
}
Matrix Symm(Matrix& A, Matrix& 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 C(B.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_csymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
Matrix Symm(Matrix& A, Matrix& 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 C(B.Rows(), B.Columns(), ComplexDouble(0,0));
cblas_zsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
Matrix Hemm(Matrix& A, Matrix& 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 C(B.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_chemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
Matrix Hemm(Matrix& A, Matrix& 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 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;
}
Matrix Symm(Matrix& A, Matrix& B, bool orderReversed)
{
if(A.Columns() != A.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix is not square!");
}
if(orderReversed)
{
if(A.Rows() != B.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
Matrix C(B.Rows(), B.Columns(), 0);
cblas_ssymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, 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!");
}
Matrix C(B.Rows(), B.Columns(), 0);
cblas_ssymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
return C;
}
}
Matrix Symm(Matrix& A, Matrix& B, bool orderReversed)
{
if(A.Columns() != A.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix is not square!");
}
if(orderReversed)
{
if(A.Rows() != B.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
Matrix C(B.Rows(), B.Columns(), 0);
cblas_dsymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, 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!");
}
Matrix C(B.Rows(), B.Columns(), 0);
cblas_dsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
return C;
}
}
Matrix Symm(Matrix& A, Matrix& B, bool orderReversed)
{
if(A.Columns() != A.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix is not square!");
}
if(orderReversed)
{
if(A.Rows() != B.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
ComplexFloat temp(1,0);
Matrix C(B.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_csymm(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!");
}
ComplexFloat temp(1,0);
Matrix C(B.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_csymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
}
Matrix Symm(Matrix& A, Matrix& B, bool orderReversed)
{
if(A.Columns() != A.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix is not square!");
}
if(orderReversed)
{
if(A.Rows() != B.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
ComplexDouble temp(1,0);
Matrix C(B.Rows(), B.Columns(), ComplexDouble(0,0));
cblas_zsymm(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 C(B.Rows(), B.Columns(), ComplexDouble(0,0));
cblas_zsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
}
Matrix Hemm(Matrix& A, Matrix& B, bool orderReversed)
{
if(A.Columns() != A.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix is not square!");
}
if(orderReversed)
{
if(A.Rows() != B.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
ComplexFloat temp(1,0);
Matrix C(B.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_chemm(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!");
}
ComplexFloat temp(1,0);
Matrix C(B.Rows(), B.Columns(), ComplexFloat(0,0));
cblas_chemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
return C;
}
}
Matrix Hemm(Matrix& A, Matrix& B, bool orderReversed)
{
if(A.Columns() != A.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix is not square!");
}
if(orderReversed)
{
if(A.Rows() != B.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
}
ComplexDouble temp(1,0);
Matrix 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 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& A, Matrix& 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& A, Matrix& 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& A, Matrix& 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& A, Matrix& 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& A, Matrix& B, Matrix& 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& A, Matrix& B, Matrix& 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& A, Matrix& B, Matrix& 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& A, Matrix& B, Matrix& 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& A, Matrix& 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& A, Matrix& 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& A, Matrix& B, Matrix& 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& A, Matrix& B, Matrix& 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 Trmm(Matrix& A, Matrix& 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 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 Trmm(Matrix& A, Matrix& 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 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 Trmm(Matrix& A, Matrix& 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 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 Trmm(Matrix& A, Matrix& 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 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 Trsm(Matrix& A, Matrix& 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 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 Trsm(Matrix& A, Matrix& 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 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 Trsm(Matrix& A, Matrix& 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 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 Trsm(Matrix& A, Matrix& 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 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;
}
};