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 00036 #ifndef MAT_MATRIXTRIANGULAR 00037 #define MAT_MATRIXTRIANGULAR 00038 #include <stdexcept> 00039 #include "MatrixBase.h" 00040 namespace mat { 00056 template<typename Treal, typename Tmatrix> 00057 class MatrixTriangular : public MatrixBase<Treal, Tmatrix> { 00058 public: 00059 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType; 00060 00061 MatrixTriangular() 00062 :MatrixBase<Treal, Tmatrix>() {} 00063 explicit MatrixTriangular(const MatrixTriangular<Treal, Tmatrix>& tri) 00064 :MatrixBase<Treal, Tmatrix>(tri) {} 00066 MatrixTriangular<Treal, Tmatrix>& 00067 operator=(const MatrixTriangular<Treal, Tmatrix>& tri) { 00068 MatrixBase<Treal, Tmatrix>::operator=(tri); 00069 return *this; 00070 } 00071 00072 inline MatrixTriangular<Treal, Tmatrix>& operator=(int const k) { 00073 *this->matrixPtr = k; 00074 return *this; 00075 } 00082 inline void assign_from_sparse 00083 (std::vector<int> const & rowind, 00084 std::vector<int> const & colind, 00085 std::vector<Treal> const & values, 00086 SizesAndBlocks const & newRows, 00087 SizesAndBlocks const & newCols) { 00088 this->resetSizesAndBlocks(newRows, newCols); 00089 this->matrixPtr->syAssignFromSparse(rowind, colind, values); 00090 } 00103 inline void add_values 00104 (std::vector<int> const & rowind, 00105 std::vector<int> const & colind, 00106 std::vector<Treal> const & values) { 00107 this->matrixPtr->syAddValues(rowind, colind, values); 00108 } 00109 00110 00111 inline void get_values 00112 (std::vector<int> const & rowind, 00113 std::vector<int> const & colind, 00114 std::vector<Treal> & values) const { 00115 this->matrixPtr->syGetValues(rowind, colind, values); 00116 } 00124 inline void get_all_values 00125 (std::vector<int> & rowind, 00126 std::vector<int> & colind, 00127 std::vector<Treal> & values) const { 00128 this->matrixPtr->syGetAllValues(rowind, colind, values); 00129 } 00135 #if 0 00136 inline void fullmatrix(Treal* const full, int const size) const { 00137 this->matrixPtr->fullmatrix(full, size, size); 00138 } /* FIXME? Should triangular matrix always have zeros below diagonal? */ 00139 #endif 00140 00141 inline void inch(const MatrixGeneral<Treal, Tmatrix>& SPD, 00142 const Treal threshold, 00143 const side looking = left, 00144 const inchversion version = unstable) { 00145 Tmatrix::inch(*SPD.matrixPtr, *this->matrixPtr, 00146 threshold, looking, version); 00147 } 00148 inline void inch(const MatrixSymmetric<Treal, Tmatrix>& SPD, 00149 const Treal threshold, 00150 const side looking = left, 00151 const inchversion version = unstable) { 00152 this->matrixPtr.haveDataStructureSet(true); 00153 Tmatrix::syInch(*SPD.matrixPtr, *this->matrixPtr, 00154 threshold, looking, version); 00155 } 00156 00157 void thresh(Treal const threshold, 00158 normType const norm); 00159 00160 inline Treal frob() const { 00161 return this->matrixPtr->frob(); 00162 } 00163 Treal eucl(Treal const requestedAccuracy, 00164 int maxIter = -1) const; 00165 00166 Treal eucl_thresh(Treal const threshold); 00167 Treal eucl_thresh_congr_trans_measure(Treal const threshold, 00168 MatrixSymmetric<Treal, Tmatrix> & trA); 00169 inline void frob_thresh(Treal threshold) { 00170 this->matrixPtr->frob_thresh(threshold); 00171 } 00172 inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */ 00173 return this->matrixPtr->nnz(); 00174 } 00175 inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */ 00176 return this->matrixPtr->nvalues(); 00177 } 00178 00179 00180 inline void write_to_buffer(void* buffer, const int n_bytes) const { 00181 this->write_to_buffer_base(buffer, n_bytes, matrix_triang); 00182 } 00183 inline void read_from_buffer(void* buffer, const int n_bytes) { 00184 this->read_from_buffer_base(buffer, n_bytes, matrix_triang); 00185 } 00186 00187 void random() { 00188 this->matrixPtr->syRandom(); 00189 } 00190 00195 template<typename TRule> 00196 void setElementsByRule(TRule & rule) { 00197 this->matrixPtr->trSetElementsByRule(rule); 00198 return; 00199 } 00200 00202 MatrixTriangular<Treal, Tmatrix>& operator+= 00203 (XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm); 00204 00205 00206 std::string obj_type_id() const {return "MatrixTriangular";} 00207 protected: 00208 inline void writeToFileProt(std::ofstream & file) const { 00209 this->writeToFileBase(file, matrix_triang); 00210 } 00211 inline void readFromFileProt(std::ifstream & file) { 00212 this->readFromFileBase(file, matrix_triang); 00213 } 00214 00215 private: 00216 00217 }; 00218 00219 /* B += alpha * A */ 00220 template<typename Treal, typename Tmatrix> 00221 inline MatrixTriangular<Treal, Tmatrix>& 00222 MatrixTriangular<Treal, Tmatrix>::operator+= 00223 (XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm) { 00224 assert(!sm.tB); 00225 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr); 00226 return *this; 00227 } 00228 00229 00230 template<typename Treal, typename Tmatrix> 00231 void MatrixTriangular<Treal, Tmatrix>:: 00232 thresh(Treal const threshold, 00233 normType const norm) { 00234 switch (norm) { 00235 case frobNorm: 00236 this->frob_thresh(threshold); 00237 break; 00238 default: 00239 throw Failure("MatrixTriangular<Treal, Tmatrix>::" 00240 "thresh(Treal const, " 00241 "normType const): " 00242 "Thresholding not imlpemented for selected norm"); 00243 } 00244 } 00245 00246 template<typename Treal, typename Tmatrix> 00247 Treal MatrixTriangular<Treal, Tmatrix>:: 00248 eucl(Treal const requestedAccuracy, 00249 int maxIter) const { 00250 VectorType guess; 00251 SizesAndBlocks cols; 00252 this->getCols(cols); 00253 guess.resetSizesAndBlocks(cols); 00254 guess.rand(); 00255 mat::ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal> ztz(*this); 00256 if (maxIter < 0) 00257 maxIter = this->get_nrows() * 100; 00258 arn::LanczosLargestMagnitudeEig 00259 <Treal, ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal>, VectorType> 00260 lan(ztz, guess, maxIter); 00261 lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() ); 00262 lan.run(); 00263 Treal eVal = 0; 00264 Treal acc = 0; 00265 lan.getLargestMagnitudeEig(eVal, acc); 00266 Interval<Treal> euclInt( template_blas_sqrt(eVal-acc), 00267 template_blas_sqrt(eVal+acc) ); 00268 if ( euclInt.low() < 0 ) 00269 euclInt = Interval<Treal>( 0, template_blas_sqrt(eVal+acc) ); 00270 if ( euclInt.length() / 2.0 > requestedAccuracy ) { 00271 std::cout << "req: " << requestedAccuracy 00272 << " obt: " << euclInt.length() / 2.0 << std::endl; 00273 throw std::runtime_error("Desired accuracy not obtained in Lanczos."); 00274 } 00275 return euclInt.midPoint(); 00276 } 00277 00278 #if 1 00279 00280 template<typename Treal, typename Tmatrix> 00281 Treal MatrixTriangular<Treal, Tmatrix>:: 00282 eucl_thresh(Treal const threshold) { 00283 EuclTruncationGeneral<MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj( *this ); 00284 return TruncObj.run( threshold ); 00285 } 00286 00287 #endif 00288 00289 template<typename Treal, typename Tmatrix> 00290 Treal MatrixTriangular<Treal, Tmatrix>:: 00291 eucl_thresh_congr_trans_measure(Treal const threshold, 00292 MatrixSymmetric<Treal, Tmatrix> & trA) { 00293 EuclTruncationCongrTransMeasure<MatrixTriangular<Treal, Tmatrix>, MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this, trA); 00294 return TruncObj.run( threshold ); 00295 } 00296 00297 00298 } /* end namespace mat */ 00299 #endif 00300 00301