ergo
MatrixTriangular.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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"
42 namespace 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  }
85  inline void assign_from_sparse
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) {
161  this->matrixPtr.haveDataStructureSet(true);
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>&
231  MatrixTriangular<Treal, Tmatrix>::operator+=
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);
262  guess.resetSizesAndBlocks(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 
void thresh(Treal const threshold, normType const norm)
Definition: MatrixTriangular.h:241
MatrixTriangular< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &tri)
Definition: MatrixTriangular.h:70
Treal eucl_thresh(Treal const threshold)
Definition: MatrixTriangular.h:291
Normal matrix.
Definition: MatrixBase.h:49
mat::SizesAndBlocks cols
Definition: test.cc:52
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
Treal low() const
Definition: Interval.h:144
inchversion
Definition: Matrix.h:76
void frob_thresh(Treal threshold)
Definition: MatrixTriangular.h:178
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:153
Definition: LanczosLargestMagnitudeEig.h:46
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:281
void inch(const MatrixGeneral< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:150
Truncation of general matrices.
Definition: truncation.h:272
Definition: MatrixBase.h:55
Treal real
Definition: MatrixTriangular.h:62
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
MatrixTriangular()
Default constructor.
Definition: MatrixTriangular.h:64
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixTriangular.h:257
Definition: Interval.h:46
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixTriangular.h:61
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixTriangular.h:189
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:236
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 writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixTriangular.h:217
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:254
size_t nnz() const
Definition: MatrixTriangular.h:181
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
MatrixTriangular(const MatrixTriangular< Treal, Tmatrix > &tri)
Copy constructor.
Definition: MatrixTriangular.h:66
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixTriangular.h:220
Treal run(Treal const threshold)
Definition: truncation.h:80
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:221
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
Treal midPoint() const
Definition: Interval.h:115
void rand()
Definition: VectorGeneral.h:102
void inch(const MatrixSymmetric< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:157
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
size_t nvalues() const
Definition: MatrixTriangular.h:184
Treal eucl_thresh_congr_trans_measure(Treal const threshold, MatrixSymmetric< Treal, Tmatrix > &trA)
Definition: MatrixTriangular.h:300
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
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
Treal frob() const
Definition: MatrixTriangular.h:169
Definition: mat_utils.h:78
void random()
Definition: MatrixTriangular.h:196
Base class for matrix API.
Definition: MatrixBase.h:69
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixTriangular.h:192
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
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
side
Definition: Matrix.h:75
Base class for matrix API.
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Definition: matInclude.h:139
Definition: Failure.h:57
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:51
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:99
Definition: Matrix.h:75
Definition: MatrixBase.h:56
MatrixTriangular< Treal, Tmatrix > & operator=(int const k)
Set matrix to zero or identity: A = 0 or A = 1.
Definition: MatrixTriangular.h:75
Symmetric matrix.
Definition: MatrixBase.h:51
std::string obj_type_id() const
Definition: MatrixTriangular.h:215
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
Truncation of general matrices with impact on matrix triple multiply as error measure.
Definition: truncation.h:338
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:139
Definition: Matrix.h:76