MatrixOps.h File Reference

Include dependency graph for MatrixOps.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	lup

Functions
void	LinearLeastSquares (Array< double, 2 > &A, Array< double, 1 > &x, Array< double, 1 > &b)
void	MatVecProd (Array< double, 2 > &A, Array< double, 1 > &x, Array< double, 1 > &Ax)
void	SVdecomp (Array< double, 2 > &A, Array< double, 2 > &U, Array< double, 1 > &S, Array< double, 2 > &V)
void	SVdecomp (Array< std::complex< double >, 2 > &A, Array< std::complex< double >, 2 > &U, Array< double, 1 > &S, Array< std::complex< double >, 2 > &V)
void	SymmEigenPairs (const Array< double, 2 > &A, int NumPairs, Array< double, 1 > &Vals, Array< double, 2 > &Vectors)
void	SymmEigenPairs (const Array< std::complex< double >, 2 > &A, int NumPairs, Array< double, 1 > &Vals, Array< std::complex< double >, 2 > &Vectors)
void	PolarOrthogonalize (Array< std::complex< double >, 2 > &A)
const Array< double, 2 >	operator* (const Array< double, 2 > &A, const Array< double, 2 > &B)
double	InnerProduct (const Array< double, 1 > &A, const Array< double, 1 > &B)
void	OuterProduct (const Array< double, 1 > &A, const Array< double, 1 > &B, Array< double, 2 > &AB)
void	MatMult (const Array< double, 2 > &A, const Array< double, 2 > &B, Array< double, 2 > &C)
void	MatMult (const Array< std::complex< double >, 2 > &A, const Array< std::complex< double >, 2 > &B, Array< std::complex< double >, 2 > &C)
double	Determinant (const Array< double, 2 > &A)
std::complex< double >	Determinant (const Array< std::complex< double >, 2 > &A)
double	DetCofactors (Array< double, 2 > &A, Array< double, 1 > &work)
	This function returns the determinant of A and replaces A with its cofactors. More...
int	DetCofactorsWorksize (int N)
	This function returns the worksize needed by the previous function. More...
std::complex< double >	ComplexDetCofactors (Array< std::complex< double >, 2 > &A, Array< std::complex< double >, 1 > &work)
	Complex versions of the functions above. More...
int	ComplexDetCofactorsWorksize (int N)
double	GJInverse (Array< double, 2 > &A)
double	GJInversePartial (Array< double, 2 > &A)
Array< double, 2 >	Inverse (Array< double, 2 > &A)
void	OutOfPlaceTranspose (Array< double, 2 > &A)
void	OutOfPlaceTranspose (Array< std::complex< double >, 2 > &A)
void	Transpose (Array< double, 2 > &A)
void	Transpose (Array< std::complex< double >, 2 > &A)
void	Print (const Mat3 &A)
void	CubicFormula (double a, double b, double c, double d, double &x1, double &x2, double &x3)
void	TestCubicFormula ()
void	EigVals (const Mat3 &M, double &lambda1, double &lambda2, double &lambda3)
void	Eig (const Mat3 &M, Mat3 &U, Vec3 &Lambda)
double	det (const Mat2 &C)
void	CholeskyBig (Array< double, 2 > &A)
Mat3	Cholesky (Mat3 C)
void	TestCholesky ()

Function Documentation

Cholesky()

Mat3 Cholesky ( Mat3  C )

inline

Definition at line 360 of file MatrixOps.h.

References qmcplusplus::Units::charge::C, and qmcplusplus::sqrt().

Referenced by TestCholesky().

{
  Mat3 L;
  double L00Inv;
  L       = 0.0;
  L(0, 0) = sqrt(C(0, 0));
  L00Inv  = 1.0 / L(0, 0);
  L(1, 0) = C(1, 0) * L00Inv;
  L(2, 0) = C(2, 0) * L00Inv;
  L(1, 1) = sqrt(C(1, 1) - L(1, 0) * L(1, 0));
  L(2, 1) = (C(2, 1) - L(2, 0) * L(1, 0)) / L(1, 1);
  L(2, 2) = sqrt(C(2, 2) - (L(2, 0) * L(2, 0) + L(2, 1) * L(2, 1)));
  return (L);
}

