MADNESS  0.10.1
tensor_lapack.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31 
32  $Id$
33 */
34 
35 
36 #ifndef MADNESS_LINALG_TENSOR_LAPACK_H__INCLUDED
37 #define MADNESS_LINALG_TENSOR_LAPACK_H__INCLUDED
38 
39 #include <madness/tensor/tensor.h>
40 #include <madness/fortran_ctypes.h>
41 
42 /*!
43  \file tensor_lapack.h
44  \brief Prototypes for a partial interface from Tensor to LAPACK
45  \ingroup linalg
46 @{
47 */
48 
49 namespace madness {
50 
51  /// Computes singular value decomposition of matrix
52 
53  /// \ingroup linalg
54  template <typename T>
55  void svd(const Tensor<T>& a, Tensor<T>& U,
56  Tensor< typename Tensor<T>::scalar_type >& s, Tensor<T>& VT);
57 
58  /// \ingroup linalg
59  template <typename T>
60  void svd_result(Tensor<T>& a, Tensor<T>& U,
61  Tensor< typename Tensor<T>::scalar_type >& s, Tensor<T>& VT, Tensor<T>& work);
62 
63  /// SVD - MATLAB syntax
64 
65  /// call as
66  /// auto [U,s,VT] = svd(A);
67  /// with A=U*S*VT
68  template <typename T>
69  std::tuple<Tensor<T>, Tensor< typename Tensor<T>::scalar_type >, Tensor<T>>
70  svd(const Tensor<T>& A) {
71  Tensor<T> U,VT;
73  svd(A,U,s,VT);
74  return std::make_tuple(U,s,VT);
75  }
76 
77  /// Solves linear equations
78 
79  /// \ingroup linalg
80  template <typename T>
81  void gesv(const Tensor<T>& a, const Tensor<T>& b, Tensor<T>& x);
82 
83  /// Solves linear equations using least squares
84 
85  /// \ingroup linalg
86  template <typename T>
87  void gelss(const Tensor<T>& a, const Tensor<T>& b, double rcond,
88  Tensor<T>& x, Tensor< typename Tensor<T>::scalar_type >& s,
89  long &rank, Tensor<typename Tensor<T>::scalar_type>& sumsq);
90 
91  /// Solves symmetric or Hermitian eigenvalue problem
92 
93  /// \ingroup linalg
94  template <typename T>
95  void syev(const Tensor<T>& A,
96  Tensor<T>& V, Tensor< typename Tensor<T>::scalar_type >& e);
97 
98  /// Solves symmetric or Hermitian eigenvalue problem - MATLAB syntax
99 
100  /// call as
101  /// auto [eval, evec] = syev(A);
102  template <typename T>
103  std::tuple<Tensor< typename Tensor<T>::scalar_type >, Tensor<T>>
104  syev(const Tensor<T>& A) {
105  Tensor<T> V;
107  syev(A,V,e);
108  return std::make_tuple(e,V);
109  }
110  // START BRYAN ADDITION
111  /// Solves non-symmetric or non-Hermitian eigenvalue problem
112 
113  template <typename T>
114  void geev(const Tensor<T>& A, Tensor<T>& VR, Tensor<std::complex<T>>& e);
115 
116  /// Solves non-symmetric or non-Hermitian generalized eigenvalue problem
117 
118  template <typename T>
119  void ggev(const Tensor<T>& A, Tensor<T>& B, Tensor<T>& VR,
120  Tensor<std::complex<T>>& e);
121  // END BRYAN ADDITIONS
122  /// Solves symmetric or Hermitian eigenvalue problem - MATLAB syntax
123 
124  /// Solves linear equations
125 
126  /// \ingroup linalg
127  template <typename T>
128  void gesv(const Tensor<T>& a, const Tensor<T>& b, Tensor<T>& x);
129 
130 
131  /// Solves symmetric or Hermitian generalized eigenvalue problem
132 
133  /// \ingroup linalg
134  template <typename T>
135  void sygv(const Tensor<T>& A, const Tensor<T>& B, int itype,
136  Tensor<T>& V, Tensor< typename Tensor<T>::scalar_type >& e);
137 
138  class World; // UGH!
139  /// Solves symmetric or Hermitian generalized eigenvalue problem
140 
141 
142  // !!!!!!!!!! sygvp and gesvp are now in the ELEMENTAL inteface
143  // /// \ingroup linalg
144  // template <typename T>
145  // void sygvp(World& world, const Tensor<T>& A, const Tensor<T>& B, int itype,
146  // Tensor<T>& V, Tensor< typename Tensor<T>::scalar_type >& e);
147 
148  // /// Solves linear equations
149 
150  // /// \ingroup linalg
151  // template <typename T>
152  // void gesvp(World& world, const Tensor<T>& a, const Tensor<T>& b, Tensor<T>& x);
153 
154  /// Cholesky factorization
155 
156  /// \ingroup linalg
157  template <typename T>
158  void cholesky(Tensor<T>& A);
159 
160  /// rank-revealing Cholesky factorization
161 
162  /// \ingroup linalg
163  template <typename T>
164  void rr_cholesky(Tensor<T>& A, typename Tensor<T>::scalar_type tol, Tensor<integer>& piv, int& rank);
165 
166  /// \ingroup linalg
167  template <typename T>
168  Tensor<T> inverse(const Tensor<T>& A);
169 
170  /// QR decomposition
171  template<typename T>
172  void qr(Tensor<T>& A, Tensor<T>& R);
173 
174  /// LQ decomposition
175  template<typename T>
176  void lq(Tensor<T>& A, Tensor<T>& L);
177  /// LQ decomposition
178  template<typename T>
179  void lq_result(Tensor<T>& A, Tensor<T>& R, Tensor<T>& tau, Tensor<T>& work,bool do_qr);
180 
181  template <typename T>
182  void geqp3(Tensor<T>& A, Tensor<T>& tau, Tensor<integer>& jpvt);
183 
184  /// orgqr generates an M-by-N complex matrix Q with orthonormal columns
185 
186  /// which is defined as the first N columns of a product of K elementary
187  /// reflectors of order M
188  /// Q = H(1) H(2) . . . H(k)
189  /// as returned by ZGEQRF.
190  template <typename T>
191  void orgqr(Tensor<T>& A, const Tensor<T>& tau);
192 
193 
194  /// Dunno
195 
196 // /// \ingroup linalg
197 // template <typename T>
198 // void triangular_solve(const Tensor<T>& L, Tensor<T>& B,
199 // const char* side, const char* transa);
200 
201  /// Runs the tensor test code, returns true on success
202 
203  /// \ingroup linalg
204  bool test_tensor_lapack();
205 
206  /// World/MRA initialization calls this before going multithreaded due to static data in \c dlamch
207 
208  /// \ingroup linalg
209  void init_tensor_lapack();
210 }
211 
212 #include <madness/tensor/elem.h>
213 
214 #endif // MADNESS_LINALG_TENSOR_LAPACK_H__INCLUDED
Definition: test_ar.cc:118
Definition: test_ar.cc:141
A tensor is a multidimension array.
Definition: tensor.h:317
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition: tensor.h:409
static const double R
Definition: csqrt.cc:46
Correspondence between C++ and Fortran types.
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
void rr_cholesky(Tensor< T > &A, typename Tensor< T >::scalar_type tol, Tensor< integer > &piv, int &rank)
Compute the rank-revealing Cholesky factorization.
Definition: lapack.cc:1203
void 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
Definition: lapack.cc:773
void 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.
Definition: lapack.cc:1052
void lq_result(Tensor< T > &A, Tensor< T > &R, Tensor< T > &tau, Tensor< T > &work, bool do_qr)
compute the LQ decomposition of the matrix A = L Q
Definition: lapack.cc:1320
void ggev(const Tensor< T > &A, Tensor< T > &B, Tensor< T > &VR, Tensor< std::complex< T >> &e)
Generalized real-non-symmetric or complex-non-Hermitian eigenproblem.
Definition: lapack.cc:1115
void qr(Tensor< T > &A, Tensor< T > &R)
compute the QR decomposition of the matrix A
Definition: lapack.cc:1273
void lq(Tensor< T > &A, Tensor< T > &R)
compute the LQ decomposition of the matrix A = L Q
Definition: lapack.cc:1297
void cholesky(Tensor< T > &A)
Compute the Cholesky factorization.
Definition: lapack.cc:1174
void init_tensor_lapack()
World/MRA initialization calls this before going multithreaded due to static data in dlamch.
Definition: lapack.cc:1622
void 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.
Definition: lapack.cc:884
void gesv(const Tensor< T > &a, const Tensor< T > &b, Tensor< T > &x)
Solve Ax = b for general A using the LAPACK *gesv routines.
Definition: lapack.cc:804
void orgqr(Tensor< T > &A, const Tensor< T > &tau)
reconstruct the orthogonal matrix Q (e.g. from QR factorization)
Definition: lapack.cc:1382
void 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.
Definition: lapack.cc:739
void geev(const Tensor< T > &A, Tensor< T > &VR, Tensor< std::complex< T >> &e)
Real non-symmetric or complex non-Hermitian eigenproblem.
Definition: lapack.cc:1008
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition: lapack.cc:832
void geqp3(Tensor< T > &A, Tensor< T > &tau, Tensor< integer > &jpvt)
Compute the QR factorization.
Definition: lapack.cc:1235
bool test_tensor_lapack()
Test the Tensor-LAPACK interface ... currently always returns true!
Definition: lapack.cc:1640
void syev(const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Real-symmetric or complex-Hermitian eigenproblem.
Definition: lapack.cc:969
static const double b
Definition: nonlinschro.cc:119
static const double a
Definition: nonlinschro.cc:118
static const double L
Definition: rk.cc:46
static double V(const coordT &r)
Definition: tdse.cc:288
Defines and implements most of Tensor.
void e()
Definition: test_sig.cc:75