MatrixBase.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 
00036 #ifndef MAT_MATRIXBASE
00037 #define MAT_MATRIXBASE
00038 #include <iostream>
00039 #include <fstream>
00040 #include <ios>
00041 #include "FileWritable.h"
00042 #include "matrix_proxy.h"
00043 #include "ValidPtr.h"
00044 #include "SizesAndBlocks.h"
00045 namespace mat {
00046   template<typename Treal, typename Tmatrix>
00047     class MatrixGeneral; 
00048   template<typename Treal, typename Tmatrix>
00049     class MatrixSymmetric;  
00050   template<typename Treal, typename Tmatrix>
00051     class MatrixTriangular; 
00052   template<typename Treal, typename Tvector>
00053     class VectorGeneral;
00054   enum matrix_type {matrix_matr, matrix_symm, matrix_triang};
00055   
00066   template<typename Treal, typename Tmatrix>
00067     class MatrixBase : public FileWritable {
00068   public:
00069     friend class MatrixGeneral<Treal, Tmatrix>;
00070     friend class MatrixSymmetric<Treal, Tmatrix>;
00071     friend class MatrixTriangular<Treal, Tmatrix>;
00072     
00073 
00074     inline void resetSizesAndBlocks(SizesAndBlocks const & newRows,
00075                                     SizesAndBlocks const & newCols) {
00076       matrixPtr.haveDataStructureSet(true);
00077       matrixPtr->resetRows(newRows);
00078       matrixPtr->resetCols(newCols);
00079     }
00080     inline void getRows(SizesAndBlocks & rowsCopy) const {
00081       matrixPtr->getRows(rowsCopy);
00082     }
00083     inline void getCols(SizesAndBlocks & colsCopy) const {
00084       matrixPtr->getCols(colsCopy);
00085     }
00086 
00091     inline bool is_empty() const {
00092       return !matrixPtr.haveDataStructureGet();
00093     }
00094 
00095     inline Treal trace() const {
00096       return matrixPtr->trace();
00097     }
00098     
00099     inline void add_identity(Treal alpha) {
00100       matrixPtr->addIdentity(alpha);
00101     }
00102     inline MatrixBase<Treal, Tmatrix>& operator*=(Treal const alpha) {
00103       *matrixPtr *= alpha;
00104       return *this;
00105     }
00106 
00107     inline bool operator==(int k) const {
00108       if (k == 0)
00109         return *matrixPtr == 0;
00110       else 
00111         throw Failure("MatrixBase::operator== only implemented for k == 0");
00112     }
00113 
00114     
00115     
00116     inline void clear() {
00117       if (is_empty())
00118         // This means that the object's data structure has not been set
00119         // There is nothing to clear and the matrixPtr is not valid either
00120         return;
00121       matrixPtr->clear();
00122     }
00123     
00124     inline size_t memory_usage() const {
00125       return matrixPtr->memory_usage();
00126     }
00127 
00128     inline void write_to_buffer_count(int& n_bytes) const {
00129       int ib_length = 3;
00130       int vb_length = 0;
00131       this->matrixPtr->write_to_buffer_count(ib_length, vb_length);
00132       n_bytes = vb_length * sizeof(Treal) + ib_length * sizeof(int);
00133     }
00134 
00135 #if 1
00136     inline int get_nrows() const {
00137       return matrixPtr->nScalarsRows();
00138     }
00139     inline int get_ncols() const {
00140       return matrixPtr->nScalarsCols();
00141     }
00142 #endif
00143     
00144     inline Tmatrix const & getMatrix() const {return *matrixPtr;}
00145     inline Tmatrix & getMatrix() {return *matrixPtr;}
00146  
00148     inline Treal maxAbsValue() const {return matrixPtr->maxAbsValue();}
00149 
00150   protected:
00151     ValidPtr<Tmatrix> matrixPtr;
00152     
00153     MatrixBase():matrixPtr(new Tmatrix) {}
00154     MatrixBase(const MatrixBase<Treal, Tmatrix>& other)
00155       :FileWritable(other), matrixPtr(new Tmatrix) {
00156       matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet());
00157       /* getConstRefForCopying() is used here to make sure it works
00158          also in the case when the matrix is written to file. */
00159       *matrixPtr = other.matrixPtr.getConstRefForCopying();
00160       matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet());
00161     }
00162       
00163     MatrixBase<Treal, Tmatrix>& 
00164       operator=(const MatrixBase<Treal, Tmatrix>& other) {
00165       FileWritable::operator=(other); /* Allows us to copy mat on file */
00166       matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet());
00167       /* getConstRefForCopying() is used here to make sure it works
00168          also in the case when the matrix is written to file. */
00169       *matrixPtr = other.matrixPtr.getConstRefForCopying();
00170       matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet());
00171       return *this;
00172     } 
00173     
00174     MatrixBase<Treal, Tmatrix>& 
00175       operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) {
00176       if (mt.A.matrixPtr.haveDataStructureGet()) {
00177         matrixPtr.haveDataStructureSet(true);
00178       }
00179       if (mt.tA)
00180         Tmatrix::transpose(*mt.A.matrixPtr, *this->matrixPtr);
00181       else
00182         *this->matrixPtr = *mt.A.matrixPtr; 
00183       return *this;
00184       // FileWritable::operator=(other);/*Could be used to copy mat on file*/
00185     } 
00186     
00187 
00188     void write_to_buffer_base(void* buffer, const int n_bytes,
00189                          const matrix_type mattype) const;     
00190     void read_from_buffer_base(void* buffer, const int n_bytes,
00191                           const matrix_type mattype);
00192 
00193     void writeToFileBase(std::ofstream & file, 
00194                          matrix_type const mattype) const;
00195     void readFromFileBase(std::ifstream & file, 
00196                           matrix_type const mattype);
00197 
00198     std::string obj_type_id() const {return "MatrixBase";}
00199     inline void inMemorySet(bool inMem) {
00200       matrixPtr.inMemorySet(inMem);
00201     }
00202 
00203     static void getPermutedIndexes(std::vector<int> const & index, 
00204                                    std::vector<int> const & permutation,
00205                                    std::vector<int> & newIndex) {
00206       newIndex.resize(index.size());
00207       for (unsigned int i = 0; i < index.size(); ++i) 
00208         newIndex[i] = permutation[index[i]];
00209     }
00210 
00211 
00212   private:
00213     
00214   };
00215 
00216    
00217   template<typename Treal, typename Tmatrix>
00218     void MatrixBase<Treal, Tmatrix>::
00219     writeToFileBase(std::ofstream & file, 
00220                 matrix_type const mattype) const {
00221     int type = (int)mattype;
00222     file.write((char*)&type,sizeof(int));
00223     
00224     if (is_empty())
00225       // This means that the object's data structure has not been set
00226       // The ValidPtr prevents setting the data structure between
00227       // calls to writeToFile and readFromFile 
00228       return;
00229     matrixPtr->writeToFile(file);
00230   }
00231   
00232   template<typename Treal, typename Tmatrix>
00233     void MatrixBase<Treal, Tmatrix>::
00234     readFromFileBase(std::ifstream & file, 
00235                  matrix_type const mattype) {
00236     char type[sizeof(int)];
00237     file.read(type, sizeof(int));
00238     if (((int)*type) != mattype)
00239       throw Failure("MatrixBase<Treal, Tmatrix>::"
00240                     "readFromFile(std::ifstream &, " 
00241                     "matrix_type const): Wrong matrix type");
00242     if (is_empty())
00243       // This means that the object's data structure has not been set
00244       return;
00245     matrixPtr->readFromFile(file);
00246   }
00247 
00248 
00249 
00250   template<typename Treal, typename Tmatrix>
00251     void MatrixBase<Treal, Tmatrix>::
00252     write_to_buffer_base(void* buffer, const int n_bytes,
00253                          const matrix_type mattype) const {
00254     int ib_length = 3; /* Length of integer buffer, at least 3: matrix_type, */
00255     /*                    ib_length and vb_length                            */
00256     int vb_length = 0; /* Length of value buffer                             */
00257     this->matrixPtr->write_to_buffer_count(ib_length, vb_length);
00258     if (n_bytes >= 
00259         (int)(vb_length * sizeof(Treal) + ib_length * sizeof(int))) {
00260       int* int_buf = (int*)buffer;
00261       int_buf[0] = mattype;
00262       int_buf[1] = ib_length;
00263       int_buf[2] = vb_length;
00264       Treal* value_buf = (Treal*)&(int_buf[ib_length]); /* Value buffer      */
00265       /* begins after integer buffer end                                     */
00266       int ib_index = 0; 
00267       int vb_index = 0;
00268       this->matrixPtr->write_to_buffer(&int_buf[3], ib_length - 3, 
00269                                        value_buf, vb_length,
00270                                        ib_index, vb_index);
00271     }
00272     else {
00273       throw Failure("MatrixBase::write_to_buffer: Buffer is too small");
00274     }
00275   }
00276     
00277   template<typename Treal, typename Tmatrix>
00278     void MatrixBase<Treal, Tmatrix>::
00279     read_from_buffer_base(void* buffer, const int n_bytes,
00280                      const matrix_type mattype) {
00281     int* int_buf = (int*)buffer;
00282     if(int_buf[0] == mattype) {
00283       int ib_length = int_buf[1];
00284       int vb_length = int_buf[2];
00285       int ib_index = 0; 
00286       int vb_index = 0;
00287       Treal* value_buf = (Treal*)&(int_buf[ib_length]);
00288       this->matrixPtr->read_from_buffer(&int_buf[3], ib_length - 3,
00289                                         value_buf, vb_length,
00290                                         ib_index, vb_index);
00291     }
00292     else {
00293       throw Failure("MatrixBase::read_from_buffer: Wrong matrix type");
00294     }
00295   }
00296 
00297 
00298 } /* end namespace mat */
00299 #endif
00300 
00301 

Generated on Wed Nov 21 09:32:01 2012 for ergo by  doxygen 1.4.7