CholeskyBig()

void CholeskyBig

(

Array< double, 2 > &

A

)

ComplexDetCofactors()

std::complex<double> ComplexDetCofactors	(	Array< std::complex< double >, 2 > &	A,
		Array< std::complex< double >, 1 > &	work
	)

Complex versions of the functions above.

ComplexDetCofactorsWorksize()

int ComplexDetCofactorsWorksize

(

int

N

)

CubicFormula()

void CubicFormula ( double  a, double  b, double  c, double  d, double &  x1, double &  x2, double &  x3  )

inline

Definition at line 226 of file MatrixOps.h.

References qmcplusplus::Units::distance::A, qmcplusplus::acos(), B(), qmcplusplus::Units::charge::C, qmcplusplus::cos(), and qmcplusplus::sqrt().

Referenced by EigVals(), and TestCubicFormula().

{
  double A = b / a;
  double B = c / a;
  double C = d / a;
  double Q = (A * A - 3.0 * B) / 9.0;
  double R = (2.0 * A * A * A - 9.0 * A * B + 27.0 * C) / 54.0;
  //cerr << "Q = " << Q << " R = " << R << "\n";
  if ((R * R) < (Q * Q * Q))
  {
    double theta    = acos(R / sqrt(Q * Q * Q));
    double twosqrtQ = 2.0 * sqrt(Q);
    double third    = 1.0 / 3.0;
    double thirdA   = third * A;
    x1              = -twosqrtQ * cos(third * theta) - thirdA;
    x2              = -twosqrtQ * cos(third * (theta + 2.0 * M_PI)) - thirdA;
    x3              = -twosqrtQ * cos(third * (theta - 2.0 * M_PI)) - thirdA;
  }
  else
  {
    std::cerr << "Complex roots detected in CubicFormula.\n";
    exit(1);
  }
}

det()

double det ( const Mat2 &  C )

inline

Definition at line 299 of file MatrixOps.h.

References qmcplusplus::Units::charge::C.

{ return (C(0, 0) * C(1, 1) - C(0, 1) * C(1, 0)); }

DetCofactors()

double DetCofactors	(	Array< double, 2 > &	A,
		Array< double, 1 > &	work
	)

This function returns the determinant of A and replaces A with its cofactors.

DetCofactorsWorksize()

int DetCofactorsWorksize

(

int

N

)

This function returns the worksize needed by the previous function.

Determinant() [1/2]

double Determinant

(

const Array< double, 2 > &

A

)

Determinant() [2/2]

std::complex<double> Determinant

(

const Array< std::complex< double >, 2 > &

A

)

Eig()

void Eig ( const Mat3 &  M, Mat3 &  U, Vec3 &  Lambda  )

inline

Definition at line 282 of file MatrixOps.h.

References EigVals(), norm(), and qmcplusplus::sqrt().

Referenced by TestCholesky().

{
  EigVals(M, Lambda(0), Lambda(1), Lambda(2));
  for (int i = 0; i < 3; i++)
  {
    double L    = Lambda(i);
    U(0, i)     = 1.0;
    U(1, i)     = (M(1, 2) / M(0, 2) * (M(0, 0) - L) - M(1, 0)) / (M(1, 1) - M(1, 2) / M(0, 2) * M(0, 1) - L);
    U(2, i)     = ((L - M(1, 1)) * U(1, i) - M(1, 0)) / M(1, 2);
    double norm = 1.0 / sqrt(U(0, i) * U(0, i) + U(1, i) * U(1, i) + U(2, i) * U(2, i));
    U(0, i) *= norm;
    U(1, i) *= norm;
    U(2, i) *= norm;
  }
}

EigVals()

void EigVals ( const Mat3 &  M, double &  lambda1, double &  lambda2, double &  lambda3  )

inline

Definition at line 270 of file MatrixOps.h.

References CubicFormula().

Referenced by Eig().

