1 /*------------------------------------------------------------------- 2 Copyright 2011 Ravishankar Sundararaman 3 4 This file is part of JDFTx. 5 6 JDFTx 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 3 of the License, or 9 (at your option) any later version. 10 11 JDFTx 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 JDFTx. If not, see <http://www.gnu.org/licenses/>. 18 -------------------------------------------------------------------*/ 19 20 #include <core/Spline.h> 21 22 namespace QuinticSpline 23 { getCoeff(const std::vector<double> & samples,bool oddExtension)24 std::vector<double> getCoeff(const std::vector<double>& samples, bool oddExtension) 25 { //Setup a penta-diagonal system of equations for quintic-Blip coefficeints: 26 std::vector<double> x = samples; //destructible copy 27 int N = x.size(); 28 std::vector<double> tmpBand[5]; std::vector<double>* band = tmpBand+2; 29 for(int m=-2; m<=2; m++) band[m].resize(N-abs(m)); 30 const double off1 = 13./33; 31 const double off2 = 1./66; 32 for(int i=0; i<N; i++) band[0][i] = 1.; 33 for(int i=0; i<N-1; i++) band[1][i] = band[-1][i] = off1; 34 for(int i=0; i<N-2; i++) band[2][i] = band[-2][i] = off2; 35 //Mirror BC at rho=0: 36 if(oddExtension) x[0] = 0.; 37 double sign = oddExtension ? -1. : 1.; 38 band[0][1] += sign * off2; 39 band[1][0] += sign * off1; 40 band[2][0] += sign * off2; 41 //Natural BC at rho=rhoMax: (third and fourth derivatives vanish) 42 const double extrap[2][3] = {{+1,-3,+3}, {+3,-8,+6}}; 43 band[-1][N-3] += off2 * extrap[0][0]; 44 band[ 0][N-2] += off2 * extrap[0][1]; 45 band[ 1][N-2] += off2 * extrap[0][2]; 46 band[-2][N-3] += off1 * extrap[0][0] + off2 * extrap[1][0]; 47 band[-1][N-2] += off1 * extrap[0][1] + off2 * extrap[1][1]; 48 band[ 0][N-1] += off1 * extrap[0][2] + off2 * extrap[1][2]; 49 //Pentadiagonal solve: Forward sweep: 50 for(int i=0; i<N; i++) 51 for(int m=1; m<=2; m++) 52 if(i+m<N) 53 { double rot = band[-m][i]/band[0][i]; //gauss elimination factor 54 for(int n=1; n<m; n++) 55 band[n-m][i+n] -= rot*band[n][i]; 56 for(int n=m; n<=2; n++) 57 if(i+n<N) 58 band[n-m][i+m] -= rot*band[n][i]; 59 x[i+m] -= rot*x[i]; 60 } 61 //Pentadiagonal solve: Reverse sweep: 62 for(int i=N-1; i>=0; i--) 63 { double invDiag = 1.0/band[0][i]; 64 for(int m=1; m<=2; m++) 65 if(i+m<N) x[i] -= band[m][i]*x[i+m]; 66 x[i] *= invDiag; 67 } 68 //Mirror two coefficients at beginning: 69 std::vector<double> coeff; 70 coeff.push_back(sign * x[2]); 71 coeff.push_back(sign * x[1]); 72 coeff.insert(coeff.end(), x.begin(), x.end()); 73 //Extrapolate two coefficients for natural BC at rhoMax: 74 coeff.push_back(extrap[0][0]*x[N-3] + extrap[0][1]*x[N-2] + extrap[0][2]*x[N-1]); 75 coeff.push_back(extrap[1][0]*x[N-3] + extrap[1][1]*x[N-2] + extrap[1][2]*x[N-1]); 76 return coeff; 77 } 78 } 79