ERKALE
ERKALE - DFT from Hel
 All Classes Functions Variables Friends Pages
linalg.h
1 /*
2  * This source code is part of
3  *
4  * E R K A L E
5  * -
6  * DFT from Hel
7  *
8  * Written by Susi Lehtola, 2010-2011
9  * Copyright (c) 2010-2011, 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 
17 
18 
19 #include "global.h"
20 class Settings;
21 
22 #ifndef ERKALE_LINALG
23 #define ERKALE_LINALG
24 
25 #include <armadillo>
26 
28 template<typename T> struct eigenvector {
30  double E;
32  arma::Col<T> c;
33 };
34 
36 template<typename T> inline bool operator<(const struct eigenvector<T> & lhs, const struct eigenvector<T> & rhs) {
37  return lhs.E < rhs.E;
38 }
39 
41 template<typename T> inline void sort_eigvec_wrk(arma::vec & eigval, arma::Mat<T> & eigvec) {
42  if(eigval.n_elem != eigvec.n_cols) {
43  ERROR_INFO();
44  throw std::runtime_error("Eigenvalue vector does not correspond to eigenvector matrix!\n");
45  }
46 
47  // Helper
48  std::vector< struct eigenvector<T> > orbs(eigval.n_elem);
49  for(size_t io=0;io<eigval.n_elem;io++) {
50  orbs[io].E=eigval(io);
51  orbs[io].c=eigvec.col(io);
52  }
53  std::stable_sort(orbs.begin(),orbs.end());
54  for(size_t io=0;io<eigval.n_elem;io++) {
55  eigval(io)=orbs[io].E;
56  eigvec.col(io)=orbs[io].c;
57  }
58 }
59 
61 void sort_eigvec(arma::vec & eigval, arma::mat & eigvec);
63 void sort_eigvec(arma::vec & eigval, arma::cx_mat & eigvec);
64 
68 template<typename T> inline void eig_sym_ordered_wrk(arma::vec & eigval, arma::Mat<T> & eigvec, const arma::Mat<T> & X) {
69  /* Solve eigenvalues of symmetric matrix with guarantee
70  of ordering of eigenvalues from smallest to biggest */
71 
72  // Solve eigenvalues and eigenvectors
73  bool ok=arma::eig_sym(eigval,eigvec,X);
74  if(!ok) {
75  ERROR_INFO();
76  printf("Unable to diagonalize matrix!\n");
77  X.print("X");
78  throw std::runtime_error("Error in eig_sym.\n");
79  }
80 
81  // Sort vectors
82  sort_eigvec_wrk<T>(eigval,eigvec);
83 }
84 
85 
89 void eig_sym_ordered(arma::vec & eigval, arma::mat & eigvec, const arma::mat & X);
93 void eig_sym_ordered(arma::vec & eigval, arma::cx_mat & eigvec, const arma::cx_mat & X);
94 
95 /* Orthogonalization routines */
96 
98 arma::mat CholeskyOrth(const arma::mat & S);
100 arma::mat SymmetricOrth(const arma::mat & S);
102 arma::mat SymmetricOrth(const arma::mat & Svec, const arma::vec & Sval);
104 arma::mat CanonicalOrth(const arma::mat & S, double cutoff=LINTHRES);
106 arma::mat CanonicalOrth(const arma::mat & Svec, const arma::vec & Sval, double cutoff);
107 
109 arma::mat BasOrth(const arma::mat & S, bool verbose);
111 arma::mat BasOrth(const arma::mat & S, const Settings & set);
112 
114 void S_half_invhalf(const arma::mat & S, arma::mat & Shalf, arma::mat & Sinvhalf, bool canonical=false, double cutoff=LINTHRES);
115 
117 arma::vec MatToVec(const arma::mat & v);
119 arma::mat VecToMat(const arma::vec & v, size_t nrows, size_t ncols);
120 
122 arma::vec slicevec(const arma::cube & c, size_t i, size_t j);
123 
125 arma::mat cosmat(const arma::mat & M);
127 arma::mat sinmat(const arma::mat & M);
129 arma::mat sincmat(const arma::mat & M);
131 arma::mat sqrtmat(const arma::mat & M);
133 arma::mat expmat(const arma::mat & M);
135 arma::cx_mat expmat(const arma::cx_mat & M);
136 
138 arma::mat orthogonalize(const arma::mat & M);
140 arma::cx_mat unitarize(const arma::cx_mat & M);
141 
143 arma::mat orthonormalize(const arma::mat & S, const arma::mat & C);
145 arma::cx_mat orthonormalize(const arma::mat & S, const arma::cx_mat & C);
146 
150 void form_NOs(const arma::mat & P, const arma::mat & S, arma::mat & AO_to_NO, arma::mat & NO_to_AO, arma::vec & occs);
151 
155 void form_NOs(const arma::mat & P, const arma::mat & S, arma::mat & AO_to_NO, arma::vec & occs);
156 
165 arma::mat pivoted_cholesky(const arma::mat & M, double thr, arma::uvec & pivot);
167 arma::mat pivoted_cholesky(const arma::mat & M, double thr);
169 arma::mat incomplete_cholesky(const arma::mat & M, size_t n);
170 
172 arma::mat B_transform(arma::mat B, const arma::mat & Cl, const arma::mat & Cr);
173 
175 void check_lapack_thread();
176 
177 #endif
double E
Energy.
Definition: linalg.h:30
Helper for eigenvector sorts.
Definition: linalg.h:28
Settings used for a calculations.
Definition: settings.h:74
arma::Col< T > c
Eigenvector.
Definition: linalg.h:32