{
  double a = -1.0;
  double b = M(0, 0) + M(1, 1) + M(2, 2);
  double c = (-M(0, 0) * M(1, 1) - M(0, 0) * M(2, 2) - M(1, 1) * M(2, 2) + M(0, 1) * M(1, 0) + M(0, 2) * M(2, 0) +
              M(1, 2) * M(2, 1));
  double d = (M(0, 0) * M(1, 1) * M(2, 2) + M(0, 1) * M(1, 2) * M(2, 0) + M(0, 2) * M(1, 0) * M(2, 1) -
              M(0, 1) * M(1, 0) * M(2, 2) - M(0, 2) * M(1, 1) * M(2, 0) - M(0, 0) * M(1, 2) * M(2, 1));
  CubicFormula(a, b, c, d, lambda1, lambda2, lambda3);
}

GJInverse()

double GJInverse

(

Array< double, 2 > &

A

)

Referenced by NonlinearFitClass< M, ModelType >::Fit(), and NonlinearFitClass< M, ModelType >::Solve().

GJInversePartial()

double GJInversePartial

(

Array< double, 2 > &

A

)

InnerProduct()

double InnerProduct	(	const Array< double, 1 > &	A,
		const Array< double, 1 > &	B
	)

Inverse()

Array<double, 2> Inverse

(

Array< double, 2 > &

A

)

LinearLeastSquares()

void LinearLeastSquares	(	Array< double, 2 > &	A,
		Array< double, 1 > &	x,
		Array< double, 1 > &	b
	)

MatMult() [1/2]

void MatMult	(	const Array< double, 2 > &	A,
		const Array< double, 2 > &	B,
		Array< double, 2 > &	C
	)

MatMult() [2/2]

void MatMult	(	const Array< std::complex< double >, 2 > &	A,
		const Array< std::complex< double >, 2 > &	B,
		Array< std::complex< double >, 2 > &	C
	)

MatVecProd()

void MatVecProd	(	Array< double, 2 > &	A,
		Array< double, 1 > &	x,
		Array< double, 1 > &	Ax
	)

operator*()

const Array<double, 2> operator*	(	const Array< double, 2 > &	A,
		const Array< double, 2 > &	B
	)

OuterProduct()

void OuterProduct	(	const Array< double, 1 > &	A,
		const Array< double, 1 > &	B,
		Array< double, 2 > &	AB
	)

OutOfPlaceTranspose() [1/2]

void OutOfPlaceTranspose ( Array< double, 2 > &  A )

inline

Definition at line 166 of file MatrixOps.h.

References qmcplusplus::Units::distance::A, qmcplusplus::Units::distance::m, and qmcplusplus::n.

Referenced by Transpose().

{
  int m = A.rows();
  int n = A.cols();
  Array<double, 2> Atrans(n, m);
  for (int i = 0; i < m; i++)
    for (int j = 0; j < n; j++)
      Atrans(j, i) = A(i, j);
  A.resize(n, m);
  A = Atrans;
}

OutOfPlaceTranspose() [2/2]

void OutOfPlaceTranspose ( Array< std::complex< double >, 2 > &  A )

inline

Definition at line 178 of file MatrixOps.h.

References qmcplusplus::Units::distance::A, qmcplusplus::Units::distance::m, and qmcplusplus::n.

{
  int m = A.rows();
  int n = A.cols();
  Array<std::complex<double>, 2> Atrans(n, m);
  for (int i = 0; i < m; i++)
    for (int j = 0; j < n; j++)
      Atrans(j, i) = A(i, j);
  A.resize(n, m);
  A = Atrans;
}

PolarOrthogonalize()

void PolarOrthogonalize

(

Array< std::complex< double >, 2 > &

A

)

Print()

void Print ( const Mat3 &  A )

inline

Definition at line 219 of file MatrixOps.h.

References qmcplusplus::Units::distance::A.

{
  fprintf(stderr, "%14.6e %14.6e %14.6e\n", A(0, 0), A(0, 1), A(0, 2));
  fprintf(stderr, "%14.6e %14.6e %14.6e\n", A(1, 0), A(1, 1), A(1, 2));
  fprintf(stderr, "%14.6e %14.6e %14.6e\n", A(2, 0), A(2, 1), A(2, 2));
}

SVdecomp() [1/2]

void SVdecomp	(	Array< double, 2 > &	A,
		Array< double, 2 > &	U,
		Array< double, 1 > &	S,
		Array< double, 2 > &	V
	)

