1 /* 2 * This source code is part of 3 * 4 * HelFEM 5 * - 6 * Finite element methods for electronic structure calculations on small systems 7 * 8 * Written by Susi Lehtola, 2018- 9 * Copyright (c) 2018- Susi Lehtola 10 * 11 * This program is free software; you can redistribute it and/or 12 * modify it under the terms of the GNU General Public License 13 * as published by the Free Software Foundation; either version 2 14 * of the License, or (at your option) any later version. 15 */ 16 #include "legendretable.h" 17 #include <algorithm> 18 #include "../legendre/Legendre_Wrapper.h" 19 20 namespace helfem { 21 namespace legendretable { operator <(const legendre_table_t & lh,const legendre_table_t & rh)22 bool operator<(const legendre_table_t & lh, const legendre_table_t & rh) { 23 return lh.xi < rh.xi; 24 } 25 LegendreTable()26 LegendreTable::LegendreTable() { 27 Lpad=-1; 28 Lmax=-1; 29 Mmax=-1; 30 } 31 LegendreTable(int Lpad_,int Lmax_,int Mmax_)32 LegendreTable::LegendreTable(int Lpad_, int Lmax_, int Mmax_) : Lpad(Lpad_), Lmax(Lmax_), Mmax(Mmax_) { 33 } 34 ~LegendreTable()35 LegendreTable::~LegendreTable() { 36 } 37 get_index(double xi,bool check) const38 size_t LegendreTable::get_index(double xi, bool check) const { 39 legendre_table_t p; 40 p.xi=xi; 41 42 std::vector<legendre_table_t>::const_iterator low(std::lower_bound(stor.begin(),stor.end(),p)); 43 if(check && low == stor.end()) { 44 std::ostringstream oss; 45 oss << "Could not find xi=" << xi << " on the list!\n"; 46 throw std::logic_error(oss.str()); 47 } 48 49 // Index is 50 size_t idx(low-stor.begin()); 51 if(check && (stor[idx].xi != xi)) { 52 std::ostringstream oss; 53 oss << "Map error: tried to get xi = " << xi << " but got xi = " << stor[idx].xi << "!\n"; 54 throw std::logic_error(oss.str()); 55 } 56 57 return idx; 58 } 59 compute(double xi)60 void LegendreTable::compute(double xi) { 61 #ifdef _OPENMP 62 #pragma omp critical 63 #endif 64 { 65 legendre_table_t entry; 66 67 // Allocate memory 68 entry.xi=xi; 69 entry.Plm.zeros(Lpad+1,Lpad+1); 70 entry.Qlm.zeros(Lpad+1,Lpad+1); 71 72 // Compute 73 if(xi!=1.0) { 74 ::calc_Plm_arr(entry.Plm.memptr(),Lpad,Lpad,xi); 75 ::calc_Qlm_arr(entry.Qlm.memptr(),Lpad,Lpad,xi); 76 } 77 78 // Store only 0 to lmax 79 entry.Plm=entry.Plm.submat(0,0,Lmax,Mmax); 80 entry.Qlm=entry.Qlm.submat(0,0,Lmax,Mmax); 81 82 // Get rid of any non-normal entries 83 for(int L=0;L<=Lmax;L++) 84 for(int M=0;M<=Mmax;M++) { 85 if(!std::isnormal(entry.Plm(L,M))) 86 entry.Plm(L,M)=0.0; 87 if(!std::isnormal(entry.Qlm(L,M))) 88 entry.Qlm(L,M)=0.0; 89 } 90 91 if(!stor.size()) 92 stor.push_back(entry); 93 else 94 // Insert at lower bound 95 stor.insert(stor.begin()+get_index(xi,false),entry); 96 } 97 } 98 get_Plm(int l,int m,double xi) const99 double LegendreTable::get_Plm(int l, int m, double xi) const { 100 #ifndef ARMA_NO_DEBUG 101 if(get_index(xi)>stor.size()) { 102 std::ostringstream oss; 103 oss << "Error in get_Plm(" << l << "," << m << "," << xi << "): index " << get_index(xi) << " greater than array size " << stor.size() << "!\n"; 104 throw std::logic_error(oss.str()); 105 } 106 #endif 107 return stor[get_index(xi)].Plm(l,m); 108 } 109 get_Qlm(int l,int m,double xi) const110 double LegendreTable::get_Qlm(int l, int m, double xi) const { 111 #ifndef ARMA_NO_DEBUG 112 if(get_index(xi)>stor.size()) { 113 std::ostringstream oss; 114 oss << "Error in get_Qlm(" << l << "," << m << "," << xi << "): index " << get_index(xi) << " greater than array size " << stor.size() << "!\n"; 115 throw std::logic_error(oss.str()); 116 } 117 #endif 118 return stor[get_index(xi)].Qlm(l,m); 119 } 120 get_Plm(int l,int m,const arma::vec & xi) const121 arma::vec LegendreTable::get_Plm(int l, int m, const arma::vec & xi) const { 122 arma::vec plm(xi.n_elem); 123 for(size_t i=0;i<xi.n_elem;i++) 124 plm(i)=get_Plm(l,m,xi(i)); 125 return plm; 126 } 127 get_Qlm(int l,int m,const arma::vec & xi) const128 arma::vec LegendreTable::get_Qlm(int l, int m, const arma::vec & xi) const { 129 arma::vec qlm(xi.n_elem); 130 for(size_t i=0;i<xi.n_elem;i++) 131 qlm(i)=get_Qlm(l,m,xi(i)); 132 return qlm; 133 } 134 } 135 } 136