ergo
MatrixTriangular.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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 
36 #ifndef MAT_MATRIXTRIANGULAR
37 #define MAT_MATRIXTRIANGULAR
38 #include <stdexcept>
39 #include "MatrixBase.h"
40 namespace mat {
56  template<typename Treal, typename Tmatrix>
57  class MatrixTriangular : public MatrixBase<Treal, Tmatrix> {
58  public:
60 
62  :MatrixBase<Treal, Tmatrix>() {}
64  :MatrixBase<Treal, Tmatrix>(tri) {}
69  return *this;
70  }
71 
73  *this->matrixPtr = k;
74  return *this;
75  }
82  inline void assign_from_sparse
83  (std::vector<int> const & rowind,
84  std::vector<int> const & colind,
85  std::vector<Treal> const & values,
86  SizesAndBlocks const & newRows,
87  SizesAndBlocks const & newCols) {
88  this->resetSizesAndBlocks(newRows, newCols);
89  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
90  }
103  inline void add_values
104  (std::vector<int> const & rowind,
105  std::vector<int> const & colind,
106  std::vector<Treal> const & values) {
107  this->matrixPtr->syAddValues(rowind, colind, values);
108  }
109 
110 
111  inline void get_values
112  (std::vector<int> const & rowind,
113  std::vector<int> const & colind,
114  std::vector<Treal> & values) const {
115  this->matrixPtr->syGetValues(rowind, colind, values);
116  }
124  inline void get_all_values
125  (std::vector<int> & rowind,
126  std::vector<int> & colind,
127  std::vector<Treal> & values) const {
128  this->matrixPtr->syGetAllValues(rowind, colind, values);
129  }
135 #if 0
136  inline void fullmatrix(Treal* const full, int const size) const {
137  this->matrixPtr->fullmatrix(full, size, size);
138  } /* FIXME? Should triangular matrix always have zeros below diagonal? */
139 #endif
140 
141  inline void inch(const MatrixGeneral<Treal, Tmatrix>& SPD,
142  const Treal threshold,
143  const side looking = left,
144  const inchversion version = unstable) {
145  Tmatrix::inch(*SPD.matrixPtr, *this->matrixPtr,
146  threshold, looking, version);
147  }
148  inline void inch(const MatrixSymmetric<Treal, Tmatrix>& SPD,
149  const Treal threshold,
150  const side looking = left,
151  const inchversion version = unstable) {
152  this->matrixPtr.haveDataStructureSet(true);
153  Tmatrix::syInch(*SPD.matrixPtr, *this->matrixPtr,
154  threshold, looking, version);
155  }
156 
157  void thresh(Treal const threshold,
158  normType const norm);
159 
160  inline Treal frob() const {
161  return this->matrixPtr->frob();
162  }
163  Treal eucl(Treal const requestedAccuracy,
164  int maxIter = -1) const;
165 
166  Treal eucl_thresh(Treal const threshold);
167  Treal eucl_thresh_congr_trans_measure(Treal const threshold,
169  inline void frob_thresh(Treal threshold) {
170  this->matrixPtr->frob_thresh(threshold);
171  }
172  inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
173  return this->matrixPtr->nnz();
174  }
175  inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
176  return this->matrixPtr->nvalues();
177  }
178 
179 
180  inline void write_to_buffer(void* buffer, const int n_bytes) const {
181  this->write_to_buffer_base(buffer, n_bytes, matrix_triang);
182  }
183  inline void read_from_buffer(void* buffer, const int n_bytes) {
184  this->read_from_buffer_base(buffer, n_bytes, matrix_triang);
185  }
186 
187  void random() {
188  this->matrixPtr->syRandom();
189  }
190 
195  template<typename TRule>
196  void setElementsByRule(TRule & rule) {
197  this->matrixPtr->trSetElementsByRule(rule);
198  return;
199  }
200 
204 
205 
206  std::string obj_type_id() const {return "MatrixTriangular";}
207  protected:
208  inline void writeToFileProt(std::ofstream & file) const {
209  this->writeToFileBase(file, matrix_triang);
210  }
211  inline void readFromFileProt(std::ifstream & file) {
212  this->readFromFileBase(file, matrix_triang);
213  }
214 
215  private:
216 
217  };
218 
219  /* B += alpha * A */
220  template<typename Treal, typename Tmatrix>
224  assert(!sm.tB);
225  Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
226  return *this;
227  }
228 
229 
230  template<typename Treal, typename Tmatrix>
232  thresh(Treal const threshold,
233  normType const norm) {
234  switch (norm) {
235  case frobNorm:
236  this->frob_thresh(threshold);
237  break;
238  default:
239  throw Failure("MatrixTriangular<Treal, Tmatrix>::"
240  "thresh(Treal const, "
241  "normType const): "
242  "Thresholding not imlpemented for selected norm");
243  }
244  }
245 
246  template<typename Treal, typename Tmatrix>
248  eucl(Treal const requestedAccuracy,
249  int maxIter) const {
250  VectorType guess;
251  SizesAndBlocks cols;
252  this->getCols(cols);
253  guess.resetSizesAndBlocks(cols);
254  guess.rand();
256  if (maxIter < 0)
257  maxIter = this->get_nrows() * 100;
259  <Treal, ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal>, VectorType>
260  lan(ztz, guess, maxIter);
261  lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
262  lan.run();
263  Treal eVal = 0;
264  Treal acc = 0;
265  lan.getLargestMagnitudeEig(eVal, acc);
266  Interval<Treal> euclInt( template_blas_sqrt(eVal-acc),
267  template_blas_sqrt(eVal+acc) );
268  if ( euclInt.low() < 0 )
269  euclInt = Interval<Treal>( 0, template_blas_sqrt(eVal+acc) );
270  if ( euclInt.length() / 2.0 > requestedAccuracy ) {
271  std::cout << "req: " << requestedAccuracy
272  << " obt: " << euclInt.length() / 2.0 << std::endl;
273  throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
274  }
275  return euclInt.midPoint();
276  }
277 
278 #if 1
279 
280  template<typename Treal, typename Tmatrix>
282  eucl_thresh(Treal const threshold) {
284  return TruncObj.run( threshold );
285  }
286 
287 #endif
288 
289  template<typename Treal, typename Tmatrix>
291  eucl_thresh_congr_trans_measure(Treal const threshold,
294  return TruncObj.run( threshold );
295  }
296 
297 
298 } /* end namespace mat */
299 #endif
300 
301 
void thresh(Treal const threshold, normType const norm)
Definition: MatrixTriangular.h:232
MatrixTriangular< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &tri)
Definition: MatrixTriangular.h:67
Treal eucl_thresh(Treal const threshold)
Definition: MatrixTriangular.h:282
Normal matrix.
Definition: MatrixBase.h:47
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:164
Treal low() const
Definition: Interval.h:142
inchversion
Definition: Matrix.h:74
void frob_thresh(Treal threshold)
Definition: MatrixTriangular.h:169
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:151
Definition: LanczosLargestMagnitudeEig.h:44
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:279
void inch(const MatrixGeneral< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:141
Truncation of general matrices.
Definition: truncation.h:270
void getCols(SizesAndBlocks &colsCopy) const
Definition: MatrixBase.h:83
Definition: MatrixBase.h:53
Definition: allocate.cc:30
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
MatrixTriangular()
Default constructor.
Definition: MatrixTriangular.h:61
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixTriangular.h:248
Definition: Interval.h:44
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixTriangular.h:59
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixTriangular.h:180
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:234
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:196
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixTriangular.h:208
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:252
size_t nnz() const
Definition: MatrixTriangular.h:172
Upper non-unit triangular matrix.
Definition: MatrixBase.h:51
MatrixTriangular(const MatrixTriangular< Treal, Tmatrix > &tri)
Copy constructor.
Definition: MatrixTriangular.h:63
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixTriangular.h:211
Treal run(Treal const threshold)
Definition: truncation.h:78
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:219
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:83
Treal midPoint() const
Definition: Interval.h:113
void rand()
Definition: VectorGeneral.h:94
void inch(const MatrixSymmetric< Treal, Tmatrix > &SPD, const Treal threshold, const side looking=left, const inchversion version=unstable)
Definition: MatrixTriangular.h:148
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:74
size_t nvalues() const
Definition: MatrixTriangular.h:175
Treal eucl_thresh_congr_trans_measure(Treal const threshold, MatrixSymmetric< Treal, Tmatrix > &trA)
Definition: MatrixTriangular.h:291
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:104
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:49
Treal frob() const
Definition: MatrixTriangular.h:160
Definition: mat_utils.h:69
void random()
Definition: MatrixTriangular.h:187
Base class for matrix API.
Definition: MatrixBase.h:67
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixTriangular.h:183
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:56
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:112
side
Definition: Matrix.h:73
Base class for matrix API.
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
Definition: matInclude.h:135
Definition: Failure.h:47
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:49
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:97
Definition: Matrix.h:73
Definition: MatrixBase.h:54
int get_nrows() const
Definition: MatrixBase.h:136
MatrixTriangular< Treal, Tmatrix > & operator=(int const k)
Set matrix to zero or identity: A = 0 or A = 1.
Definition: MatrixTriangular.h:72
Symmetric matrix.
Definition: MatrixBase.h:49
std::string obj_type_id() const
Definition: MatrixTriangular.h:206
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:125
Truncation of general matrices with impact on matrix triple multiply as error measure.
Definition: truncation.h:336
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:135
Definition: Matrix.h:74