ergo
MatrixTriangular.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
38#ifndef MAT_MATRIXTRIANGULAR
39#define MAT_MATRIXTRIANGULAR
40#include <stdexcept>
41#include "MatrixBase.h"
42namespace mat {
58 template<typename Treal, typename Tmatrix>
59 class MatrixTriangular : public MatrixBase<Treal, Tmatrix> {
60 public:
62 typedef Treal real;
63
65 :MatrixBase<Treal, Tmatrix>() {}
67 :MatrixBase<Treal, Tmatrix>(tri) {}
72 return *this;
73 }
74
76 *this->matrixPtr = k;
77 return *this;
78 }
86 (std::vector<int> const & rowind,
87 std::vector<int> const & colind,
88 std::vector<Treal> const & values,
89 SizesAndBlocks const & newRows,
90 SizesAndBlocks const & newCols) {
91 this->resetSizesAndBlocks(newRows, newCols);
92 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
93 }
106 inline void add_values
107 (std::vector<int> const & rowind,
108 std::vector<int> const & colind,
109 std::vector<Treal> const & values) {
110 this->matrixPtr->syAddValues(rowind, colind, values);
111 }
112
113
114 inline void get_values
115 (std::vector<int> const & rowind,
116 std::vector<int> const & colind,
117 std::vector<Treal> & values) const {
118 this->matrixPtr->syGetValues(rowind, colind, values);
119 }
127 inline void get_all_values
128 (std::vector<int> & rowind,
129 std::vector<int> & colind,
130 std::vector<Treal> & values) const {
131 rowind.resize(0);
132 colind.resize(0);
133 values.resize(0);
134 rowind.reserve(nnz());
135 colind.reserve(nnz());
136 values.reserve(nnz());
137 this->matrixPtr->syGetAllValues(rowind, colind, values);
138 }
144#if 0
145 inline void fullmatrix(Treal* const full, int const size) const {
146 this->matrixPtr->fullmatrix(full, size, size);
147 } /* FIXME? Should triangular matrix always have zeros below diagonal? */
148#endif
149
150 inline void inch(const MatrixGeneral<Treal, Tmatrix>& SPD,
151 const Treal threshold,
152 const side looking = left,
153 const inchversion version = unstable) {
154 Tmatrix::inch(*SPD.matrixPtr, *this->matrixPtr,
155 threshold, looking, version);
156 }
157 inline void inch(const MatrixSymmetric<Treal, Tmatrix>& SPD,
158 const Treal threshold,
159 const side looking = left,
160 const inchversion version = unstable) {
162 Tmatrix::syInch(*SPD.matrixPtr, *this->matrixPtr,
163 threshold, looking, version);
164 }
165
166 void thresh(Treal const threshold,
167 normType const norm);
168
169 inline Treal frob() const {
170 return this->matrixPtr->frob();
171 }
172 Treal eucl(Treal const requestedAccuracy,
173 int maxIter = -1) const;
174
175 Treal eucl_thresh(Treal const threshold);
176 Treal eucl_thresh_congr_trans_measure(Treal const threshold,
178 inline void frob_thresh(Treal threshold) {
179 this->matrixPtr->frob_thresh(threshold);
180 }
181 inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
182 return this->matrixPtr->nnz();
183 }
184 inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
185 return this->matrixPtr->nvalues();
186 }
187
188
189 inline void write_to_buffer(void* buffer, const int n_bytes) const {
190 this->write_to_buffer_base(buffer, n_bytes, matrix_triang);
191 }
192 inline void read_from_buffer(void* buffer, const int n_bytes) {
193 this->read_from_buffer_base(buffer, n_bytes, matrix_triang);
194 }
195
196 void random() {
197 this->matrixPtr->syRandom();
198 }
199
204 template<typename TRule>
205 void setElementsByRule(TRule & rule) {
206 this->matrixPtr->trSetElementsByRule(rule);
207 return;
208 }
209
213
214
215 std::string obj_type_id() const {return "MatrixTriangular";}
216 protected:
217 inline void writeToFileProt(std::ofstream & file) const {
218 this->writeToFileBase(file, matrix_triang);
219 }
220 inline void readFromFileProt(std::ifstream & file) {
221 this->readFromFileBase(file, matrix_triang);
222 }
223
224 private:
225
226 };
227
228 /* B += alpha * A */
229 template<typename Treal, typename Tmatrix>
230 inline MatrixTriangular<Treal, Tmatrix>&
233 assert(!sm.tB);
234 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
235 return *this;
236 }
237
238
239 template<typename Treal, typename Tmatrix>
241 thresh(Treal const threshold,
242 normType const norm) {
243 switch (norm) {
244 case frobNorm:
245 this->frob_thresh(threshold);
246 break;
247 default:
248 throw Failure("MatrixTriangular<Treal, Tmatrix>::"
249 "thresh(Treal const, "
250 "normType const): "
251 "Thresholding not imlpemented for selected norm");
252 }
253 }
254
255 template<typename Treal, typename Tmatrix>
257 eucl(Treal const requestedAccuracy,
258 int maxIter) const {
259 VectorType guess;
261 this->getCols(cols);
263 guess.rand();
265 if (maxIter < 0)
266 maxIter = this->get_nrows() * 100;
269 lan(ztz, guess, maxIter);
270 lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
271 lan.run();
272 Treal eVal = 0;
273 Treal acc = 0;
274 lan.getLargestMagnitudeEig(eVal, acc);
275 Interval<Treal> euclInt( template_blas_sqrt(eVal-acc),
276 template_blas_sqrt(eVal+acc) );
277 if ( euclInt.low() < 0 )
278 euclInt = Interval<Treal>( 0, template_blas_sqrt(eVal+acc) );
279 if ( euclInt.length() / 2.0 > requestedAccuracy ) {
280 std::cout << "req: " << (double)requestedAccuracy
281 << " obt: " << (double)(euclInt.length() / 2.0) << std::endl;
282 throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
283 }
284 return euclInt.midPoint();
285 }
286
287#if 1
288
289 template<typename Treal, typename Tmatrix>
291 eucl_thresh(Treal const threshold) {
293 return TruncObj.run( threshold );
294 }
295
296#endif
297
298 template<typename Treal, typename Tmatrix>
300 eucl_thresh_congr_trans_measure(Treal const threshold,
303 return TruncObj.run( threshold );
304 }
305
306
307} /* end namespace mat */
308#endif
309
310
Base class for matrix API.
Treal run(Treal const threshold)
Definition: truncation.h:80
Truncation of general matrices with impact on matrix triple multiply as error measure.
Definition: truncation.h:338
Truncation of general matrices.
Definition: truncation.h:272
Definition: Failure.h:57
Definition: Interval.h:46
Treal low() const
Definition: Interval.h:144
Treal midPoint() const
Definition: Interval.h:115
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Base class for matrix API.
Definition: MatrixBase.h:69
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:281
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:221
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:254
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:153
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:236
Normal matrix.
Definition: MatrixGeneral.h:59
Symmetric matrix.
Definition: MatrixSymmetric.h:68
Upper non-unit triangular matrix.
Definition: MatrixTriangular.h:59
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Assign from sparse matrix given by three vectors.
Definition: MatrixTriangular.h:86
void frob_thresh(Treal threshold)
Definition: MatrixTriangular.h:178
void setElementsByRule(TRule &rule)
Uses rule depending on the row and column indexes to set matrix elements The Trule class should have ...
Definition: MatrixTriangular.h:205
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixTriangular.h:220
MatrixTriangular< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &tri)
Definition: MatrixTriangular.h:70
void random()
Definition: MatrixTriangular.h:196
void add_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Add given set of values to the matrix (+=).
Definition: MatrixTriangular.h:107
std::string obj_type_id() const
Definition: MatrixTriangular.h:215
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixTriangular.h:61
void inch(const MatrixGeneral< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:150
Treal frob() const
Definition: MatrixTriangular.h:169
void inch(const MatrixSymmetric< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:157
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Get values given by row and column index lists.
Definition: MatrixTriangular.h:115
Treal eucl_thresh_congr_trans_measure(Treal const threshold, MatrixSymmetric< Treal, Tmatrix > &trA)
Definition: MatrixTriangular.h:300
Treal real
Definition: MatrixTriangular.h:62
void thresh(Treal const threshold, normType const norm)
Definition: MatrixTriangular.h:241
Treal eucl_thresh(Treal const threshold)
Definition: MatrixTriangular.h:291
void get_all_values(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &values) const
Get all values and corresponding row and column index lists, in matrix.
Definition: MatrixTriangular.h:128
MatrixTriangular()
Default constructor
Definition: MatrixTriangular.h:64
MatrixTriangular(const MatrixTriangular< Treal, Tmatrix > &tri)
Copy constructor
Definition: MatrixTriangular.h:66
MatrixTriangular< Treal, Tmatrix > & operator=(int const k)
Set matrix to zero or identity: A = 0 or A = 1.
Definition: MatrixTriangular.h:75
size_t nvalues() const
Definition: MatrixTriangular.h:184
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixTriangular.h:189
size_t nnz() const
Definition: MatrixTriangular.h:181
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixTriangular.h:257
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixTriangular.h:192
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixTriangular.h:217
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:99
Definition: VectorGeneral.h:48
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
void rand()
Definition: VectorGeneral.h:108
Definition: LanczosLargestMagnitudeEig.h:47
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
mat::SizesAndBlocks cols
Definition: test.cc:52
Definition: allocate.cc:39
side
Definition: Matrix.h:75
@ left
Definition: Matrix.h:75
@ matrix_triang
Definition: MatrixBase.h:56
normType
Definition: matInclude.h:139
@ frobNorm
Definition: matInclude.h:139
inchversion
Definition: Matrix.h:76
@ unstable
Definition: Matrix.h:76
Definition: mat_utils.h:78
This proxy expresses the result of multiplication of two objects, of possibly different types,...
Definition: matrix_proxy.h:51
Treal template_blas_sqrt(Treal x)