ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00034 #ifndef MAT_BISECTION 00035 #define MAT_BISECTION 00036 #include <cmath> 00037 namespace mat { 00045 template<typename Treal> 00046 inline int sign(Treal value) { 00047 if (value > 0) 00048 return 1; 00049 else if (value < 0) 00050 return -1; 00051 else 00052 return 0; 00053 } 00054 00055 00067 template<typename Treal, typename Tfun> 00068 Treal bisection(Tfun const & fun, Treal min, Treal max, Treal const tol) { 00069 int sign_min = sign(fun.eval(min)); 00070 int sign_max = sign(fun.eval(max)); 00071 if (sign_min == sign_max) 00072 throw Failure("bisection(Tfun&, Treal, Treal, Treal): interval " 00073 "incorrect"); 00074 Treal middle = (max + min) / 2; 00075 int sign_middle = sign(fun.eval(middle)); 00076 while (template_blas_fabs(max - min) > tol * 2 && sign_middle != 0) { 00077 if (sign_middle == sign_min) { 00078 min = middle; 00079 sign_min = sign_middle; 00080 } 00081 else { /* (sign_middle == sign_max) */ 00082 max = middle; 00083 sign_max = sign_middle; 00084 } 00085 middle = (max + min) / 2; 00086 sign_middle = sign(fun.eval(middle)); 00087 } 00088 return middle; 00089 } 00090 00091 } /* end namespace mat */ 00092 #endif