portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-14 05:07:05 +01:00
Code Releases Activity

Adding several iMatrix utilities

Browse Source
This commit is contained in:
neo
2015-06-10 14:16:33 +09:00
parent e79a57b423
commit f3dd829459
4 changed files with 64 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

48
lib/tensors/Tensor_Ta.h
View File

@ -63,26 +63,70 @@ namespace Grid {
for(int c1=0;c1<N;c1++){
nrm = 0.0;
for(int c2=0;c2<N;c2++)
nrm = real(innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]));
nrm += real(innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]));
nrm = 1.0/sqrt(nrm);
std::cout << "norm : "<< nrm << "\n";
for(int c2=0;c2<N;c2++)
ret._internal[c1][c2]*= nrm;
for (int b=c1+1; b<N; ++b){
decltype(ret._internal[b][b]*ret._internal[b][b]) pr = 0.0;
for(int c=0; c<N; ++c)
pr += ret._internal[c1][c]*ret._internal[b][c];
pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
std::cout << "pr : "<< pr << "\n";
for(int c=0; c<N; ++c){
ret._internal[b][c] -= pr * ret._internal[c1][c];
}
}
}
// assuming the determinant is ok
return ret;
}
///////////////////////////////////////////////
// Determinant function for scalar, vector, matrix
///////////////////////////////////////////////
inline ComplexF Determinant( const ComplexF &arg){ return arg;}
inline ComplexD Determinant( const ComplexD &arg){ return arg;}
inline RealF Determinant( const RealF &arg){ return arg;}
inline RealD Determinant( const RealD &arg){ return arg;}
template<class vtype> inline auto Determinant(const iScalar<vtype>&r) -> iScalar<decltype(Determinant(r._internal))>
{
iScalar<decltype(Determinant(r._internal))> ret;
ret._internal = Determinant(r._internal);
return ret;
}
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
inline auto Determinant(const iMatrix<vtype,N> &arg)-> iScalar<decltype(Determinant(arg._internal[0][0]))>
{
iMatrix<vtype,N> ret(arg);
iScalar<decltype(Determinant(arg._internal[0][0]))> det = 1.0;
/* Conversion of matrix to upper triangular */
for(int i = 0; i < N; i++){
for(int j = 0; j < N; j++){
if(j>i){
vtype ratio = ret._internal[j][i]/ret._internal[i][i];
for(int k = 0; k < N; k++){
ret._internal[j][k] -= ratio * ret._internal[i][k];
}
}
}
}
for(int i = 0; i < N; i++)
det *= ret._internal[i][i];
return det;
}
///////////////////////////////////////////////
// Exponentiate function for scalar, vector, matrix
///////////////////////////////////////////////