ergo
sparse_matrix.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28 #if !defined(_SPARSE_MATRIX_H_)
29 #define _SPARSE_MATRIX_H_ 1
30 
36 #include <stdio.h>
37 #include <vector>
38 #include <algorithm>
39 
40 
41 #include "realtype.h"
42 #include "matrix_typedefs.h"
43 #include "basisinfo.h"
44 #include "sparse_pattern.h"
45 
46 #if !defined(BEGIN_NAMESPACE)
47 #define BEGIN_NAMESPACE(x) namespace x {
48 #define END_NAMESPACE(x) } /* x */
49 #endif
50 
51 BEGIN_NAMESPACE(Dft)
52 
53 
54 class SparseMatrix {
55  class Exception : public std::exception {
56  const char *msg;
57  public:
58  explicit Exception(const char *msg_) : msg(msg_) {}
59  virtual const char *what() const throw() { return msg; }
60  };
61 
64  int **offsets;
65  int **his;
66  int *cnt;
67  int n;
69  void createOffsets(const SparsePattern& pattern);
70  public:
73  explicit SparseMatrix(const SparsePattern& pattern_);
74  SparseMatrix(const SparsePattern& pattern_,
75  const symmMatrix& m, const int *aoMap,
76  std::vector<int> const & permutationHML);
77 
79  for(int i=0; i<n; i++) {
80  delete [](columns[i]);
81  delete [](offsets[i]);
82  delete [](his[i]);
83  }
84  delete []columns;
85  delete []offsets;
86  delete []his;
87  delete []cnt;
88  }
89 
90  void print(const char *title) const;
91 
93  void addSymmetrizedTo(symmMatrix& sMat,
94  const int *aoMap,
95  std::vector<int> const & permutationHML) const;
96 
100  void add(int row, int col, ergo_real val) {
101  ergo_real *columnData = columns[col];
102  const int *hi = his[col];
103  int idx;
104  for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
105  //int idx = std::upper_bound(hi, hi+cnt[col], row)-hi;
106  if(idx >= cnt[col])
107  throw Exception("SparseMatrix::add called with incorrect args");
108  int offset = offsets[col][idx];
109  /* Add it... */
110  //assert(row-offset>=0);
111  //assert(row-offset<pattern.getColumnSize(col));
112  columnData[row-offset] += val;
113  }
114 
115  /* This operator[] syntax glue that we could in principle use can be
116  expensive performance-wise, do it the old-fashioned way.
117  Checking against intervals.end() is *terribly* expensive!!!
118  */
119  ergo_real at(int row, int col) const {
120  const ergo_real *columnData = columns[col];
121  const int *hi = his[col];
122  int idx; for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
123  if(idx >= cnt[col])
124  throw Exception("SparseMatrix::at called with incorrect args");
125  //int idx = std::upper_bound(hi, hi+cnt[col], row)-hi;
126  int offset = offsets[col][idx];
127  /* return it... */
128  //assert(row-offset>=0);
129  //assert(row-offset<pattern.getColumnSize(col));
130  return columnData[row-offset];
131  }
132 };
133 
134 
135 END_NAMESPACE(Dft)
136 
137 void
138 getrho_blocked_lda(int nbast, const Dft::SparseMatrix& dmat,
139  const ergo_real * gao,
140  const int* nblocks, const int (*iblocks)[2],
141  int ldaib, ergo_real *tmp, int nvclen, ergo_real *rho);
142 void
143 getrho_blocked_gga(int nbast, const Dft::SparseMatrix& dmat,
144  const ergo_real * gao,
145  const int* nblocks, const int (*iblocks)[2],
146  int ldaib, ergo_real *tmp, int nvclen,
147  ergo_real *rho, ergo_real (*grad)[3]);
148 
149 #endif /* _SPARSE_MATRIX_H_ */