Skip to content

Commit 784f9a3

Browse files
PeizeLinPeizeLin
committed
1. add Lapack_Interface
2. add Tensor_Algorithm::inverse_matrix() 3. add Exx_Stress_Cell_Nearest
1 parent 02f2e1f commit 784f9a3

16 files changed

+807
-11
lines changed

‎include/RI/global/Blas_Interface.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,6 @@
1010
#include <string>
1111
#include <stdexcept>
1212

13-
1413
#ifdef __MKL_RI
1514
#include <mkl_trans.h>
1615
#endif

‎include/RI/global/Lapack-Fortran.h‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
// ===================
2+
// Author: Peize Lin
3+
// date: 2022.12.25
4+
// ===================
5+
6+
#pragma once
7+
8+
#include <complex>
9+
10+
namespace RI
11+
{
12+
13+
extern "C"
14+
{
15+
// potrf computes the Cholesky factorization of a real symmetric positive definite matrix
16+
void spotrf_(const char*const uplo, const int*const n, float*const A, const int*const lda, int*const info);
17+
void dpotrf_(const char*const uplo, const int*const n, double*const A, const int*const lda, int*const info);
18+
void cpotrf_(const char*const uplo, const int*const n, std::complex<float>*const A, const int*const lda, int*const info);
19+
void zpotrf_(const char*const uplo, const int*const n, std::complex<double>*const A, const int*const lda, int*const info);
20+
21+
// potri takes potrf's output to perform matrix inversion
22+
void spotri_(const char*const uplo, const int*const n, float*const A, const int*const lda, int*const info);
23+
void dpotri_(const char*const uplo, const int*const n, double*const A, const int*const lda, int*const info);
24+
void cpotri_(const char*const uplo, const int*const n, std::complex<float>*const A, const int*const lda, int*const info);
25+
void zpotri_(const char*const uplo, const int*const n, std::complex<double>*const A, const int*const lda, int*const info);
26+
27+
// solve the eigenproblem Ax=ex, where A is Symmetric
28+
void ssyev_(const char*const jobz, const char*const uplo,
29+
const int*const n, float*const A, const int*const lda, float*const W,
30+
float*const WORK, const int*const lwork, int*const info);
31+
void dsyev_(const char*const jobz, const char*const uplo,
32+
const int*const n, double*const A, const int*const lda, double*const W,
33+
double*const WORK, const int*const lwork, int*const info);
34+
// solve the eigenproblem Ax=ex, where A is Hermitian
35+
void cheev_(const char*const jobz, const char*const uplo,
36+
const int*const n, std::complex<float>*const A, const int*const lda, float*const W,
37+
std::complex<float>*const WORK, const int*const lwork, float*const RWORK, int*const info);
38+
void zheev_(const char*const jobz, const char*const uplo,
39+
const int*const n, std::complex<double>*const A, const int*const lda, double*const W,
40+
std::complex<double>*const WORK, const int*const lwork, double*const RWORK, int*const info);
41+
}
42+
43+
}

‎include/RI/global/Lapack_Interface-Contiguous.h‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
// ===================
2+
// Author: Peize Lin
3+
// date: 2022.12.25
4+
// ===================
5+
6+
#pragma once
7+
8+
#include "Lapack_Interface.h"
9+
#include "Global_Func-2.h"
10+
11+
namespace RI
12+
{
13+
14+
namespace Lapack_Interface
15+
{
16+
// potrf computes the Cholesky factorization of a real symmetric positive definite matrix
17+
template<typename T>
18+
inline int potrf( const char &uplo, const int &n, T*const A )
19+
{
20+
return potrf(uplo, n, A, n);
21+
}
22+
23+
// potri takes potrf's output to perform matrix inversion
24+
template<typename T>
25+
inline int potri( const char &uplo, const int &n, T*const A )
26+
{
27+
return potri(uplo, n, A, n);
28+
}
29+
30+
// solve the eigenproblem Ax=ex, where A is Hermitian
31+
template<typename T>
32+
inline int heev(const char &jobz, const char &uplo,
33+
const int &n, T*const A, Global_Func::To_Real_t<T>*const W)
34+
{
35+
return heev(jobz, uplo, n, A, n, W);
36+
}
37+
}
38+
39+
}

