MatrixHierarchicBase.h

Go to the documentation of this file.
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

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1