1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
/* Ergo, version 3.8.2, a program for linear scaling electronic structure
* calculations.
* Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
* and Anastasia Kruchinina.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Primary academic reference:
* Ergo: An open-source program for linear-scaling electronic structure
* calculations,
* Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
* Kruchinina,
* SoftwareX 7, 107 (2018),
* <http://dx.doi.org/10.1016/j.softx.2018.03.005>
*
* For further information about Ergo, see <http://www.ergoscf.org>.
*/
#if !defined(_SPARSE_MATRIX_H_)
#define _SPARSE_MATRIX_H_ 1
/** @file sparse_matrix.h Declares a sparse matrix optimized for the
XC code. The object provides methods for fast preallocation of
the matrix elements, and some matrix elements iterators as needed
for the numerical matrix element integration schemes. */
#include <stdio.h>
#include <vector>
#include <algorithm>
#include "realtype.h"
#include "matrix_typedefs.h"
#include "basisinfo.h"
#include "sparse_pattern.h"
#if !defined(BEGIN_NAMESPACE)
#define BEGIN_NAMESPACE(x) namespace x {
#define END_NAMESPACE(x) } /* x */
#endif
BEGIN_NAMESPACE(Dft)
/** Sparse matrix structure optimized for XC data access pattern. */
class SparseMatrix {
class Exception : public std::exception {
const char *msg;
public:
explicit Exception(const char *msg_) : msg(msg_) {}
virtual const char *what() const throw() { return msg; }
};
const SparsePattern& pattern;
ergo_real **columns;
int **offsets; /**< for accelerated at() and add() methods. */
int **his; /**< for accelerated at() and add() methods. */
int *cnt; /**< for accelerated at() and add() methods. */
int n;
/** Fills in offsets and his based on pattern. */
void createOffsets(const SparsePattern& pattern);
public:
/** Constructs a square matrix and preallocate according to the
specified pattern. */
explicit SparseMatrix(const SparsePattern& pattern_);
SparseMatrix(const SparsePattern& pattern_,
const symmMatrix& m, const int *aoMap,
std::vector<int> const & permutationHML);
~SparseMatrix() {
for(int i=0; i<n; i++) {
delete [](columns[i]);
delete [](offsets[i]);
delete [](his[i]);
}
delete []columns;
delete []offsets;
delete []his;
delete []cnt;
}
void print(const char *title) const;
/** Assigns itself to a given hierarchic matrix. */
void addSymmetrizedTo(symmMatrix& sMat,
const int *aoMap,
std::vector<int> const & permutationHML) const;
/** Adds given value to an element in given row and column.
Checking against intervals.end() is *terribly* expensive!!!
Luckily, we do not have to do it.*/
void add(int row, int col, ergo_real val) {
ergo_real *columnData = columns[col];
const int *hi = his[col];
int idx;
for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
//int idx = std::upper_bound(hi, hi+cnt[col], row)-hi;
if(idx >= cnt[col])
throw Exception("SparseMatrix::add called with incorrect args");
int offset = offsets[col][idx];
/* Add it... */
//assert(row-offset>=0);
//assert(row-offset<pattern.getColumnSize(col));
columnData[row-offset] += val;
}
/* This operator[] syntax glue that we could in principle use can be
expensive performance-wise, do it the old-fashioned way.
Checking against intervals.end() is *terribly* expensive!!!
*/
ergo_real at(int row, int col) const {
const ergo_real *columnData = columns[col];
const int *hi = his[col];
int idx; for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
if(idx >= cnt[col])
throw Exception("SparseMatrix::at called with incorrect args");
//int idx = std::upper_bound(hi, hi+cnt[col], row)-hi;
int offset = offsets[col][idx];
/* return it... */
//assert(row-offset>=0);
//assert(row-offset<pattern.getColumnSize(col));
return columnData[row-offset];
}
};
END_NAMESPACE(Dft)
void
getrho_blocked_lda(int nbast, const Dft::SparseMatrix& dmat,
const ergo_real * gao,
const int* nblocks, const int (*iblocks)[2],
int ldaib, ergo_real *tmp, int nvclen, ergo_real *rho);
void
getrho_blocked_gga(int nbast, const Dft::SparseMatrix& dmat,
const ergo_real * gao,
const int* nblocks, const int (*iblocks)[2],
int ldaib, ergo_real *tmp, int nvclen,
ergo_real *rho, ergo_real (*grad)[3]);
#endif /* _SPARSE_MATRIX_H_ */
|