‎include/RI/global/Lapack_Interface-Tensor.h‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
// ===================
2+
// Author: Peize Lin
3+
// date: 2022.12.25
4+
// ===================
5+
6+
#pragma once
7+
8+
#include "Lapack_Interface-Contiguous.h"
9+
#include "Tensor.h"
10+
#include <cassert>
11+
#include <string>
12+
13+
namespace RI
14+
{
15+
16+
namespace Lapack_Interface
17+
{
18+
// potrf computes the Cholesky factorization of a real symmetric positive definite matrix
19+
template<typename T>
20+
inline int potrf( const char &uplo, Tensor<T> &A )
21+
{
22+
assert(A.shape.size()==2);
23+
assert(A.shape[0]==A.shape[1]);
24+
return potrf(uplo, A.shape[0], A.ptr());
25+
}
26+
27+
// potri takes potrf's output to perform matrix inversion
28+
template<typename T>
29+
inline int potri( const char &uplo, Tensor<T> &A )
30+
{
31+
assert(A.shape.size()==2);
32+
assert(A.shape[0]==A.shape[1]);
33+
return potri(uplo, A.shape[0], A.ptr());
34+
}
35+
36+
// solve the eigenproblem Ax=ex, where A is Hermitian
37+
template<typename T>
38+
inline int heev(const char &jobz, const char &uplo,
39+
Tensor<T> &A, std::vector<Global_Func::To_Real_t<T>> &W)
40+
{
41+
assert(A.shape.size()==2);
42+
assert(A.shape[0]==A.shape[1]);
43+
assert(A.shape[0]==W.size());
44+
return heev(jobz, uplo, A.shape[0], A.ptr(), W.data());
45+
}
46+
}
47+
48+
}

‎include/RI/global/Lapack_Interface.h‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
// ===================
2+
// Author: Peize Lin
3+
// date: 2022.12.25
4+
// ===================
5+
6+
#pragma once
7+
8+
#include "Lapack-Fortran.h"
9+
10+
#include <string>
11+
#include <stdexcept>
12+
13+
14+
#ifdef __MKL_RI
15+
#include <mkl_trans.h>
16+
#endif
17+
18+
#define LAPACK_INFO_CHECK(x) if(const int info=(x)) throw std::runtime_error("info="+std::to_string(info)+".\n"+std::string(__FILE__)+" line "+std::to_string(__LINE__));
19+
20+
namespace RI
21+
{
22+
23+
namespace Lapack_Interface
24+
{
25+
// potrf computes the Cholesky factorization of a real symmetric positive definite matrix
26+
inline int potrf( const char &uplo, const int &n, float*const A, const int &lda )
27+
{
28+
int info;
29+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
30+
spotrf_( &uplo_changed, &n, A, &lda, &info );
31+
return info;
32+
}
33+
inline int potrf( const char &uplo, const int &n, double*const A, const int &lda )
34+
{
35+
int info;
36+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
37+
dpotrf_( &uplo_changed, &n, A, &lda, &info );
38+
return info;
39+
}
40+
inline int potrf( const char &uplo, const int &n, std::complex<float>*const A, const int &lda )
41+
{
42+
int info;
43+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
44+
cpotrf_( &uplo_changed, &n, A, &lda, &info );
45+
return info;
46+
}
47+
inline int potrf( const char &uplo, const int &n, std::complex<double>*const A, const int &lda )
48+
{
49+
int info;
50+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
51+
zpotrf_( &uplo_changed, &n, A, &lda, &info );
52+
return info;
53+
}
54+
55+
// potri takes potrf's output to perform matrix inversion
56+
inline int potri( const char &uplo, const int &n, float*const A, const int &lda )
57+
{
58+
int info;
59+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
60+
spotri_( &uplo_changed, &n, A, &lda, &info);
61+
return info;
62+
}
63+
inline int potri( const char &uplo, const int &n, double*const A, const int &lda )
64+
{
65+
int info;
66+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
67+
dpotri_( &uplo_changed, &n, A, &lda, &info);
68+
return info;
69+
}
70+
inline int potri( const char &uplo, const int &n, std::complex<float>*const A, const int &lda )
71+
{
72+
int info;
73+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
74+
cpotri_( &uplo_changed, &n, A, &lda, &info);
75+
return info;
76+
}
77+
inline int potri( const char &uplo, const int &n, std::complex<double>*const A, const int &lda )
78+
{
79+
int info;
80+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
81+
zpotri_( &uplo_changed, &n, A, &lda, &info);
82+
return info;
83+
}
84+
85+
// solve the eigenproblem Ax=ex, where A is Symmetric
86+
inline int syev(const char &jobz, const char &uplo,
87+
const int &n, float*const A, const int &lda, float*const W,
88+
float*const WORK, const int &lwork)
89+
{
90+
int info;
91+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
92+
ssyev_(&jobz, &uplo_changed, &n, A, &lda, W, WORK, &lwork, &info);
93+
return info;
94+
}
95+
inline int syev(const char &jobz, const char &uplo,
96+
const int &n, double*const A, const int &lda, double*const W,
97+
double*const WORK, const int &lwork)
98+
{
99+
int info;
100+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
101+
dsyev_(&jobz, &uplo_changed, &n, A, &lda, W, WORK, &lwork, &info);
102+
return info;
103+
}
104+
// solve the eigenproblem Ax=ex, where A is Hermitian
105+
inline int heev(const char &jobz, const char &uplo,
106+
const int &n, std::complex<float>*const A, const int &lda, float*const W,
107+
std::complex<float>*const WORK, const int &lwork, float*const RWORK)
108+
{
109+
int info;
110+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
111+
cheev_(&jobz, &uplo_changed, &n, A, &lda, W, WORK, &lwork, RWORK, &info);
112+
return info;
113+
}
114+
inline int heev(const char &jobz, const char &uplo,
115+
const int &n, std::complex<double>*const A, const int &lda, double*const W,
116+
std::complex<double>*const WORK, const int &lwork, double*const RWORK)
117+
{
118+
int info;
119+
const char uplo_changed = Blas_Interface::change_uplo(uplo);
120+
zheev_(&jobz, &uplo_changed, &n, A, &lda, W, WORK, &lwork, RWORK, &info);
121+
return info;
122+
}
123+
124+
// solve the eigenproblem Ax=ex, where A is Hermitian
125+
template<typename T,
126+
typename std::enable_if< std::is_arithmetic<T>::value,int>::type =0>
127+
inline int heev(const char &jobz, const char &uplo,
128+
const int &n, T*const A, const int &lda, T*const W)
129+
{
130+
T work_tmp=100;
131+
constexpr int minus_one = -1;
132+
LAPACK_INFO_CHECK(syev(jobz, uplo, n, A, lda, W, &work_tmp, minus_one)); // get best lwork
133+
134+
const int lwork = work_tmp;
135+
std::vector<T> WORK(std::max(1,lwork));
136+
return syev(jobz, uplo, n, A, lda, W, WORK.data(), lwork);
137+
}
138+
template<typename T,
139+
typename std::enable_if< std::is_arithmetic<T>::value,int>::type =0>
140+
inline int heev(const char &jobz, const char &uplo,
141+
const int &n, std::complex<T>*const A, const int &lda, T*const W)
142+
{
143+
std::vector<T> RWORK(std::max(1,3*n-2));
144+
145+
std::complex<T> work_tmp;
146+
constexpr int minus_one = -1;
147+
LAPACK_INFO_CHECK(heev(jobz, uplo, n, A, lda, W, &work_tmp, minus_one, RWORK.data())); // get best lwork
148+
149+
const int lwork = std::real(work_tmp);
150+
std::vector<std::complex<T>> WORK(std::max(1,lwork));
151+
return heev(jobz, uplo, n, A, lda, W, WORK.data(), lwork, RWORK.data());
152+
}
153+
}
154+
155+
}
156+
157+
#undef LAPACK_INFO_CHECK

