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 00038 #ifndef MAT_MATRIXHIERARCHICBASE 00039 #define MAT_MATRIXHIERARCHICBASE 00040 #include "matInclude.h" 00041 namespace mat{ 00048 template<class Treal, class Telement = Treal> 00049 class MatrixHierarchicBase { 00050 public: 00051 /* No public constructors (!) */ 00052 inline bool operator==(int k) const { 00053 if (k == 0) 00054 return this->is_zero(); 00055 else 00056 throw Failure("Matrix::operator== only implemented for k == 0"); 00057 } 00058 /* Check if matrix is zero (k = 0) */ 00059 #if 1 00060 inline const int& nScalarsRows() const 00061 {return rows.getNScalars();} 00062 inline const int& nScalarsCols() const 00063 {return cols.getNScalars();} 00064 #endif 00065 00066 inline const int& nrows() const /* Number of rows in Telement matrix */ 00067 {return rows.getNBlocks();} 00068 inline const int& ncols() const /* Number of cols in Telement matrix */ 00069 {return cols.getNBlocks();} 00070 00071 inline Telement& operator() /* Returns the element A(row, col) */ 00072 (int row, int col) { 00073 assert(elements); 00074 assert(row >= 0); 00075 assert(col >= 0); 00076 assert(row < nrows()); 00077 assert(col < ncols()); 00078 return elements[row + col * nrows()]; 00079 } 00080 inline const Telement& operator() /*Write protected reference returned*/ 00081 (int row, int col) const { 00082 assert(elements); 00083 assert(row >= 0); 00084 assert(col >= 0); 00085 assert(row < nrows()); 00086 assert(col < ncols()); 00087 return elements[row + col * nrows()]; 00088 } 00089 00090 inline Telement& operator[] 00091 (int index) { 00092 assert(elements); 00093 assert(index >= 0); 00094 assert(index < nElements()); 00095 return elements[index]; 00096 } 00097 inline Telement const & operator[] 00098 (int index) const { 00099 assert(elements); 00100 assert(index >= 0); 00101 assert(index < nElements()); 00102 return elements[index]; 00103 } 00104 00105 inline bool is_zero() const {return !elements;} 00106 00107 inline int nElements() const { 00108 return rows.getNBlocks() * cols.getNBlocks(); 00109 } 00110 00111 inline void resetRows(SizesAndBlocks const & newRows) { 00112 delete[] elements; 00113 elements = 0; 00114 rows = newRows; 00115 } 00116 inline void resetCols(SizesAndBlocks const & newCols) { 00117 delete[] elements; 00118 elements = 0; 00119 cols = newCols; 00120 } 00121 00122 inline void getRows(SizesAndBlocks & rowsCopy) const { 00123 rowsCopy = rows; 00124 } 00125 inline void getCols(SizesAndBlocks & colsCopy) const { 00126 colsCopy = cols; 00127 } 00128 00129 00130 inline bool highestLevel() const { 00131 return (rows.getNTotalScalars() == rows.getNScalars() && 00132 cols.getNTotalScalars() == cols.getNScalars()); 00133 } 00134 00140 inline bool is_empty() const { 00141 return rows.is_empty() || cols.is_empty(); 00142 } 00143 protected: 00144 00145 00146 MatrixHierarchicBase() 00147 : elements(0) {} 00148 MatrixHierarchicBase(SizesAndBlocks const & rowsInp, 00149 SizesAndBlocks const & colsInp) 00150 : rows(rowsInp), cols(colsInp), elements(0) {} 00151 MatrixHierarchicBase(const MatrixHierarchicBase<Treal, Telement>& mat); 00152 00153 00154 MatrixHierarchicBase<Treal, Telement>& 00155 operator=(const MatrixHierarchicBase<Treal, Telement>& mat); 00156 00157 static void swap(MatrixHierarchicBase<Treal, Telement>& A, 00158 MatrixHierarchicBase<Treal, Telement>& B); 00159 00160 virtual ~MatrixHierarchicBase(); 00161 SizesAndBlocks rows; 00162 SizesAndBlocks cols; 00163 Telement* elements; /* Length is nRows * nCols unless 0 */ 00164 00165 00166 00167 #if 0 00168 inline void assert_alloc() { 00169 if (this->cap < this->nel) { 00170 delete[] this->elements; 00171 this->cap = this->nel; 00172 this->elements = new Telement[this->cap]; 00173 for (int ind = 0; ind < this->cap; ind++) 00174 this->elements[ind] = 0; 00175 } 00176 } 00177 #endif 00178 private: 00179 00180 }; /* end class */ 00181 00182 00183 00184 template<class Treal, class Telement> /* Copy constructor */ 00185 MatrixHierarchicBase<Treal, Telement>:: 00186 MatrixHierarchicBase(const MatrixHierarchicBase<Treal, Telement>& mat) 00187 : rows(mat.rows), cols(mat.cols), elements(0) { 00188 if (!mat.is_zero()) { 00189 elements = new Telement[nElements()]; 00190 for (int i = 0; i < nElements(); i++) 00191 elements[i] = mat.elements[i]; 00192 } 00193 } 00194 00195 00196 template<class Treal, class Telement> /*Assignment operator*/ 00197 MatrixHierarchicBase<Treal, Telement>& 00198 MatrixHierarchicBase<Treal, Telement>:: 00199 operator=(const MatrixHierarchicBase<Treal, Telement>& mat) { 00200 if (mat.is_zero()) { /* Condition also matches empty matrices. */ 00201 rows = mat.rows; 00202 cols = mat.cols; 00203 delete[] elements; 00204 elements = 0; 00205 return *this; 00206 } 00207 if (is_zero() || (nElements() != mat.nElements())) { 00208 delete[] elements; 00209 elements = new Telement[mat.nElements()]; 00210 } 00211 rows = mat.rows; 00212 cols = mat.cols; 00213 for (int i = 0; i < nElements(); i++) 00214 elements[i] = mat.elements[i]; 00215 return *this; 00216 } 00217 00218 template<class Treal, class Telement> 00219 void MatrixHierarchicBase<Treal, Telement>:: 00220 swap(MatrixHierarchicBase<Treal, Telement>& A, 00221 MatrixHierarchicBase<Treal, Telement>& B) { 00222 assert(A.rows == B.rows && A.cols == B.cols); 00223 Telement* elementsTmp = A.elements; 00224 A.elements = B.elements; 00225 B.elements = elementsTmp; 00226 } 00227 00228 00229 template<class Treal, class Telement> 00230 MatrixHierarchicBase<Treal, Telement>::~MatrixHierarchicBase() { 00231 delete[] elements; 00232 elements = 0; 00233 } 00234 00235 00236 } /* end namespace mat */ 00237 00238 #endif