MADNESS  version 0.9
Files | Functions
Linear algebra (interface to LAPACK)
Collaboration diagram for Linear algebra (interface to LAPACK):

Files

file  tensor_lapack.h
 Prototypes for a partial interface from Tensor to LAPACK.
 

Functions

template<typename T >
void madness::svd (const Tensor< T > &a, Tensor< T > &U, Tensor< typename Tensor< T >::scalar_type > &s, Tensor< T > &VT)
 Compute the singluar value decomposition of an n-by-m matrix using *gesvd. More...
 
template<typename T >
void madness::svd_result (Tensor< T > &a, Tensor< T > &U, Tensor< typename Tensor< T >::scalar_type > &s, Tensor< T > &VT, Tensor< T > &work)
 same as svd, but it optimizes away the tensor construction: a = U * diag(s) * VT More...
 
template<typename T >
void madness::gesv (const Tensor< T > &a, const Tensor< T > &b, Tensor< T > &x)
 Solve Ax = b for general A using the LAPACK *gesv routines. More...
 
template<typename T >
void madness::gelss (const Tensor< T > &a, const Tensor< T > &b, double rcond, Tensor< T > &x, Tensor< typename Tensor< T >::scalar_type > &s, long &rank, Tensor< typename Tensor< T >::scalar_type > &sumsq)
 Solve Ax = b for general A using the LAPACK *gelss routines. More...
 
