sparse_matrix.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 
00028 #if !defined(_SPARSE_MATRIX_H_)
00029 #define _SPARSE_MATRIX_H_ 1
00030 
00036 #include <stdio.h>
00037 #include <vector>
00038 #include <algorithm>
00039 
00040 
00041 #include "realtype.h"
00042 #include "matrix_typedefs.h"
00043 #include "basisinfo.h"
00044 #include "sparse_pattern.h"
00045 
00046 #if !defined(BEGIN_NAMESPACE)
00047 #define BEGIN_NAMESPACE(x) namespace x {
00048 #define END_NAMESPACE(x)   } /* x */
00049 #endif
00050 
00051 BEGIN_NAMESPACE(Dft)
00052 
00053 
00054 class SparseMatrix {
00055   class Exception : public std::exception {
00056     const char *msg;
00057   public:
00058   explicit Exception(const char *msg_) : msg(msg_) {}
00059     virtual const char *what() const throw() { return msg; }
00060   };
00061 
00062   const SparsePattern& pattern;
00063   ergo_real **columns;
00064   int **offsets; 
00065   int **his;     
00066   int *cnt;      
00067   int n;
00069   void createOffsets(const SparsePattern& pattern);
00070  public:
00073   explicit SparseMatrix(const SparsePattern& pattern_);
00074   SparseMatrix(const SparsePattern& pattern_,
00075                const symmMatrix& m, const int *aoMap,
00076                std::vector<int> const & permutationHML);
00077 
00078   ~SparseMatrix() {
00079     for(int i=0; i<n; i++) {
00080       delete [](columns[i]);
00081       delete [](offsets[i]);
00082       delete [](his[i]);
00083     }
00084     delete []columns;
00085     delete []offsets;
00086     delete []his;
00087     delete []cnt;
00088   }
00089 
00090   void print(const char *title) const;
00091 
00093   void addSymmetrizedTo(symmMatrix& sMat,
00094                         const int *aoMap,
00095                         std::vector<int> const & permutationHML) const;
00096 
00100   void add(int row, int col, ergo_real val) {
00101     ergo_real *columnData = columns[col];
00102     const int *hi = his[col];
00103     int idx;
00104     for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
00105     //int idx = std::upper_bound(hi, hi+cnt[col], row)-hi;
00106     if(idx >= cnt[col])
00107       throw Exception("SparseMatrix::add called with incorrect args");
00108     int offset = offsets[col][idx];
00109     /* Add it... */
00110     //assert(row-offset>=0);
00111     //assert(row-offset<pattern.getColumnSize(col));
00112     columnData[row-offset] += val;
00113   }
00114 
00115   /* This operator[] syntax glue that we could in principle use can be
00116      expensive performance-wise, do it the old-fashioned way.
00117      Checking against intervals.end() is *terribly* expensive!!!
00118   */
00119   ergo_real at(int row, int col) const {
00120     const ergo_real *columnData = columns[col];
00121     const int *hi = his[col]; 
00122     int idx; for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
00123     if(idx >= cnt[col])
00124       throw Exception("SparseMatrix::at called with incorrect args");
00125     //int idx = std::upper_bound(hi, hi+cnt[col], row)-hi;
00126     int offset = offsets[col][idx];
00127     /* return it... */
00128     //assert(row-offset>=0);
00129     //assert(row-offset<pattern.getColumnSize(col));
00130     return columnData[row-offset];
00131   }
00132 };
00133 
00134 
00135 END_NAMESPACE(Dft)
00136 
00137 void
00138 getrho_blocked_lda(int nbast, const Dft::SparseMatrix& dmat,
00139                    const ergo_real * gao,
00140                    const int* nblocks, const int (*iblocks)[2],
00141                    int ldaib, ergo_real *tmp, int nvclen, ergo_real *rho);
00142 void
00143 getrho_blocked_gga(int nbast, const Dft::SparseMatrix& dmat,
00144                    const ergo_real * gao,
00145                    const int* nblocks, const int (*iblocks)[2],
00146                    int ldaib, ergo_real *tmp, int nvclen,
00147                    ergo_real *rho, ergo_real (*grad)[3]);
00148 
00149 #endif /* _SPARSE_MATRIX_H_ */

Generated on Mon Sep 17 14:30:39 2012 for ergo by  doxygen 1.4.7