ergo
MatrixTriangular.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 
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>
221  inline MatrixTriangular<Treal, Tmatrix>&
222  MatrixTriangular<Treal, Tmatrix>::operator+=
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;
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