template<typename T >
void madness::syev (const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
 Real-symmetric or complex-Hermitian eigenproblem. More...
 
template<typename T >
void madness::sygv (const Tensor< T > &A, const Tensor< T > &B, int itype, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
 Generalized real-symmetric or complex-Hermitian eigenproblem. More...
 
template<typename T >
void madness::cholesky (Tensor< T > &A)
 Compute the Cholesky factorization. More...
 
bool madness::test_tensor_lapack ()
 Test the Tensor-LAPACK interface ... currently always returns true! More...
 
void madness::init_tensor_lapack ()
 World/MRA initialization calls this before going multithreaded due to static data in dlamch. More...
 

Detailed Description

Function Documentation

template<typename T >
void madness::cholesky ( Tensor< T > &  A)

Compute the Cholesky factorization.

Cholesky factorization.

Compute the Cholesky factorization of the symmetric positive definite matrix A

For memory efficiency A is modified inplace. Its upper triangle will hold the result and the lower trianlge will be zeroed such that input = inner(transpose(output),output).

References dpotrf_(), and TENSOR_ASSERT.

Referenced by madness::test_cholesky().

template<typename T >
void madness::gelss ( const Tensor< T > &  a,
const Tensor< T > &  b,
double  rcond,
Tensor< T > &  x,
Tensor< typename Tensor< T >::scalar_type > &  s,
long &  rank,
Tensor< typename Tensor< T >::scalar_type > &  sumsq 
)

Solve Ax = b for general A using the LAPACK *gelss routines.

Solves linear equations using least squares.

A should be a matrix (float, double, float_complex, double_complex) and b should be either a vector, or a matrix with each vector stored in a column (i.e., b[n,nrhs]).

If the LAPACK routine fails, it throws a TensorException with the LAPACK info as the value. Otherwise, it returns the solution(s). The input A and b are unchanged. There is no need to worry about Python/C/Fortran ordering issues. It will solve Ax=b as written.

This from the LAPACK documentation

RCOND   (input) REAL
RCOND is used to determine the effective  rank  of A.
Singular values S(i) <= RCOND*S(1) are treated
as zero.  If RCOND < 0, machine precision is  used
instead.

RANK    (output) INTEGER
The  effective rank of A, i.e., the number of singular
values which are greater than RCOND*S(1).

Finally, the optional vector sumsq will store the sum-of-squares residual in the case of a rectangular matrix (least squares regression).

References a(), b(), madness::copy(), dgelss_(), m, max, mpfr::min(), TENSOR_ASSERT, and madness::transpose().

Referenced by madness::GMRES(), KAIN(), madness::KAIN(), optimize_coeffs(), and madness::test_gelss().

template<typename T >
void madness::gesv ( const Tensor< T > &  a,
const Tensor< T > &  b,
Tensor< T > &  x 
)

Solve Ax = b for general A using the LAPACK *gesv routines.

Solves linear equations.

A should be a square matrix (float, double, float_complex, double_complex) and b should be either a vector, or a matrix with each vector stored in a column (i.e., b[n,nrhs]).

If the LAPACK routine fails, it throws a TensorException with the LAPACK info as the value. Otherwise, it returns the solution(s). The input A and b are unchanged. There is no need to worry about Python/C/Fortran ordering issues. It will solve Ax=b as written.

References a(), b(), madness::copy(), dgesv_(), m, TENSOR_ASSERT, and madness::transpose().

Referenced by madness::gesvp(), and madness::test_gesv().

void madness::init_tensor_lapack ( )

World/MRA initialization calls this before going multithreaded due to static data in dlamch.

References dlamch_(), and slamch_().

Referenced by madness::startup().

template<typename T >
void madness::svd ( const Tensor< T > &  a,
Tensor< T > &  U,
Tensor< typename Tensor< T >::scalar_type > &  s,
Tensor< T > &  VT 
)

Compute the singluar value decomposition of an n-by-m matrix using *gesvd.

Computes singular value decomposition of matrix.

Returns via arguments U, s, VT where

A = U * diag(s) * VT for A real A = U * diag(s) * VH for A complex

or

UT * A * V = diag(s) for A real UH * A * V = diag(s) for A complex

If A is [m,n] and r=min(m,n) then we have U[m,r], s[r], and VT[r,n]

On failure, throws TensorException with value set to Lapack's info.

References a(), madness::copy(), dgesvd_(), m, max, mpfr::min(), and TENSOR_ASSERT.

Referenced by madness::ConvolutionData1D< Q >::make_approx(), madness::ortho3(), madness::ortho5(), madness::test_svd(), and madness::TensorTrain< T >::truncate().

template<typename T >
void madness::svd_result ( Tensor< T > &  a,
Tensor< T > &  U,
Tensor< typename Tensor< T >::scalar_type > &  s,
Tensor< T > &  VT,
Tensor< T > &  work 
)

same as svd, but it optimizes away the tensor construction: a = U * diag(s) * VT

note that S and VT are swapped in the calling list for c/fortran consistency!

Parameters
[in,out]aa (m,n) matrix to be svd'ed; upon return will hold VT the first min(m,n) rows of VT, stored rowwise
[in,out]Uleft singular vectors, stored columnwise
[in,out]sthe singular values
[in,out]VTnot referenced
[in,out]workwork array; optimial length lwork: lwork = max<integer>(3*min(m,n)+max(m,n),5*min(m,n)-4)*32;

References a(), dgesvd_(), m, and TENSOR_ASSERT.

Referenced by madness::TensorTrain< T >::decompose().

template<typename T >
void madness::syev ( const Tensor< T > &  A,
Tensor< T > &  V,
Tensor< typename Tensor< T >::scalar_type > &  e 
)

Real-symmetric or complex-Hermitian eigenproblem.

Solves symmetric or Hermitian eigenvalue problem.

A is a real symmetric or complex Hermitian matrix. Return V and e where V is a matrix whose columns are the eigenvectors and e is a vector containing the corresponding eigenvalues. If the LAPACK routine fails, it raises a TensorException with value=infor. The input matrix is unchanged. The eigenvalues are sorted into ascending order. s/dsyev are used for real symmetric matrices; c/zheev are used for complex Hermitian.

The reults will satisfy A*V(_,i) = V(_,i)*e(i).

References dsyev_(), max, TENSOR_ASSERT, and madness::transpose().

Referenced by madness::Solver< T, NDIM >::csqrt(), madness::Solver< T, NDIM >::initial_guess(), Molecule::orient(), madness::Molecule::orient(), madness::ortho3(), madness::ortho5(), madness::Solver< T, NDIM >::print_fock_matrix_eigs(), madness::Solver< T, NDIM >::solve(), sqrt(), and madness::test_syev().

template<typename T >
void madness::sygv ( const Tensor< T > &  A,
const Tensor< T > &  B,
int  itype,
Tensor< T > &  V,
Tensor< typename Tensor< T >::scalar_type > &  e 
)

Generalized real-symmetric or complex-Hermitian eigenproblem.

Solves symmetric or Hermitian generalized eigenvalue problem.

This from the LAPACK documentation

S/DSYGV computes all the eigenvalues, and optionally, the eigenvectors
of a real generalized symmetric-definite eigenproblem, of the form
A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.  Here A and B
are assumed to be symmetric and B is also positive definite.

C/ZHEGV computes all the eigenvalues, and optionally, the eigenvectors
of a complex generalized Hermitian-definite eigenproblem, of the form
A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and B
are assumed to be Hermitian and B is also positive definite.

ITYPE   (input) INTEGER
Specifies the problem type to be solved:
= 1:  A*x = (lambda)*B*x
= 2:  A*B*x = (lambda)*x
= 3:  B*A*x = (lambda)*x

References b(), dsygv_(), max, TENSOR_ASSERT, and madness::transpose().

Referenced by SCF::diag_fock_matrix(), madness::Solver< T, NDIM >::do_rhs(), madness::Solver< T, NDIM >::do_rhs_simple(), doit(), madness::Solver< T, NDIM >::initial_guess(), madness::Solver< T, NDIM >::print_fock_matrix_eigs(), madness::Solver< T, NDIM >::print_potential_matrix_eigs(), madness::sygvp(), and madness::test_sygv().

bool madness::test_tensor_lapack ( )

Test the Tensor-LAPACK interface ... currently always returns true!

Dunno.

Runs the tensor test code, returns true on success