‎include/RI/global/Tensor.h‎

Lines changed: 20 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -60,16 +60,16 @@ class Tensor
6060

6161

6262
template<typename T>
63-
Tensor<T> operator+ (const Tensor<T> &t1, const Tensor<T> &t2);
63+
extern Tensor<T> operator+ (const Tensor<T> &t1, const Tensor<T> &t2);
6464
template<typename T>
65-
Tensor<T> operator- (const Tensor<T> &t1, const Tensor<T> &t2);
65+
extern Tensor<T> operator- (const Tensor<T> &t1, const Tensor<T> &t2);
6666

6767
template<typename T>
68-
Tensor<T> operator* (const Tensor<T> &t1, const Tensor<T> &t2);
68+
extern Tensor<T> operator* (const Tensor<T> &t1, const Tensor<T> &t2);
6969
template<typename T>
70-
Tensor<T> operator* (const T &t1, const Tensor<T> &t2);
70+
extern Tensor<T> operator* (const T &t1, const Tensor<T> &t2);
7171
template<typename T>
72-
Tensor<T> operator* (const Tensor<T> &t1, const T &t2);
72+
extern Tensor<T> operator* (const Tensor<T> &t1, const T &t2);
7373

7474

7575
namespace Global_Func
@@ -78,6 +78,21 @@ namespace Global_Func
7878
Tensor<Tout> convert(const Tensor<Tin> &t);
7979
}
8080

81+
82+
template<typename T, std::size_t N0>
83+
extern Tensor<T> to_Tensor(const std::array<T,N0> &a);
84+
template<typename T, std::size_t N0, std::size_t N1>
85+
extern Tensor<T> to_Tensor(const std::array<std::array<T,N1>,N0> &a);
86+
template<typename T, std::size_t N0, std::size_t N1, std::size_t N2>
87+
extern Tensor<T> to_Tensor(const std::array<std::array<std::array<T,N2>,N1>,N0> &a);
88+
89+
template<typename T, std::size_t N0>
90+
extern std::array<T,N0> to_array(const Tensor<T> &t);
91+
template<typename T, std::size_t N0, std::size_t N1>
92+
extern std::array<std::array<T,N1>,N0> to_array(const Tensor<T> &t);
93+
template<typename T, std::size_t N0, std::size_t N1, std::size_t N2>
94+
extern std::array<std::array<std::array<T,N2>,N1>,N0> to_array(const Tensor<T> &t);
95+
8196
}
8297

8398
#include "Blas_Interface-Tensor.h"

0 commit comments

Comments
 (0)