1139 lines
28 KiB
C++
1139 lines
28 KiB
C++
// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
|
||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||
//
|
||
// This file is part of the MFEM library. For more information and source code
|
||
// availability visit https://mfem.org.
|
||
//
|
||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||
// CONTRIBUTING.md for details.
|
||
|
||
#include "complex_densemat.hpp"
|
||
#include <complex>
|
||
|
||
#ifdef MFEM_USE_LAPACK
|
||
#ifdef MFEM_USE_SINGLE
|
||
extern "C" void
|
||
cgetrf_(int *, int *, std::complex<float> *, int *, int *, int *);
|
||
extern "C" void
|
||
cgetrs_(char *, int *, int *, std::complex<float> *, int *, int *,
|
||
std::complex<float> *, int *, int *);
|
||
extern "C" void
|
||
cgetri_(int *, std::complex<float> *, int *, int *,
|
||
std::complex<float> *, int *, int *);
|
||
extern "C" void
|
||
ctrsm_(char *, char *, char *, char *, int *, int *, std::complex<float> *,
|
||
std::complex<float> *, int *, std::complex<float> *, int *);
|
||
extern "C" void
|
||
cpotrf_(char *, int *, std::complex<float> *, int *, int *);
|
||
|
||
extern "C" void
|
||
ctrtrs_(char *, char*, char *, int *, int *, std::complex<float> *, int *,
|
||
std::complex<float> *, int *, int *);
|
||
extern "C" void
|
||
cpotri_(char *, int *, std::complex<float> *, int*, int *);
|
||
|
||
extern "C" void
|
||
cpotrs_(char *, int *, int *, std::complex<float> *, int *,
|
||
std::complex<float> *, int *, int *);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
extern "C" void
|
||
zgetrf_(int *, int *, std::complex<double> *, int *, int *, int *);
|
||
extern "C" void
|
||
zgetrs_(char *, int *, int *, std::complex<double> *, int *, int *,
|
||
std::complex<double> *, int *, int *);
|
||
extern "C" void
|
||
zgetri_(int *, std::complex<double> *, int *, int *,
|
||
std::complex<double> *, int *, int *);
|
||
extern "C" void
|
||
ztrsm_(char *, char *, char *, char *, int *, int *, std::complex<double> *,
|
||
std::complex<double> *, int *, std::complex<double> *, int *);
|
||
extern "C" void
|
||
zpotrf_(char *, int *, std::complex<double> *, int *, int *);
|
||
|
||
extern "C" void
|
||
ztrtrs_(char *, char*, char *, int *, int *, std::complex<double> *, int *,
|
||
std::complex<double> *, int *, int *);
|
||
extern "C" void
|
||
zpotri_(char *, int *, std::complex<double> *, int*, int *);
|
||
|
||
extern "C" void
|
||
zpotrs_(char *, int *, int *, std::complex<double> *, int *,
|
||
std::complex<double> *, int *, int *);
|
||
#endif
|
||
#endif
|
||
|
||
namespace mfem
|
||
{
|
||
|
||
DenseMatrix & ComplexDenseMatrix::real()
|
||
{
|
||
MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
|
||
return dynamic_cast<DenseMatrix &>(*Op_Real_);
|
||
}
|
||
|
||
DenseMatrix & ComplexDenseMatrix::imag()
|
||
{
|
||
MFEM_ASSERT(Op_Imag_, "ComplexDenseMatrix has no imaginary part!");
|
||
return dynamic_cast<DenseMatrix &>(*Op_Imag_);
|
||
}
|
||
|
||
const DenseMatrix & ComplexDenseMatrix::real() const
|
||
{
|
||
MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
|
||
return dynamic_cast<const DenseMatrix &>(*Op_Real_);
|
||
}
|
||
|
||
const DenseMatrix & ComplexDenseMatrix::imag() const
|
||
{
|
||
MFEM_ASSERT(Op_Imag_, "ComplexDenseMatrix has no imaginary part!");
|
||
return dynamic_cast<const DenseMatrix &>(*Op_Imag_);
|
||
}
|
||
|
||
DenseMatrix * ComplexDenseMatrix::GetSystemMatrix() const
|
||
{
|
||
int h = height/2;
|
||
int w = width/2;
|
||
DenseMatrix * A = new DenseMatrix(2*h,2*w);
|
||
real_t * data = A->Data();
|
||
real_t * data_r = nullptr;
|
||
real_t * data_i = nullptr;
|
||
|
||
const real_t factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
|
||
*A = 0.;
|
||
if (hasRealPart())
|
||
{
|
||
data_r = real().Data();
|
||
for (int j = 0; j<w; j++)
|
||
{
|
||
for (int i = 0; i<h; i++)
|
||
{
|
||
data[i+j*height] = data_r[i+j*h];
|
||
data[i+h+(j+h)*height] = factor*data_r[i+j*h];
|
||
}
|
||
}
|
||
}
|
||
if (hasImagPart())
|
||
{
|
||
data_i = imag().Data();
|
||
for (int j = 0; j<w; j++)
|
||
{
|
||
for (int i = 0; i<h; i++)
|
||
{
|
||
data[i+h+j*height] = factor*data_i[i+j*h];
|
||
data[i+(j+h)*height] = -data_i[i+j*h];
|
||
}
|
||
}
|
||
}
|
||
return A;
|
||
}
|
||
|
||
ComplexDenseMatrix * ComplexDenseMatrix::ComputeInverse()
|
||
{
|
||
MFEM_VERIFY(height == width, "Matrix has to be square");
|
||
|
||
// complex data
|
||
int h = height/2;
|
||
int w = width/2;
|
||
std::complex<real_t> * data = new std::complex<real_t>[h*w];
|
||
|
||
// copy data
|
||
if (hasRealPart() && hasImagPart())
|
||
{
|
||
real_t * data_r = real().Data();
|
||
real_t * data_i = imag().Data();
|
||
for (int i = 0; i < h*w; i++)
|
||
{
|
||
data[i] = std::complex<real_t> (data_r[i], data_i[i]);
|
||
}
|
||
}
|
||
else if (hasRealPart())
|
||
{
|
||
real_t * data_r = real().Data();
|
||
for (int i = 0; i < h*w; i++)
|
||
{
|
||
data[i] = std::complex<real_t> (data_r[i], 0.);
|
||
}
|
||
}
|
||
else if (hasImagPart())
|
||
{
|
||
real_t * data_i = imag().Data();
|
||
for (int i = 0; i < h*w; i++)
|
||
{
|
||
data[i] = std::complex<real_t> (0., data_i[i]);
|
||
}
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("ComplexDenseMatrix has neither only a real nor an imag part");
|
||
}
|
||
|
||
#ifdef MFEM_USE_LAPACK
|
||
int *ipiv = new int[w];
|
||
int lwork = -1;
|
||
std::complex<real_t> qwork, *work;
|
||
int info;
|
||
|
||
#ifdef MFEM_USE_SINGLE
|
||
cgetrf_(&w, &w, data, &w, ipiv, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
zgetrf_(&w, &w, data, &w, ipiv, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
if (info)
|
||
{
|
||
mfem_error("DenseMatrix::Invert() : Error in ZGETRF");
|
||
}
|
||
|
||
#ifdef MFEM_USE_SINGLE
|
||
cgetri_(&w, data, &w, ipiv, &qwork, &lwork, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
zgetri_(&w, data, &w, ipiv, &qwork, &lwork, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
lwork = (int) qwork.real();
|
||
work = new std::complex<real_t>[lwork];
|
||
|
||
#ifdef MFEM_USE_SINGLE
|
||
cgetri_(&w, data, &w, ipiv, work, &lwork, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
zgetri_(&w, data, &w, ipiv, work, &lwork, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
if (info)
|
||
{
|
||
mfem_error("DenseMatrix::Invert() : Error in ZGETRI");
|
||
}
|
||
|
||
delete [] work;
|
||
delete [] ipiv;
|
||
|
||
#else
|
||
// compiling without LAPACK
|
||
int c, i, j, n = w;
|
||
real_t a, b;
|
||
Array<int> piv(n);
|
||
std::complex<real_t> ac,bc;
|
||
|
||
for (c = 0; c < n; c++)
|
||
{
|
||
a = std::abs(data[c+c*h]);
|
||
i = c;
|
||
for (j = c + 1; j < n; j++)
|
||
{
|
||
b = std::abs(data[j+c*h]);
|
||
if (a < b)
|
||
{
|
||
a = b;
|
||
i = j;
|
||
}
|
||
}
|
||
if (a == 0.0)
|
||
{
|
||
mfem_error("DenseMatrix::Invert() : singular matrix");
|
||
}
|
||
piv[c] = i;
|
||
for (j = 0; j < n; j++)
|
||
{
|
||
mfem::Swap<std::complex<real_t>>(data[c+j*h], data[i+j*h]);
|
||
}
|
||
|
||
real_t fone = 1.0; // TODO: why?
|
||
ac = data[c+c*h] = fone / data[c+c*h];
|
||
for (j = 0; j < c; j++)
|
||
{
|
||
data[c+j*h] *= ac;
|
||
}
|
||
for (j++; j < n; j++)
|
||
{
|
||
data[c+j*h] *= ac;
|
||
}
|
||
for (i = 0; i < c; i++)
|
||
{
|
||
data[i+c*h] = ac * (bc = -data[i+c*h]);
|
||
for (j = 0; j < c; j++)
|
||
{
|
||
data[i+j*h] += bc * data[c+j*h];
|
||
}
|
||
for (j++; j < n; j++)
|
||
{
|
||
data[i+j*h] += bc * data[c+j*h];
|
||
}
|
||
}
|
||
for (i++; i < n; i++)
|
||
{
|
||
data[i+c*h] = ac * (bc = -data[i+c*h]);
|
||
for (j = 0; j < c; j++)
|
||
{
|
||
data[i+j*h] += bc * data[c+j*h];
|
||
}
|
||
for (j++; j < n; j++)
|
||
{
|
||
data[i+j*h] += bc * data[c+j*h];
|
||
}
|
||
}
|
||
}
|
||
|
||
for (c = n - 1; c >= 0; c--)
|
||
{
|
||
j = piv[c];
|
||
for (i = 0; i < n; i++)
|
||
{
|
||
mfem::Swap<std::complex<real_t>>(data[i+c*h], data[i+j*h]);
|
||
}
|
||
}
|
||
|
||
#endif
|
||
|
||
DenseMatrix * C_r = new DenseMatrix(h);
|
||
DenseMatrix * C_i = new DenseMatrix(h);
|
||
|
||
real_t * datac_r = C_r->Data();
|
||
real_t * datac_i = C_i->Data();
|
||
|
||
for (int k = 0; k < h*w; k++)
|
||
{
|
||
datac_r[k] = data[k].real();
|
||
datac_i[k] = data[k].imag();
|
||
}
|
||
delete [] data;
|
||
return new ComplexDenseMatrix(C_r,C_i,true,true);
|
||
|
||
}
|
||
|
||
ComplexDenseMatrix * Mult(const ComplexDenseMatrix &A,
|
||
const ComplexDenseMatrix &B)
|
||
{
|
||
// C = C_r + i C_i = (A_r + i * A_i) * (B_r + i * B_i)
|
||
// = A_r * B_r - A_i B_i + i (A_r * B_i + A_i * B_r)
|
||
|
||
int h = A.Height()/2;
|
||
int w = B.Width()/2;
|
||
|
||
MFEM_VERIFY(A.Width() == B.Height(), "Incompatible matrix dimensions");
|
||
|
||
//only real case (imag is null)
|
||
DenseMatrix * C_r = nullptr;
|
||
DenseMatrix * C_i = nullptr;
|
||
if ((A.hasRealPart() && B.hasRealPart()) ||
|
||
(A.hasImagPart() && B.hasImagPart()))
|
||
{
|
||
C_r = new DenseMatrix(h,w);
|
||
}
|
||
if ((A.hasRealPart() && B.hasImagPart()) ||
|
||
(A.hasImagPart() && B.hasRealPart()))
|
||
{
|
||
C_i = new DenseMatrix(h,w);
|
||
}
|
||
|
||
MFEM_VERIFY(C_r || C_i, "Both real and imag parts are null");
|
||
|
||
if (A.hasRealPart() && B.hasRealPart())
|
||
{
|
||
Mult(A.real(), B.real(),*C_r);
|
||
}
|
||
if (A.hasImagPart() && B.hasImagPart())
|
||
{
|
||
if (A.hasRealPart() && B.hasRealPart())
|
||
{
|
||
AddMult_a(-1.,A.imag(), B.imag(),*C_r);
|
||
}
|
||
else
|
||
{
|
||
Mult(A.imag(), B.imag(),*C_r);
|
||
}
|
||
}
|
||
|
||
if (A.hasRealPart() && B.hasImagPart())
|
||
{
|
||
Mult(A.real(), B.imag(),*C_i);
|
||
}
|
||
|
||
if (A.hasImagPart() && B.hasRealPart())
|
||
{
|
||
if (A.hasRealPart() && B.hasImagPart())
|
||
{
|
||
AddMult(A.imag(), B.real(),*C_i);
|
||
}
|
||
else
|
||
{
|
||
Mult(A.imag(), B.real(),*C_i);
|
||
}
|
||
}
|
||
|
||
return new ComplexDenseMatrix(C_r,C_i,true,true);
|
||
}
|
||
|
||
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A,
|
||
const ComplexDenseMatrix &B)
|
||
{
|
||
// C = C_r + i C_i = (A_r^t - i * A_i^t) * (B_r + i * B_i)
|
||
// = A_r^t * B_r + A_i^t * B_i + i (A_r^t * B_i - A_i^t * B_r)
|
||
|
||
int h = A.Width()/2;
|
||
int w = B.Width()/2;
|
||
|
||
MFEM_VERIFY(A.Height() == B.Height(), "Incompatible matrix dimensions");
|
||
|
||
//only real case (imag is null)
|
||
DenseMatrix * C_r = nullptr;
|
||
DenseMatrix * C_i = nullptr;
|
||
if ((A.hasRealPart() && B.hasRealPart()) ||
|
||
(A.hasImagPart() && B.hasImagPart()))
|
||
{
|
||
C_r = new DenseMatrix(h,w);
|
||
}
|
||
if ((A.hasRealPart() && B.hasImagPart()) ||
|
||
(A.hasImagPart() && B.hasRealPart()))
|
||
{
|
||
C_i = new DenseMatrix(h,w);
|
||
}
|
||
|
||
MFEM_VERIFY(C_r || C_i, "Both real and imag parts are null");
|
||
|
||
if (A.hasRealPart() && B.hasRealPart())
|
||
{
|
||
MultAtB(A.real(), B.real(),*C_r);
|
||
}
|
||
if (A.hasImagPart() && B.hasImagPart())
|
||
{
|
||
if (A.hasRealPart() && B.hasRealPart())
|
||
{
|
||
DenseMatrix tempC_r(h,w);
|
||
MultAtB(A.imag(), B.imag(),tempC_r);
|
||
(*C_r) += tempC_r;
|
||
}
|
||
else
|
||
{
|
||
MultAtB(A.imag(), B.imag(),*C_r);
|
||
}
|
||
}
|
||
|
||
if (A.hasRealPart() && B.hasImagPart())
|
||
{
|
||
MultAtB(A.real(), B.imag(),*C_i);
|
||
}
|
||
|
||
if (A.hasImagPart() && B.hasRealPart())
|
||
{
|
||
if (A.hasRealPart() && B.hasImagPart())
|
||
{
|
||
DenseMatrix tempC_i(h,w);
|
||
MultAtB(A.imag(), B.real(),tempC_i);
|
||
(*C_i) -= tempC_i;
|
||
}
|
||
else
|
||
{
|
||
MultAtB(A.imag(), B.real(),*C_i);
|
||
}
|
||
}
|
||
|
||
return new ComplexDenseMatrix(C_r,C_i,true,true);
|
||
}
|
||
|
||
std::complex<real_t> * ComplexFactors::RealToComplex
|
||
(int m, const real_t * x_r, const real_t * x_i) const
|
||
{
|
||
std::complex<real_t> * x = new std::complex<real_t>[m];
|
||
if (x_r && x_i)
|
||
{
|
||
for (int i = 0; i<m; i++)
|
||
{
|
||
x[i] = std::complex<real_t>(x_r[i], x_i[i]);
|
||
}
|
||
}
|
||
else if (data_r)
|
||
{
|
||
for (int i = 0; i<m; i++)
|
||
{
|
||
x[i] = std::complex<real_t>(x_r[i], 0.);
|
||
}
|
||
}
|
||
else if (data_i)
|
||
{
|
||
for (int i = 0; i<m; i++)
|
||
{
|
||
x[i] = std::complex<real_t>(0., x_i[i]);
|
||
}
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("ComplexFactors::RealToComplex:both real and imag part are null");
|
||
return nullptr;
|
||
}
|
||
return x;
|
||
}
|
||
|
||
void ComplexFactors::ComplexToReal(int m, const std::complex<real_t> * x,
|
||
real_t * x_r, real_t * x_i) const
|
||
{
|
||
for (int i = 0; i<m; i++)
|
||
{
|
||
x_r[i] = x[i].real();
|
||
x_i[i] = x[i].imag();
|
||
}
|
||
}
|
||
|
||
|
||
void ComplexFactors::SetComplexData(int m)
|
||
{
|
||
if (data) { return; }
|
||
MFEM_VERIFY(data_r || data_i, "ComplexFactors data not set");
|
||
data = RealToComplex(m,data_r,data_i);
|
||
}
|
||
|
||
|
||
bool ComplexLUFactors::Factor(int m, real_t TOL)
|
||
{
|
||
SetComplexData(m*m);
|
||
#ifdef MFEM_USE_LAPACK
|
||
int info = 0;
|
||
MFEM_VERIFY(data, "Matrix data not set");
|
||
#ifdef MFEM_USE_SINGLE
|
||
if (m) { cgetrf_(&m, &m, data, &m, ipiv, &info); }
|
||
#elif defined MFEM_USE_DOUBLE
|
||
if (m) { zgetrf_(&m, &m, data, &m, ipiv, &info); }
|
||
#endif
|
||
return info == 0;
|
||
#else
|
||
// compiling without LAPACK
|
||
std::complex<real_t> *data_ptr = this->data;
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
// pivoting
|
||
{
|
||
int piv = i;
|
||
real_t a = std::abs(data_ptr[piv+i*m]);
|
||
for (int j = i+1; j < m; j++)
|
||
{
|
||
const real_t b = std::abs(data_ptr[j+i*m]);
|
||
if (b > a)
|
||
{
|
||
a = b;
|
||
piv = j;
|
||
}
|
||
}
|
||
ipiv[i] = piv;
|
||
if (piv != i)
|
||
{
|
||
// swap rows i and piv in both L and U parts
|
||
for (int j = 0; j < m; j++)
|
||
{
|
||
mfem::Swap<std::complex<real_t>>(data_ptr[i+j*m], data_ptr[piv+j*m]);
|
||
}
|
||
}
|
||
}
|
||
|
||
if (abs(data_ptr[i + i*m]) <= TOL)
|
||
{
|
||
return false; // failed
|
||
}
|
||
|
||
real_t fone = 1.0; // TODO: why?
|
||
const std::complex<real_t> a_ii_inv = fone / data_ptr[i+i*m];
|
||
for (int j = i+1; j < m; j++)
|
||
{
|
||
data_ptr[j+i*m] *= a_ii_inv;
|
||
}
|
||
for (int k = i+1; k < m; k++)
|
||
{
|
||
const std::complex<real_t> a_ik = data_ptr[i+k*m];
|
||
for (int j = i+1; j < m; j++)
|
||
{
|
||
data_ptr[j+k*m] -= a_ik * data_ptr[j+i*m];
|
||
}
|
||
}
|
||
}
|
||
#endif
|
||
|
||
return true; // success
|
||
}
|
||
|
||
std::complex<real_t> ComplexLUFactors::Det(int m) const
|
||
{
|
||
std::complex<real_t> det(1.0,0.);
|
||
for (int i=0; i<m; i++)
|
||
{
|
||
if (ipiv[i] != i-ipiv_base)
|
||
{
|
||
det *= -data[m * i + i];
|
||
}
|
||
else
|
||
{
|
||
det *= data[m * i + i];
|
||
}
|
||
}
|
||
return det;
|
||
}
|
||
|
||
void ComplexLUFactors::Mult(int m, int n, real_t *X_r, real_t * X_i) const
|
||
{
|
||
std::complex<real_t> * x = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
// X <- U X
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
std::complex<real_t> x_i = x[i] * data[i+i*m];
|
||
for (int j = i+1; j < m; j++)
|
||
{
|
||
x_i += x[j] * data[i+j*m];
|
||
}
|
||
x[i] = x_i;
|
||
}
|
||
// X <- L X
|
||
for (int i = m-1; i >= 0; i--)
|
||
{
|
||
std::complex<real_t> x_i = x[i];
|
||
for (int j = 0; j < i; j++)
|
||
{
|
||
x_i += x[j] * data[i+j*m];
|
||
}
|
||
x[i] = x_i;
|
||
}
|
||
// X <- P^{-1} X
|
||
for (int i = m-1; i >= 0; i--)
|
||
{
|
||
mfem::Swap<std::complex<real_t>>(x[i], x[ipiv[i]-ipiv_base]);
|
||
}
|
||
x += m;
|
||
}
|
||
ComplexFactors::ComplexToReal(m*n,x,X_r,X_i);
|
||
delete [] x;
|
||
}
|
||
|
||
void ComplexLUFactors::LSolve(int m, int n, real_t *X_r, real_t * X_i) const
|
||
{
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
// X <- P X
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
mfem::Swap<std::complex<real_t>>(x[i], x[ipiv[i]-ipiv_base]);
|
||
}
|
||
// X <- L^{-1} X
|
||
for (int j = 0; j < m; j++)
|
||
{
|
||
const std::complex<real_t> x_j = x[j];
|
||
for (int i = j+1; i < m; i++)
|
||
{
|
||
x[i] -= data[i+j*m] * x_j;
|
||
}
|
||
}
|
||
x += m;
|
||
}
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexLUFactors::USolve(int m, int n, real_t *X_r, real_t * X_i) const
|
||
{
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
// X <- U^{-1} X
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = m-1; j >= 0; j--)
|
||
{
|
||
const std::complex<real_t> x_j = ( x[j] /= data[j+j*m] );
|
||
for (int i = 0; i < j; i++)
|
||
{
|
||
x[i] -= data[i+j*m] * x_j;
|
||
}
|
||
}
|
||
x += m;
|
||
}
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexLUFactors::Solve(int m, int n, real_t *X_r, real_t * X_i) const
|
||
{
|
||
#ifdef MFEM_USE_LAPACK
|
||
std::complex<real_t> * x = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
char trans = 'N';
|
||
int info = 0;
|
||
#ifdef MFEM_USE_SINGLE
|
||
if (m > 0 && n > 0) { cgetrs_(&trans, &m, &n, data, &m, ipiv, x, &m, &info); }
|
||
#elif defined MFEM_USE_DOUBLE
|
||
if (m > 0 && n > 0) { zgetrs_(&trans, &m, &n, data, &m, ipiv, x, &m, &info); }
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
MFEM_VERIFY(!info, "LAPACK: error in ZGETRS");
|
||
ComplexFactors::ComplexToReal(m*n,x,X_r,X_i);
|
||
delete [] x;
|
||
#else
|
||
// compiling without LAPACK
|
||
LSolve(m, n, X_r, X_i);
|
||
USolve(m, n, X_r, X_i);
|
||
#endif
|
||
}
|
||
|
||
void ComplexLUFactors::RightSolve(int m, int n, real_t *X_r, real_t * X_i) const
|
||
{
|
||
std::complex<real_t> * X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> * x = X;
|
||
#ifdef MFEM_USE_LAPACK
|
||
char n_ch = 'N', side = 'R', u_ch = 'U', l_ch = 'L';
|
||
if (m > 0 && n > 0)
|
||
{
|
||
std::complex<real_t> alpha(1.0,0.0);
|
||
#ifdef MFEM_USE_SINGLE
|
||
ctrsm_(&side,&u_ch,&n_ch,&n_ch,&n,&m,&alpha,data,&m,X,&n);
|
||
ctrsm_(&side,&l_ch,&n_ch,&u_ch,&n,&m,&alpha,data,&m,X,&n);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
ztrsm_(&side,&u_ch,&n_ch,&n_ch,&n,&m,&alpha,data,&m,X,&n);
|
||
ztrsm_(&side,&l_ch,&n_ch,&u_ch,&n,&m,&alpha,data,&m,X,&n);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
}
|
||
#else
|
||
// compiling without LAPACK
|
||
// X <- X U^{-1}
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = 0; j < m; j++)
|
||
{
|
||
const std::complex<real_t> x_j = ( x[j*n] /= data[j+j*m]);
|
||
for (int i = j+1; i < m; i++)
|
||
{
|
||
x[i*n] -= data[j + i*m] * x_j;
|
||
}
|
||
}
|
||
++x;
|
||
}
|
||
|
||
// X <- X L^{-1}
|
||
x = X;
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = m-1; j >= 0; j--)
|
||
{
|
||
const std::complex<real_t> x_j = x[j*n];
|
||
for (int i = 0; i < j; i++)
|
||
{
|
||
x[i*n] -= data[j + i*m] * x_j;
|
||
}
|
||
}
|
||
++x;
|
||
}
|
||
#endif
|
||
// X <- X P
|
||
x = X;
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int i = m-1; i >= 0; --i)
|
||
{
|
||
mfem::Swap<std::complex<real_t>>(x[i*n], x[(ipiv[i]-ipiv_base)*n]);
|
||
}
|
||
++x;
|
||
}
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexLUFactors::GetInverseMatrix(int m, real_t *X_r, real_t *X_i) const
|
||
{
|
||
// A^{-1} = U^{-1} L^{-1} P
|
||
// X <- U^{-1} (set only the upper triangular part of X)
|
||
std::complex<real_t> * X = ComplexFactors::RealToComplex(m*m,X_r,X_i);
|
||
std::complex<real_t> * x = X;
|
||
real_t fone = 1.0; // TODO: why?
|
||
for (int k = 0; k < m; k++)
|
||
{
|
||
const std::complex<real_t> minus_x_k = -( x[k] = fone/data[k+k*m] );
|
||
for (int i = 0; i < k; i++)
|
||
{
|
||
x[i] = data[i+k*m] * minus_x_k;
|
||
}
|
||
for (int j = k-1; j >= 0; j--)
|
||
{
|
||
const std::complex<real_t> x_j = ( x[j] /= data[j+j*m] );
|
||
for (int i = 0; i < j; i++)
|
||
{
|
||
x[i] -= data[i+j*m] * x_j;
|
||
}
|
||
}
|
||
x += m;
|
||
}
|
||
// X <- X L^{-1} (use input only from the upper triangular part of X)
|
||
{
|
||
int k = m-1;
|
||
for (int j = 0; j < k; j++)
|
||
{
|
||
const std::complex<real_t> minus_L_kj = -data[k+j*m];
|
||
for (int i = 0; i <= j; i++)
|
||
{
|
||
X[i+j*m] += X[i+k*m] * minus_L_kj;
|
||
}
|
||
for (int i = j+1; i < m; i++)
|
||
{
|
||
X[i+j*m] = X[i+k*m] * minus_L_kj;
|
||
}
|
||
}
|
||
}
|
||
for (int k = m-2; k >= 0; k--)
|
||
{
|
||
for (int j = 0; j < k; j++)
|
||
{
|
||
const std::complex<real_t> L_kj = data[k+j*m];
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
X[i+j*m] -= X[i+k*m] * L_kj;
|
||
}
|
||
}
|
||
}
|
||
// X <- X P
|
||
for (int k = m-1; k >= 0; k--)
|
||
{
|
||
const int piv_k = ipiv[k]-ipiv_base;
|
||
if (k != piv_k)
|
||
{
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
Swap<std::complex<real_t>>(X[i+k*m], X[i+piv_k*m]);
|
||
}
|
||
}
|
||
}
|
||
ComplexFactors::ComplexToReal(m*m,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
|
||
bool ComplexCholeskyFactors::Factor(int m, real_t TOL)
|
||
{
|
||
SetComplexData(m*m);
|
||
#ifdef MFEM_USE_LAPACK
|
||
int info = 0;
|
||
char uplo = 'L';
|
||
MFEM_VERIFY(data, "Matrix data not set");
|
||
#ifdef MFEM_USE_SINGLE
|
||
if (m) {cpotrf_(&uplo, &m, data, &m, &info);}
|
||
#elif defined MFEM_USE_DOUBLE
|
||
if (m) {zpotrf_(&uplo, &m, data, &m, &info);}
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
return info == 0;
|
||
#else
|
||
// Cholesky–Crout algorithm
|
||
real_t fone = 1.0; // TODO: why?
|
||
for (int j = 0; j<m; j++)
|
||
{
|
||
std::complex<real_t> a(0.,0.);
|
||
for (int k = 0; k<j; k++)
|
||
{
|
||
a+=data[j+k*m]*std::conj(data[j+k*m]);
|
||
}
|
||
|
||
MFEM_VERIFY((data[j+j*m] - a).real() > 0.,
|
||
"CholeskyFactors::Factor: The matrix is not SPD");
|
||
|
||
data[j+j*m] = std::sqrt((data[j+j*m] - a).real());
|
||
|
||
if (data[j + j*m].real() <= TOL) { return false; }
|
||
|
||
for (int i = j+1; i<m; i++)
|
||
{
|
||
a = std::complex<real_t>(0.,0.);
|
||
for (int k = 0; k<j; k++)
|
||
{
|
||
a+= data[i+k*m]*std::conj(data[j+k*m]);
|
||
}
|
||
data[i+j*m] = fone/data[j+m*j]*(data[i+j*m] - a);
|
||
}
|
||
}
|
||
return true; // success
|
||
#endif
|
||
}
|
||
|
||
std::complex<real_t> ComplexCholeskyFactors::Det(int m) const
|
||
{
|
||
std::complex<real_t> det(1.0,0.0);
|
||
for (int i=0; i<m; i++)
|
||
{
|
||
det *= data[i + i*m];
|
||
}
|
||
return det;
|
||
}
|
||
|
||
void ComplexCholeskyFactors::LMult(int m, int n, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
// X <- L X
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = m-1; j >= 0; j--)
|
||
{
|
||
std::complex<real_t> x_j = x[j] * data[j+j*m];
|
||
for (int i = 0; i < j; i++)
|
||
{
|
||
x_j += x[i] * data[j+i*m];
|
||
}
|
||
x[j] = x_j;
|
||
}
|
||
x += m;
|
||
}
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexCholeskyFactors::UMult(int m, int n, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
std::complex<real_t> x_i = x[i] * data[i+i*m];
|
||
for (int j = i+1; j < m; j++)
|
||
{
|
||
x_i += x[j] * std::conj(data[j+i*m]);
|
||
}
|
||
x[i] = x_i;
|
||
}
|
||
x += m;
|
||
}
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexCholeskyFactors::LSolve(int m, int n, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
#ifdef MFEM_USE_LAPACK
|
||
char uplo = 'L';
|
||
char trans = 'N';
|
||
char diag = 'N';
|
||
int info = 0;
|
||
|
||
#ifdef MFEM_USE_SINGLE
|
||
ctrtrs_(&uplo, &trans, &diag, &m, &n, data, &m, x, &m, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
ztrtrs_(&uplo, &trans, &diag, &m, &n, data, &m, x, &m, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
MFEM_VERIFY(!info, "ComplexCholeskyFactors:LSolve:: info");
|
||
#else
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
// X <- L^{-1} X
|
||
for (int j = 0; j < m; j++)
|
||
{
|
||
const std::complex<real_t> x_j = (x[j] /= data[j+j*m]);
|
||
for (int i = j+1; i < m; i++)
|
||
{
|
||
x[i] -= data[i+j*m] * x_j;
|
||
}
|
||
}
|
||
x += m;
|
||
}
|
||
#endif
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexCholeskyFactors::USolve(int m, int n, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
#ifdef MFEM_USE_LAPACK
|
||
|
||
char uplo = 'L';
|
||
char trans = 'C';
|
||
char diag = 'N';
|
||
int info = 0;
|
||
|
||
#ifdef MFEM_USE_SINGLE
|
||
ctrtrs_(&uplo, &trans, &diag, &m, &n, data, &m, x, &m, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
ztrtrs_(&uplo, &trans, &diag, &m, &n, data, &m, x, &m, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
MFEM_VERIFY(!info, "ComplexCholeskyFactors:USolve:: info");
|
||
#else
|
||
// X <- L^{-t} X
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = m-1; j >= 0; j--)
|
||
{
|
||
const std::complex<real_t> x_j = ( x[j] /= data[j+j*m] );
|
||
for (int i = 0; i < j; i++)
|
||
{
|
||
x[i] -= std::conj(data[j+i*m]) * x_j;
|
||
}
|
||
}
|
||
x += m;
|
||
}
|
||
#endif
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexCholeskyFactors::Solve(int m, int n, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
#ifdef MFEM_USE_LAPACK
|
||
char uplo = 'L';
|
||
int info = 0;
|
||
std::complex<real_t> *x = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
#ifdef MFEM_USE_SINGLE
|
||
cpotrs_(&uplo, &m, &n, data, &m, x, &m, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
zpotrs_(&uplo, &m, &n, data, &m, x, &m, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
MFEM_VERIFY(!info, "ComplexCholeskyFactors:Solve:: info");
|
||
ComplexFactors::ComplexToReal(m*n,x,X_r,X_i);
|
||
delete x;
|
||
#else
|
||
LSolve(m, n, X_r,X_i);
|
||
USolve(m, n, X_r,X_i);
|
||
#endif
|
||
}
|
||
|
||
void ComplexCholeskyFactors::RightSolve(int m, int n, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
|
||
std::complex<real_t> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
|
||
std::complex<real_t> *x = X;
|
||
#ifdef MFEM_USE_LAPACK
|
||
char side = 'R';
|
||
char uplo = 'L';
|
||
char transt = 'C';
|
||
char trans = 'N';
|
||
char diag = 'N';
|
||
|
||
std::complex<real_t> alpha(1.0,0.0);
|
||
if (m > 0 && n > 0)
|
||
{
|
||
#ifdef MFEM_USE_SINGLE
|
||
ctrsm_(&side,&uplo,&transt,&diag,&n,&m,&alpha,data,&m,x,&n);
|
||
ctrsm_(&side,&uplo,&trans,&diag,&n,&m,&alpha,data,&m,x,&n);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
ztrsm_(&side,&uplo,&transt,&diag,&n,&m,&alpha,data,&m,x,&n);
|
||
ztrsm_(&side,&uplo,&trans,&diag,&n,&m,&alpha,data,&m,x,&n);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
}
|
||
#else
|
||
// X <- X L^{-H}
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = 0; j < m; j++)
|
||
{
|
||
const std::complex<real_t> x_j = ( x[j*n] /= data[j+j*m]);
|
||
for (int i = j+1; i < m; i++)
|
||
{
|
||
x[i*n] -= std::conj(data[i + j*m]) * x_j;
|
||
}
|
||
}
|
||
++x;
|
||
}
|
||
// X <- X L^{-1}
|
||
x = X;
|
||
for (int k = 0; k < n; k++)
|
||
{
|
||
for (int j = m-1; j >= 0; j--)
|
||
{
|
||
const std::complex<real_t> x_j = (x[j*n] /= data[j+j*m]);
|
||
for (int i = 0; i < j; i++)
|
||
{
|
||
x[i*n] -= data[j + i*m] * x_j;
|
||
}
|
||
}
|
||
++x;
|
||
}
|
||
#endif
|
||
ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
void ComplexCholeskyFactors::GetInverseMatrix(int m, real_t * X_r,
|
||
real_t * X_i) const
|
||
{
|
||
// A^{-1} = L^{-t} L^{-1}
|
||
std::complex<real_t> * X = new std::complex<real_t>[m*m];
|
||
#ifdef MFEM_USE_LAPACK
|
||
// copy the lower triangular part of L to X
|
||
for (int i = 0; i<m; i++)
|
||
{
|
||
for (int j = i; j<m; j++)
|
||
{
|
||
X[j+i*m] = data[j+i*m];
|
||
}
|
||
}
|
||
char uplo = 'L';
|
||
int info = 0;
|
||
#ifdef MFEM_USE_SINGLE
|
||
cpotri_(&uplo, &m, X, &m, &info);
|
||
#elif defined MFEM_USE_DOUBLE
|
||
zpotri_(&uplo, &m, X, &m, &info);
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
#endif
|
||
MFEM_VERIFY(!info, "ComplexCholeskyFactors:GetInverseMatrix:: info");
|
||
// fill in the upper triangular part
|
||
for (int i = 0; i<m; i++)
|
||
{
|
||
for (int j = i+1; j<m; j++)
|
||
{
|
||
X[i+j*m] = std::conj(X[j+i*m]);
|
||
}
|
||
}
|
||
#else
|
||
// L^-t * L^-1 (in place)
|
||
real_t fone = 1.0; // TODO: why?
|
||
for (int k = 0; k<m; k++)
|
||
{
|
||
X[k+k*m] = fone/data[k+k*m];
|
||
for (int i = k+1; i < m; i++)
|
||
{
|
||
std::complex<real_t> s(0.,0.);
|
||
for (int j=k; j<i; j++)
|
||
{
|
||
s -= data[i+j*m] * X[j+k*m]/data[i+i*m];
|
||
}
|
||
X[i+k*m] = s;
|
||
}
|
||
}
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
for (int j = i; j < m; j++)
|
||
{
|
||
std::complex<real_t> s(0.,0.);
|
||
for (int k=j; k<m; k++)
|
||
{
|
||
s += X[k+i*m] * std::conj(X[k+j*m]);
|
||
}
|
||
X[j+i*m] = s;
|
||
X[i+j*m] = std::conj(s);
|
||
}
|
||
}
|
||
#endif
|
||
ComplexFactors::ComplexToReal(m*m,X,X_r,X_i);
|
||
delete [] X;
|
||
}
|
||
|
||
} // mfem namespace
|