ergo
bisection.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program 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  * This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
36 #ifndef MAT_BISECTION
37 #define MAT_BISECTION
38 #include <cmath>
39 namespace mat {
47  template<typename Treal>
48  inline int sign(Treal value) {
49  if (value > 0)
50  return 1;
51  else if (value < 0)
52  return -1;
53  else
54  return 0;
55  }
56 
57 
69  template<typename Treal, typename Tfun>
70  Treal bisection(Tfun const & fun, Treal min, Treal max, Treal const tol) {
71  int sign_min = sign(fun.eval(min));
72  int sign_max = sign(fun.eval(max));
73  if (sign_min == sign_max)
74  throw Failure("bisection(Tfun&, Treal, Treal, Treal): interval "
75  "incorrect");
76  Treal middle = (max + min) / 2;
77  int sign_middle = sign(fun.eval(middle));
78  while (template_blas_fabs(max - min) > tol * 2 && sign_middle != 0) {
79  if (sign_middle == sign_min) {
80  min = middle;
81  sign_min = sign_middle;
82  }
83  else { /* (sign_middle == sign_max) */
84  max = middle;
85  sign_max = sign_middle;
86  }
87  middle = (max + min) / 2;
88  sign_middle = sign(fun.eval(middle));
89  }
90  return middle;
91  }
92 
93 } /* end namespace mat */
94 #endif
Vec_1
Vector< real, real > Vec_1
Definition: test_LanczosSeveralLargestEig.cc:63
mat::MatrixSymmetric::randomZeroStructure
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixSymmetric.h:588
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
mat::maxdiff
static Treal maxdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:44
mat::FileWritable::getStatsCountRead
static std::string getStatsCountRead()
Definition: FileWritable.cc:364
MatrixTriangular.h
testAccumulation
static void testAccumulation(const Tmatrix &syFock, int size, Treal *fockfull)
Definition: API_test.cc:88
Vec_3
Vector< real, Vec_2 > Vec_3
Definition: test_LanczosSeveralLargestEig.cc:65
mat::symv
static void symv(const char *uplo, const int *n, const T *alpha, const T *A, const int *lda, const T *x, const int *incx, const T *beta, T *y, const int *incy)
Definition: mat_gblas.h:400
Vector.h
mat::frobdiff
static Treal frobdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:69
main
int main(int argc, char *argv[])
Definition: API_test.cc:2413
mat::zero
@ zero
Definition: matInclude.h:138
mat::MatrixSymmetric::gershgorin
void gershgorin(Treal &lmin, Treal &lmax) const
Definition: MatrixSymmetric.h:476
get_wall_seconds
static double get_wall_seconds()
Definition: bench_gemm_only.cc:49
mat::trmm
static void trmm(const char *side, const char *uplo, const char *transa, const char *diag, const int *m, const int *n, const T *alpha, const T *A, const int *lda, T *B, const int *ldb)
Definition: mat_gblas.h:281
Mat_1
Matrix< real, real > Mat_1
Definition: test_LanczosSeveralLargestEig.cc:60
normalVector
VectorGeneral< real, vect > normalVector
Definition: test_LanczosSeveralLargestEig.cc:72
mat::MatrixTriangular::get_values
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Get values given by row and column index lists.
Definition: MatrixTriangular.h:115
mat::MatrixSymmetric::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Assign from sparse matrix given by three vectors.
Definition: MatrixSymmetric.h:176
mat::MatrixGeneral::frob_thresh
void frob_thresh(Treal threshold)
Definition: MatrixGeneral.h:299
mat::MatrixSymmetric::get_all_values
void get_all_values(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &values) const
Get all values and corresponding row and column index lists, in matrix.
Definition: MatrixSymmetric.h:303
Mat_2
Matrix< real, Mat_1 > Mat_2
Definition: test_LanczosSeveralLargestEig.cc:61
mat::Params::setNProcs
static void setNProcs(unsigned int const nP)
Definition: matInclude.h:112
mat::MatrixTriangular::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Assign from sparse matrix given by three vectors.
Definition: MatrixTriangular.h:86
Sum
Definition: API_test.cc:71
template_blas_exp
Treal template_blas_exp(Treal x)
symmMatrix
MatrixSymmetric< real, matri > symmMatrix
Definition: test_LanczosSeveralLargestEig.cc:69
main
int main(int argc, char *argv[])
Definition: bench.cc:134
mat::Params::getMatrixParallelLevel
static unsigned int getMatrixParallelLevel()
Definition: matInclude.h:120
frobdiff
static Treal frobdiff(const Treal *f1, const Treal *f2, int size)
Definition: API_test.cc:2464
Matrix.h
mat::VectorGeneral::eucl
Treal eucl() const
Definition: VectorGeneral.h:178
mat::FileWritable::getStatsTimeCopyAndAssign
static std::string getStatsTimeCopyAndAssign()
Definition: FileWritable.cc:358
mat::frobNorm
@ frobNorm
Definition: matInclude.h:139
Lanczos.h
Vec_2
Vector< real, Vec_1 > Vec_2
Definition: test_LanczosSeveralLargestEig.cc:64
expRule
Definition: API_test.cc:78
mat::gemv
static void gemv(const char *ta, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *x, const int *incx, const T *beta, T *y, const int *incy)
Definition: mat_gblas.h:391
mat::MatrixTriangular::eucl_thresh_congr_trans_measure
Treal eucl_thresh_congr_trans_measure(Treal const threshold, MatrixSymmetric< Treal, Tmatrix > &trA)
Definition: MatrixTriangular.h:300
real
ergo_real real
Definition: test.cc:46
mat_gblas.h
rows
mat::SizesAndBlocks rows
Definition: test.cc:51
mat::packedtofull
static void packedtofull(const Treal *packed, Treal *full, const int size)
Definition: mat_gblas.h:1148
mat::MatrixSymmetric::frob
Treal frob() const
Definition: MatrixSymmetric.h:360
maxdiff_tri
static Treal maxdiff_tri(const Treal *f1, const Treal *f2, int size)
Definition: API_test.cc:2451
template_blas_fabs
Treal template_blas_fabs(Treal x)
mat::Params::setMatrixParallelLevel
static void setMatrixParallelLevel(unsigned int const mPL)
Definition: matInclude.h:129
mat::MatrixGeneral::eucl
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixGeneral.h:476
Sum::accumulate
Treal accumulate(const Treal &a, int row, int col)
Definition: API_test.cc:73
Mat_3
Matrix< real, Mat_2 > Mat_3
Definition: test_LanczosSeveralLargestEig.cc:62
Mat_4
mat::Matrix< ergo_real, Mat_3 > Mat_4
Definition: matrix_typedefs.h:55
B
#define B
mat::MatrixTriangular::inch
void inch(const MatrixGeneral< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:150
mat::sign
int sign(Treal value)
Sign function returns the sign of the input.
Definition: bisection.h:48
mat::pptrf
static void pptrf(const char *uplo, const int *n, T *ap, int *info)
Definition: mat_gblas.h:263
normalMatrix
MatrixGeneral< real, matri > normalMatrix
Definition: test_LanczosSeveralLargestEig.cc:71
triangMatrix
MatrixTriangular< real, matri > triangMatrix
Definition: test_LanczosSeveralLargestEig.cc:70
mat::accumulate
Treal accumulate(MatrixSymmetric< Treal, Tmatrix > &A, Top &op)
Performs operation specified in 'op' on all nonzero matrix elements and sums up the result and return...
Definition: MatrixSymmetric.h:1251
mat::MatrixTriangular::frob
Treal frob() const
Definition: MatrixTriangular.h:169
mat::MatrixSymmetric::setElementsByRule
void setElementsByRule(TRule &rule)
Uses rule depending on the row and column indexes to set matrix elements The Trule class should have ...
Definition: MatrixSymmetric.h:597
mainFun
int mainFun(int argc, char *argv[])
Definition: bench.cc:56
mat::MatrixSymmetric
Symmetric matrix.
Definition: MatrixBase.h:51
mat::MatrixTriangular
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
mat::VectorGeneral::resetSizesAndBlocks
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
matri
Mat_3 matri
Definition: test_LanczosSeveralLargestEig.cc:67
mat::MatrixSymmetric::eucl
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:673
mat::arn::LanczosLargestMagnitudeEig
Definition: LanczosLargestMagnitudeEig.h:47
mat::MatrixSymmetric::eucl_thresh
Treal eucl_thresh(Treal const threshold, MatrixTriangular< Treal, Tmatrix > const *const Zptr=NULL)
Definition: MatrixSymmetric.h:857
mat::Failure
Definition: Failure.h:57
mat::MatrixTriangular::eucl
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixTriangular.h:257
mat::Params::getNProcs
static unsigned int getNProcs()
Definition: matInclude.h:103
expRule::set
Treal set(int const row, int const col)
Definition: API_test.cc:80
mat::MatrixSymmetric::mixed_norm_thresh
Treal mixed_norm_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:913
dotIsBroken
bool dotIsBroken()
Definition: API_test.cc:105
mat::tripackedtofull
static void tripackedtofull(const Treal *packed, Treal *full, const int size)
Definition: mat_gblas.h:1170
mat::FileWritable::getStatsTimeRead
static std::string getStatsTimeRead()
Definition: FileWritable.cc:355
realIsSingle< float >
bool realIsSingle< float >()
Definition: API_test.cc:57
mat::Vector
Vector class.
Definition: Matrix.h:78
mat::gemm
static void gemm(const char *ta, const char *tb, const int *n, const int *k, const int *l, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: mat_gblas.h:232
mat::MatrixGeneral::get_values
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Get values given by row and column index lists.
Definition: MatrixGeneral.h:173
mat::MatrixGeneral::setElementsByRule
void setElementsByRule(TRule &rule)
Definition: MatrixGeneral.h:457
mat::MatrixTriangular::eucl_thresh
Treal eucl_thresh(Treal const threshold)
Definition: MatrixTriangular.h:291
mat::tptri
static void tptri(const char *uplo, const char *diag, const int *n, T *ap, int *info)
Definition: mat_gblas.h:275
mat::VectorGeneral
Definition: MatrixBase.h:55
mat::MatrixGeneral
Normal matrix.
Definition: MatrixBase.h:49
mat::MatrixSymmetric::random
void random()
Definition: MatrixSymmetric.h:584
mat::MatrixGeneral::frob
Treal frob() const
Definition: MatrixGeneral.h:284
mat::transpose
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
cols
mat::SizesAndBlocks cols
Definition: test.cc:52
mat
Definition: allocate.cc:39
A
#define A
mat::FileWritable::getStatsCountCopyAndAssign
static std::string getStatsCountCopyAndAssign()
Definition: FileWritable.cc:367
mat::Failure::what
virtual const char * what() const
Definition: Failure.h:67
mat::Matrix
Matrix class and heart of the matrix library.
Definition: Matrix.h:115
mat::FileWritable::resetStats
static void resetStats()
Definition: FileWritable.cc:308
max
#define max(a, b)
Definition: integrator.cc:87
mat::symm
static void symm(const char *side, const char *uplo, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: mat_gblas.h:342
mat::VectorGeneral::rand
void rand()
Definition: VectorGeneral.h:108
min
int min(int a, int b)
Definition: lin_trans.cc:66
mat::MatrixGeneral::randomZeroStructure
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixGeneral.h:452
mat::VectorGeneral::fullvector
void fullvector(std::vector< Treal > &fullVector) const
Definition: VectorGeneral.h:88
mat::trmv
static void trmv(const char *uplo, const char *trans, const char *diag, const int *n, const T *A, const int *lda, T *x, const int *incx)
Definition: mat_gblas.h:409
mat::MatrixSymmetric::get_values
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Get values given by row and column index lists.
Definition: MatrixSymmetric.h:273
mat::MatrixGeneral::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Assign from sparse matrix given by three arrays.
Definition: MatrixGeneral.h:122
vect
Vec_3 vect
Definition: test_LanczosSeveralLargestEig.cc:68
mat::MatrixTriangular::setElementsByRule
void setElementsByRule(TRule &rule)
Uses rule depending on the row and column indexes to set matrix elements The Trule class should have ...
Definition: MatrixTriangular.h:205
mat::FileWritable::getStatsCountWrite
static std::string getStatsCountWrite()
Definition: FileWritable.cc:361
mat::MatrixTriangular::get_all_values
void get_all_values(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &values) const
Get all values and corresponding row and column index lists, in matrix.
Definition: MatrixTriangular.h:128
mat::SizesAndBlocks
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
mat::syrk
static void syrk(const char *uplo, const char *trans, const int *n, const int *k, const T *alpha, const T *A, const int *lda, const T *beta, T *C, const int *ldc)
Definition: mat_gblas.h:334
mainFun
int mainFun(int argc, char *argv[])
Definition: API_test.cc:122
SizesAndBlocks.h
Class used to keep track of the block sizes used at different levels in the hierarchical matrix data ...
mat::bisection
Treal bisection(Tfun const &fun, Treal min, Treal max, Treal const tol)
Bisection algorithm for root finding.
Definition: bisection.h:70
mainFun
int mainFun(int size, int parallel_level, int block_size, int block_factor_1, int block_factor_2, int block_factor_3)
Definition: bench_gemm_only.cc:60
mat::MatrixSymmetric::assignFromFull
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixSymmetric.h:110
mat::maxdiff_tri
static Treal maxdiff_tri(const Treal *f1, const Treal *f2, int size)
Definition: general.h:56
mat::DiffMatrix
Definition: mat_utils.h:43
mat::FileWritable::getStatsTimeWrite
static std::string getStatsTimeWrite()
Definition: FileWritable.cc:352
MatrixSymmetric.h
realIsSingle
bool realIsSingle()
Definition: API_test.cc:55
mat::MatrixTriangular::random
void random()
Definition: MatrixTriangular.h:196
mat::MatrixSymmetric::fullMatrix
void fullMatrix(std::vector< Treal > &fullMat) const
Save matrix as full matrix.
Definition: MatrixSymmetric.h:138
MatrixGeneral.h
VectorGeneral.h
mat::MatrixGeneral::eucl_thresh
Treal eucl_thresh(Treal const threshold)
Definition: MatrixGeneral.h:525
mat::dot
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition: mat_gblas.h:425
maxdiff
static Treal maxdiff(std::vector< Treal > const &f1, std::vector< Treal > const &f2)
Definition: API_test.cc:2436
main
int main(int argc, char *argv[])
Definition: bench_gemm_only.cc:110