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