SVdecomp() [2/2]

void SVdecomp	(	Array< std::complex< double >, 2 > &	A,
		Array< std::complex< double >, 2 > &	U,
		Array< double, 1 > &	S,
		Array< std::complex< double >, 2 > &	V
	)

SymmEigenPairs() [1/2]

void SymmEigenPairs	(	const Array< double, 2 > &	A,
		int	NumPairs,
		Array< double, 1 > &	Vals,
		Array< double, 2 > &	Vectors
	)

SymmEigenPairs() [2/2]

void SymmEigenPairs	(	const Array< std::complex< double >, 2 > &	A,
		int	NumPairs,
		Array< double, 1 > &	Vals,
		Array< std::complex< double >, 2 > &	Vectors
	)

TestCholesky()

void TestCholesky ( )

inline

Definition at line 376 of file MatrixOps.h.

References qmcplusplus::Units::charge::C, Cholesky(), Eig(), and Transpose().

{
  Mat3 V, D, C, L;
  V(0, 0) = 1.0;
  V(0, 1) = 2.0;
  V(0, 2) = 3.0;
  V(1, 0) = 6.0;
  V(1, 1) = 1.0;
  V(1, 2) = 9.0;
  V(2, 0) = 7.0;
  V(2, 1) = 5.0;
  V(2, 2) = 7.0;

  D(0, 0) = 1.0;
  D(0, 1) = 0.0;
  D(0, 2) = 0.0;
  D(1, 0) = 0.0;
  D(1, 1) = 2.0;
  D(1, 2) = 0.0;
  D(2, 0) = 0.0;
  D(2, 1) = 0.0;
  D(2, 2) = 3.0;

  C = V * D * Transpose(V);
  std::cerr << "C = " << C << "\n";
  std::cerr << "V = " << V << "\n";
  L = Cholesky(C);
  std::cerr << "L = " << L << "\n";
  std::cerr << "L L^T = " << L * Transpose(L) << "\n";

  Mat3 E, U;
  Vec3 Lambda;
  for (int i = 0; i < 10000000; i++)
  {
    Eig(C, U, Lambda);
    //E = Inverse(C);
    //L = Cholesky (E);
  }
}

TestCubicFormula()

void TestCubicFormula ( )

inline

Definition at line 252 of file MatrixOps.h.

References CubicFormula().

{
  double x1 = 2.3;
  double x2 = 1.1;
  double x3 = 9.2;
  double a  = 1.0;
  double b  = -(x1 + x2 + x3);
  double c  = (x1 * x2 + x1 * x3 + x2 * x3);
  double d  = -x1 * x2 * x3;

  double y1, y2, y3;
  CubicFormula(a, b, c, d, y1, y2, y3);
  std::cerr << "y1 = " << y1 << "\n";
  std::cerr << "y2 = " << y2 << "\n";
  std::cerr << "y3 = " << y3 << "\n";
}

Transpose() [1/2]

void Transpose ( Array< double, 2 > &  A )

inline

Definition at line 190 of file MatrixOps.h.

References qmcplusplus::Units::distance::A, qmcplusplus::Units::distance::m, qmcplusplus::n, and OutOfPlaceTranspose().

Referenced by TestCholesky().

{
  int m = A.rows();
  int n = A.cols();
  if (m != n)
    OutOfPlaceTranspose(A);
  else
  {
    for (int i = 0; i < m; i++)
      for (int j = i + 1; j < m; j++)
        std::swap(A(i, j), A(j, i));
  }
}

Transpose() [2/2]

void Transpose ( Array< std::complex< double >, 2 > &  A )

inline

Definition at line 204 of file MatrixOps.h.

References qmcplusplus::Units::distance::A, qmcplusplus::Units::distance::m, qmcplusplus::n, and OutOfPlaceTranspose().

{
  int m = A.rows();
  int n = A.cols();
  if (m != n)
    OutOfPlaceTranspose(A);
  else
  {
    for (int i = 0; i < m; i++)
      for (int j = i + 1; j < m; j++)
        std::swap(A(i, j), A(j, i));
  }
}