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
21 #ifndef ERKALE_LINALG
22 #define ERKALE_LINALG
23
24 #include <armadillo>
25
26 /// Helper for eigenvector sorts
27 template<typename T> struct eigenvector {
28 /// Energy
29 double E;
30 /// Eigenvector
31 arma::Col<T> c;
32 };
33
34 /// Helper for sorts
35 template<typename T> inline bool operator<(const struct eigenvector<T> & lhs, const struct eigenvector<T> & rhs) {
36 return lhs.E < rhs.E;
37 }
38
39 /// Sort vectors in order of increasing eigenvalue
sort_eigvec_wrk(arma::vec & eigval,arma::Mat<T> & eigvec)40 template<typename T> inline void sort_eigvec_wrk(arma::vec & eigval, arma::Mat<T> & eigvec) {
41 if(eigval.n_elem != eigvec.n_cols) {
42 ERROR_INFO();
43 throw std::runtime_error("Eigenvalue vector does not correspond to eigenvector matrix!\n");
44 }
45
46 // Helper
47 std::vector< struct eigenvector<T> > orbs(eigval.n_elem);
48 for(size_t io=0;io<eigval.n_elem;io++) {
49 orbs[io].E=eigval(io);
50 orbs[io].c=eigvec.col(io);
51 }
52 std::stable_sort(orbs.begin(),orbs.end());
53 for(size_t io=0;io<eigval.n_elem;io++) {
54 eigval(io)=orbs[io].E;
55 eigvec.col(io)=orbs[io].c;
56 }
57 }
58
59 /// Wrapper
60 void sort_eigvec(arma::vec & eigval, arma::mat & eigvec);
61 /// Sort vectors in order of increasing eigenvalue
62 void sort_eigvec(arma::vec & eigval, arma::cx_mat & eigvec);
63
64 /**
65 * Solve eigenvalues of symmetric matrix with guarantee
66 * of eigenvalue ordering from smallest to biggest */
eig_sym_ordered_wrk(arma::vec & eigval,arma::Mat<T> & eigvec,const arma::Mat<T> & X)67 template<typename T> inline void eig_sym_ordered_wrk(arma::vec & eigval, arma::Mat<T> & eigvec, const arma::Mat<T> & X) {
68 /* Solve eigenvalues of symmetric matrix with guarantee
69 of ordering of eigenvalues from smallest to biggest */
70
71 // Solve eigenvalues and eigenvectors
72 bool ok=arma::eig_sym(eigval,eigvec,X);
73 if(!ok) {
74 //ERROR_INFO();
75 //printf("Unable to diagonalize matrix!\n");
76 //X.print("X");
77 throw std::runtime_error("Error in eig_sym.\n");
78 }
79
80 // Sort vectors
81 sort_eigvec_wrk<T>(eigval,eigvec);
82 }
83
84
85 /**
86 * Solve eigenvalues of symmetric matrix with guarantee
87 * of eigenvalue ordering from smallest to biggest */
88 void eig_sym_ordered(arma::vec & eigval, arma::mat & eigvec, const arma::mat & X);
89 /**
90 * Solve eigenvalues of Hermitian matrix with guarantee
91 * of eigenvalue ordering from smallest to biggest */
92 void eig_sym_ordered(arma::vec & eigval, arma::cx_mat & eigvec, const arma::cx_mat & X);
93
94 /* Orthogonalization routines */
95
96 /// Cholesky orthogonalization of basis set
97 arma::mat CholeskyOrth(const arma::mat & S);
98 /// Symmetric orthogonalization of basis set
99 arma::mat SymmetricOrth(const arma::mat & S);
100 /// Same, but use computed decomoposition
101 arma::mat SymmetricOrth(const arma::mat & Svec, const arma::vec & Sval);
102 /// Canonical orthogonalization of basis set
103 arma::mat CanonicalOrth(const arma::mat & S, double cutoff);
104 /// Same, but use computed decomposition
105 arma::mat CanonicalOrth(const arma::mat & Svec, const arma::vec & Sval, double cutoff);
106
107 /// Automatic orthonormalization.
108 arma::mat BasOrth(const arma::mat & S, bool verbose);
109 /// Orthogonalize basis
110 arma::mat BasOrth(const arma::mat & S);
111
112 /// Form matrices S^1/2 and S^-1/2.
113 void S_half_invhalf(const arma::mat & S, arma::mat & Shalf, arma::mat & Sinvhalf, double cutoff);
114
115 /// Transform matrix to vector
116 arma::vec MatToVec(const arma::mat & v);
117 /// Transform vector to matrix
118 arma::mat VecToMat(const arma::vec & v, size_t nrows, size_t ncols);
119
120 /// Get vector from cube: c(i,j,:)
121 arma::vec slicevec(const arma::cube & c, size_t i, size_t j);
122
123 /// Compute cosine of matrix
124 arma::mat cosmat(const arma::mat & M);
125 /// Compute sine of matrix
126 arma::mat sinmat(const arma::mat & M);
127 /// Compute sin(x)/x of matrix
128 arma::mat sincmat(const arma::mat & M);
129 /// Compute square root of matrix
130 arma::mat sqrtmat(const arma::mat & M);
131 /// Compute exponential of matrix
132 arma::mat expmat(const arma::mat & M);
133 /// Compute exponential of matrix
134 arma::cx_mat expmat(const arma::cx_mat & M);
135
136 /// Check that the matrix is unitary
137 void check_unitarity(const arma::cx_mat & W);
138 /// Check that the matrix is orthogonal
139 void check_orthogonality(const arma::mat & W);
140
141 /// Orthogonalize
142 arma::mat orthogonalize(const arma::mat & M);
143 /// Unitarize
144 arma::cx_mat unitarize(const arma::cx_mat & M);
145
146 /// Orthonormalize vectors
147 arma::mat orthonormalize(const arma::mat & S, const arma::mat & C);
148 /// Orthonormalize vectors
149 arma::cx_mat orthonormalize(const arma::mat & S, const arma::cx_mat & C);
150
151 /**
152 * Find natural orbitals from P. Orbitals returned in decreasing occupation number.
153 */
154 void form_NOs(const arma::mat & P, const arma::mat & S, arma::mat & AO_to_NO, arma::mat & NO_to_AO, arma::vec & occs);
155
156 /**
157 * Find natural orbitals from P. Orbitals returned in decreasing occupation number.
158 */
159 void form_NOs(const arma::mat & P, const arma::mat & S, arma::mat & AO_to_NO, arma::vec & occs);
160
161 /**
162 * Does a pivoted Cholesky decomposition. The algorithm is adapted
163 * from the reference
164 *
165 * H. Harbrecht, M. Peters, and R. Schneider, "On the low-rank
166 * approximation by the pivoted Cholesky decomposition",
167 * Appl. Num. Math. 62, 428 (2012).
168 */
169 arma::mat pivoted_cholesky(const arma::mat & M, double thr, arma::uvec & pivot);
170 /// Same as above but doesn't return pivot
171 arma::mat pivoted_cholesky(const arma::mat & M, double thr);
172 /// Incomplete Cholesky factorization of matrix M, use n vectors
173 arma::mat incomplete_cholesky(const arma::mat & M, size_t n);
174
175 /// Transform B matrix to MO basis
176 arma::mat B_transform(arma::mat B, const arma::mat & Cl, const arma::mat & Cr);
177
178 /// Check thread safety of LAPACK library
179 void check_lapack_thread();
180
181 #endif
182