ergo
matrix_utilities.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 
28 #ifndef MATRIX_UTILITIES_HEADER
29 #define MATRIX_UTILITIES_HEADER
30 
31 #include "matrix_typedefs.h"
32 #include "basisinfo.h"
33 
34 #if 0
35 
37 Perm* prepare_matrix_permutation(const BasisInfoStruct& basisInfo,
38  int sparse_block_size,
39  int factor1, int factor2, int factor3);
40 #else
41 
43  int sparse_block_size,
44  int factor1,
45  int factor2,
46  int factor3);
47 
48 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
49  int sparse_block_size,
50  int factor1,
51  int factor2,
52  int factor3,
53  std::vector<int> & permutation);
54 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
55  int sparse_block_size,
56  int factor1,
57  int factor2,
58  int factor3,
59  std::vector<int> & permutation,
60  std::vector<int> & inversePermutation);
61 void getMatrixPermutationOnlyFactor2(const std::vector<ergo_real> & xcoords,
62  const std::vector<ergo_real> & ycoords,
63  const std::vector<ergo_real> & zcoords,
64  int sparse_block_size_lowest,
65  int first_factor, // this factor may be different from 2, all other factors are always 2.
66  std::vector<int> & permutation,
67  std::vector<int> & inversePermutation);
69  int sparse_block_size_lowest,
70  int first_factor, // this factor may be different from 2, all other factors are always 2.
71  std::vector<int> & permutation,
72  std::vector<int> & inversePermutation);
73 
74 #endif
76 
78  symmMatrix & M,
79  ergo_real eps);
80 
82  std::vector<int> const & inversePermutationHML);
83 
84 void output_matrix(int n, const ergo_real* matrix, const char* matrixName);
85 
86 template<class Tmatrix>
87 ergo_real compute_maxabs_sparse(const Tmatrix & M)
88 {
89  return M.maxAbsValue();
90 
91 }
92 
93 template<typename RandomAccessIterator>
95  RandomAccessIterator first;
96  explicit matrix_utilities_CompareClass(RandomAccessIterator firstel)
97  : first(firstel){}
98  bool operator() (int i,int j) { return (*(first + i) < *(first + j));}
99 };
100 
101 template<typename Tmatrix>
102 void get_all_nonzeros( Tmatrix const & A,
103  std::vector<int> const & inversePermutation,
104  std::vector<int> & rowind,
105  std::vector<int> & colind,
106  std::vector<ergo_real> & values) {
107  rowind.resize(0);
108  colind.resize(0);
109  values.resize(0);
110  size_t nvalues = 0;
111  size_t nvalues_tmp = A.nvalues();
112  std::vector<int> rowind_tmp; rowind_tmp.reserve(nvalues_tmp);
113  std::vector<int> colind_tmp; colind_tmp.reserve(nvalues_tmp);
114  std::vector<ergo_real> values_tmp; values_tmp.reserve(nvalues_tmp);
115  A.get_all_values(rowind_tmp,
116  colind_tmp,
117  values_tmp,
118  inversePermutation,
119  inversePermutation);
120  // Count the number of nonzeros
121  for(size_t i = 0; i < nvalues_tmp; i++) {
122  nvalues += ( values_tmp[i] != 0 );
123  }
124  rowind.reserve(nvalues);
125  colind.reserve(nvalues);
126  values.reserve(nvalues);
127  // Extract all nonzeros
128  for(size_t i = 0; i < nvalues_tmp; i++) {
129  if ( values_tmp[i] != 0 ) {
130  rowind.push_back( rowind_tmp[i] );
131  colind.push_back( colind_tmp[i] );
132  values.push_back( values_tmp[i] );
133  }
134  }
135 } // end get_all_nonzeros(...)
136 
137 
138 
139 template<typename Tmatrix>
141  Tmatrix const & A,
142  std::vector<int> const & inversePermutation,
143  std::string name,
144  int resolution_r,
145  int resolution_m
146  ) {
147  std::string m_name = name + ".m";
148  std::ofstream os(m_name.c_str());
149  // Get xyz coords for all indices
150  int n = basisInfo.noOfBasisFuncs;
151  std::vector<ergo_real> x(n);
152  std::vector<ergo_real> y(n);
153  std::vector<ergo_real> z(n);
154  for(int i = 0; i < n; i++) {
155  x[i] = basisInfo.basisFuncList[i].centerCoords[0];
156  y[i] = basisInfo.basisFuncList[i].centerCoords[1];
157  z[i] = basisInfo.basisFuncList[i].centerCoords[2];
158  }
159 
160  size_t number_of_stored_zeros = 0;
161  ergo_real minAbsValue = 1e22;
162  ergo_real maxAbsValue = 0;
163 
164  // Get all matrix elements
165  size_t nvalues = 0;
166  std::vector<int> rowind;
167  std::vector<int> colind;
168  std::vector<ergo_real> values;
169  {
170  std::vector<int> rowind_tmp;
171  std::vector<int> colind_tmp;
172  std::vector<ergo_real> values_tmp;
173  get_all_nonzeros( A, inversePermutation, rowind_tmp, colind_tmp, values_tmp);
174 
175  bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
176  if (matrixIsSymmetric) {
177  // Also include lower triangle
178  size_t nvalues_tmp = values_tmp.size();
179  rowind.reserve(nvalues_tmp*2);
180  colind.reserve(nvalues_tmp*2);
181  values.reserve(nvalues_tmp*2);
182  for(size_t i = 0; i < nvalues_tmp; i++) {
183  rowind.push_back( rowind_tmp[i] );
184  colind.push_back( colind_tmp[i] );
185  values.push_back( values_tmp[i] );
186  if ( rowind_tmp[i] != colind_tmp[i] ) {
187  rowind.push_back( colind_tmp[i] );
188  colind.push_back( rowind_tmp[i] );
189  values.push_back( values_tmp[i] );
190  }
191  } // end for
192  } // end if
193  else {
194  rowind = rowind_tmp;
195  colind = colind_tmp;
196  values = values_tmp;
197  } // end else
198 
199  nvalues = values.size();
200  // Take absolute value
201  for(size_t i = 0; i < nvalues; i++) {
202  ergo_real fabsVal = fabs( values[i] );
203  values[i] = fabsVal;
204  minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
205  maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
206  }
207  }
208 
209  os << "%% Run for example like this: matlab -nosplash -nodesktop -r " << name << std::endl;
210  os << "number_of_stored_zeros = " << number_of_stored_zeros << ";" << std::endl;
211  os << "number_of_stored_nonzeros = " << nvalues << ";" << std::endl;
212 
213  // Get distances for all matrix elements
214  std::vector<ergo_real> distances(nvalues);
215  for(size_t i = 0; i < nvalues; i++) {
216  ergo_real diff_x = x[ rowind[i] ] - x[ colind[i] ];
217  ergo_real diff_y = y[ rowind[i] ] - y[ colind[i] ];
218  ergo_real diff_z = z[ rowind[i] ] - z[ colind[i] ];
219  distances[i] = std::sqrt( diff_x * diff_x + diff_y * diff_y + diff_z * diff_z );
220  }
221 
222  // Index vector
223  std::vector<size_t> index(nvalues);
224  for ( size_t ind = 0; ind < index.size(); ++ind ) {
225  index[ind] = ind;
226  }
227 
228  // Sort based on distance
230  compareDist( distances.begin() );
231  std::sort ( index.begin(), index.end(), compareDist );
232 
233  // Min and max distances
234  ergo_real minDistance = *std::min_element( distances.begin(), distances.end() );
235  ergo_real maxDistance = *std::max_element( distances.begin(), distances.end() );
236  // Size of box in r direction
237  ergo_real rbox_length = (maxDistance - minDistance) / resolution_r;
238 
239  // Get max absolute value of A
240  ergo_real maxMagLog10 = std::log10(maxAbsValue);
241  ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
242  // Size of box in m direction
243  ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
244 
245  os << "A = [ " << std::endl;
246  // Loop over r boxes
247  size_t start_ind = 0;
248  ergo_real r_low = minDistance;
249  for ( int rbox = 0; rbox < resolution_r; rbox++ ) {
250  ergo_real r_upp = r_low + rbox_length;
251  // Find end index
252  size_t end_ind = start_ind;
253  while ( end_ind < nvalues && distances[index[end_ind]] < r_upp )
254  end_ind++;
255  // Now we have the bounds for box in r-direction
256  // Sort based on magnitude
258  compareMagnitude( values.begin() );
259  std::sort ( index.begin() + start_ind, index.begin() + end_ind, compareMagnitude );
260  // Loop over m boxes
261  ergo_real m_low = minMagLog10;
262  size_t ind_m = start_ind;
263 
264  // Skip very small values
265  while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_low )
266  ind_m++;
267  size_t skipped_small = ind_m - start_ind;
268  os << r_low << " "
269  << r_upp << " "
270  << 0 << " "
271  << std::pow(10,m_low) << " "
272  << skipped_small
273  << std::endl;
274 
275  for ( int mbox = 0; mbox < resolution_m; mbox++ ) {
276  ergo_real m_upp = m_low + mbox_length;
277  size_t count = 0;
278  while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_upp ) {
279  ind_m++;
280  count++;
281  }
282  // Now we have r_low r_upp m_low m_upp count
283  // Write to stream
284  os << r_low << " "
285  << r_upp << " "
286  << std::pow(10,m_low) << " "
287  << std::pow(10,m_upp) << " "
288  << count
289  << std::endl;
290  m_low = m_upp;
291  }
292 
293  r_low = r_upp;
294  start_ind = end_ind;
295  }
296  os << "];" << std::endl;
297  os << "B=[];" << std::endl;
298  os << "for ind = 1 : size(A,1)" << std::endl;
299  os << " if (A(ind,3) ~= 0)" << std::endl;
300  os << " B = [B; A(ind,:)];" << std::endl;
301  os << " end" << std::endl;
302  os << "end" << std::endl;
303  os << "%col = jet(101);" << std::endl;
304  os << "col = gray(101);col=col(end:-1:1,:);" << std::endl;
305  os << "maxCount = max(B(:,5));" << std::endl;
306  os << "ax = [0 30 1e-12 1e3]" << std::endl;
307 
308  os << "fighandle = figure;" << std::endl;
309  os << "for ind = 1 : size(B,1)" << std::endl;
310  os << " rmin = B(ind, 1); rmax = B(ind, 2);" << std::endl;
311  os << " mmin = B(ind, 3); mmax = B(ind, 4);" << std::endl;
312  os << " colind = round( 1+100 * B(ind,5) / maxCount);" << std::endl;
313  os << " fill([rmin rmin rmax rmax rmin], [mmin mmax mmax mmin mmin], col(colind,:), 'EdgeColor', col(colind,:) )" << std::endl;
314  os << " hold on" << std::endl;
315  os << "end" << std::endl;
316  os << "set(gca,'YScale','log')" << std::endl;
317  os << "axis(ax)" << std::endl;
318  os << "set(gca,'FontSize',16)" << std::endl;
319  os << "xlabel('Distance')" << std::endl;
320  os << "ylabel('Magnitude')" << std::endl;
321  os << "print( fighandle, '-depsc2', '" << name << "')" << std::endl;
322 
323  os << "fighandle = figure;" << std::endl;
324  os << "for ind = 1 : size(B,1)" << std::endl;
325  os << " if (B(ind,5) ~= 0)" << std::endl;
326  os << " rmin = B(ind, 1); rmax = B(ind, 2);" << std::endl;
327  os << " mmin = B(ind, 3); mmax = B(ind, 4);" << std::endl;
328  os << " msize = 3+1*ceil(20 * B(ind,5) / maxCount);" << std::endl;
329  os << " plot((rmin+rmax)/2,(mmin+mmax)/2,'k.','MarkerSize',msize)" << std::endl;
330  os << " hold on" << std::endl;
331  os << " end" << std::endl;
332  os << "end" << std::endl;
333  os << "set(gca,'YScale','log')" << std::endl;
334  os << "axis(ax)" << std::endl;
335  os << "set(gca,'FontSize',16)" << std::endl;
336  os << "xlabel('Distance')" << std::endl;
337  os << "ylabel('Magnitude')" << std::endl;
338  os << "print( fighandle, '-depsc2', '" << name << "_dots')" << std::endl;
339  os << "exit(0);" << std::endl;
340  os.close();
341 } // end output_distance_vs_magnitude(...)
342 
343 template<typename Tmatrix>
344 void output_magnitude_histogram( Tmatrix const & A,
345  std::string name,
346  int resolution_m
347  ) {
348  std::string m_name = name + ".m";
349  std::ofstream os(m_name.c_str());
350 
351  size_t number_of_stored_zeros = 0;
352  ergo_real minAbsValue = 1e22;
353  ergo_real maxAbsValue = 0;
354 
355  // Get all matrix elements
356  size_t nvalues = 0;
357  std::vector<int> rowind;
358  std::vector<int> colind;
359  std::vector<ergo_real> values;
360  {
361  // Get all nonzeros
362  rowind.resize(0);
363  colind.resize(0);
364  values.resize(0);
365  size_t nvalues_tmp = A.nvalues();
366  std::vector<int> rowind_tmp; rowind_tmp.reserve(nvalues_tmp);
367  std::vector<int> colind_tmp; colind_tmp.reserve(nvalues_tmp);
368  std::vector<ergo_real> values_tmp; values_tmp.reserve(nvalues_tmp);
369  A.get_all_values(rowind_tmp,
370  colind_tmp,
371  values_tmp);
372  // Count the number of nonzeros
373  for(size_t i = 0; i < nvalues_tmp; i++) {
374  nvalues += ( values_tmp[i] != 0 );
375  }
376 
377  bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
378  if (matrixIsSymmetric) {
379  // Also include lower triangle
380  rowind.reserve(nvalues*2);
381  colind.reserve(nvalues*2);
382  values.reserve(nvalues*2);
383  // Extract all nonzeros
384  for(size_t i = 0; i < nvalues_tmp; i++) {
385  if ( values_tmp[i] != 0 ) {
386  rowind.push_back( rowind_tmp[i] );
387  colind.push_back( colind_tmp[i] );
388  values.push_back( values_tmp[i] );
389  if ( rowind_tmp[i] != colind_tmp[i] ) {
390  rowind.push_back( colind_tmp[i] );
391  colind.push_back( rowind_tmp[i] );
392  values.push_back( values_tmp[i] );
393  }
394  }
395  }
396  nvalues = values.size();
397  } // end if
398  else {
399  rowind.reserve(nvalues);
400  colind.reserve(nvalues);
401  values.reserve(nvalues);
402  // Extract all nonzeros
403  for(size_t i = 0; i < nvalues_tmp; i++) {
404  if ( values_tmp[i] != 0 ) {
405  rowind.push_back( rowind_tmp[i] );
406  colind.push_back( colind_tmp[i] );
407  values.push_back( values_tmp[i] );
408  }
409  }
410  assert( nvalues == values.size() );
411  } // end else
412  // Take absolute value
413  for(size_t i = 0; i < nvalues; i++) {
414  ergo_real fabsVal = fabs( values[i] );
415  values[i] = fabsVal;
416  minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
417  maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
418  }
419  }
420 
421  os << "%% Run for example like this: matlab -nosplash -nodesktop -r " << name << std::endl;
422  os << "number_of_stored_zeros = " << number_of_stored_zeros << ";" << std::endl;
423  os << "number_of_stored_nonzeros = " << nvalues << ";" << std::endl;
424 
425  // Index vector
426  std::vector<size_t> index(nvalues);
427  for ( size_t ind = 0; ind < index.size(); ++ind ) {
428  index[ind] = ind;
429  }
430 
431  // Get min and max absolute value of A
432  ergo_real maxMagLog10 = std::log10(maxAbsValue);
433  ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
434  // Size of box in m direction
435  ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
436 
437  os << "A = [ " << std::endl;
438  // Sort based on magnitude
440  compareMagnitude( values.begin() );
441  std::sort ( index.begin(), index.end(), compareMagnitude );
442  // Loop over m boxes
443  ergo_real m_low = minMagLog10;
444  size_t ind_m = 0;
445 
446  // Skip very small values
447  while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_low )
448  ind_m++;
449  size_t skipped_small = ind_m;
450  os << 0 << " "
451  << std::pow(10,m_low) << " "
452  << skipped_small
453  << std::endl;
454 
455  for ( int mbox = 0; mbox < resolution_m; mbox++ ) {
456  ergo_real m_upp = m_low + mbox_length;
457  size_t count = 0;
458  while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_upp ) {
459  ind_m++;
460  count++;
461  }
462  // Now we have m_low m_upp count
463  // Write to stream
464  os << std::pow(10,m_low) << " "
465  << std::pow(10,m_upp) << " "
466  << count
467  << std::endl;
468  m_low = m_upp;
469  }
470  os << "];" << std::endl;
471 
472  os.close();
473 } // end output_magnitude_histogram
474 
475 template<typename Tmatrix>
477  std::vector<int> const & inversePermutation,
478  std::string filename,
479  std::string identifier,
480  std::string method_and_basis)
481 {
482 
483  // Get all matrix elements
484  size_t nvalues = 0;
485  std::vector<int> rowind;
486  std::vector<int> colind;
487  std::vector<ergo_real> values;
488  get_all_nonzeros( A, inversePermutation, rowind, colind, values);
489  nvalues = values.size();
490  // Now we have all matrix elements
491  // Open file stream
492  std::string mtx_filename = filename + ".mtx";
493  std::ofstream os(mtx_filename.c_str());
494 
495  time_t rawtime;
496  struct tm * timeinfo;
497  time ( &rawtime );
498  timeinfo = localtime ( &rawtime );
499 
500  std::string matrix_market_matrix_type = "general";
501  bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
502  if (matrixIsSymmetric)
503  matrix_market_matrix_type = "symmetric";
504  os << "%%MatrixMarket matrix coordinate real " << matrix_market_matrix_type << std::endl
505  << "%===============================================================================" << std::endl
506  << "% Generated by the Ergo quantum chemistry program version " << VERSION << " (www.ergoscf.org)" << std::endl
507  << "% Date : " << asctime (timeinfo) // newline added by asctime
508  << "% ID-string : " << identifier << std::endl
509  << "% Method : " << method_and_basis << std::endl
510  << "%" << std::endl
511  << "% MatrixMarket file format:" << std::endl
512  << "% +-----------------" << std::endl
513  << "% | % comments" << std::endl
514  << "% | nrows ncols nentries" << std::endl
515  << "% | i_1 j_1 A(i_1,j_1)" << std::endl
516  << "% | i_2 j_2 A(i_1,j_1)" << std::endl
517  << "% | ..." << std::endl
518  << "% | i_nentries j_nentries A(i_nentries,j_nentries) " << std::endl
519  << "% +----------------" << std::endl
520  << "% Note that indices are 1-based, i.e. A(1,1) is the first element." << std::endl
521  << "%" << std::endl
522  << "%===============================================================================" << std::endl;
523  os << A.get_nrows() << " " << A.get_ncols() << " " << nvalues << std::endl;
524  if (matrixIsSymmetric)
525  for(size_t i = 0; i < nvalues; i++) {
526  // Output lower triangle
527  if ( rowind[i] < colind[i] )
528  os << colind[i]+1 << " " << rowind[i]+1 << " " << std::setprecision(10) << values[i] << std::endl;
529  else
530  os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << values[i] << std::endl;
531  }
532  else
533  for(size_t i = 0; i < nvalues; i++) {
534  os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << values[i] << std::endl;
535  }
536  os.close();
537 } // end write_matrix_in_matrix_market_format(...)
538 
539 
540 #endif
#define A
void fill_matrix_with_random_numbers(int n, symmMatrix &M)
Definition: matrix_utilities.cc:275
double ergo_real
Definition: realtype.h:53
bool check_if_matrix_contains_strange_elements(const symmMatrix &M, std::vector< int > const &inversePermutationHML)
This function is supposed to check if a matrix contains any strange numbers such as "inf" or "nan"...
Definition: matrix_utilities.cc:327
void output_magnitude_histogram(Tmatrix const &A, std::string name, int resolution_m)
Definition: matrix_utilities.h:344
void output_distance_vs_magnitude(BasisInfoStruct const &basisInfo, Tmatrix const &A, std::vector< int > const &inversePermutation, std::string name, int resolution_r, int resolution_m)
Definition: matrix_utilities.h:140
#define VERSION
Definition: config.h:202
mat::SizesAndBlocks prepareMatrixSizesAndBlocks(int n_basis_functions, int sparse_block_size, int factor1, int factor2, int factor3)
Definition: matrix_utilities.cc:37
BasisFuncStruct * basisFuncList
Definition: basisinfo.h:120
void add_random_diag_perturbation(int n, symmMatrix &M, ergo_real eps)
Definition: matrix_utilities.cc:306
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
Vector3D centerCoords
Definition: basisinfo.h:86
ergo_real compute_maxabs_sparse(const Tmatrix &M)
Definition: matrix_utilities.h:87
void get_all_nonzeros(Tmatrix const &A, std::vector< int > const &inversePermutation, std::vector< int > &rowind, std::vector< int > &colind, std::vector< ergo_real > &values)
Definition: matrix_utilities.h:102
int noOfBasisFuncs
Definition: basisinfo.h:119
Definition: matrix_utilities.h:94
void output_matrix(int n, const ergo_real *matrix, const char *matrixName)
Definition: matrix_utilities.cc:355
Definition: basisinfo.h:111
void getMatrixPermutationOnlyFactor2(const std::vector< ergo_real > &xcoords, const std::vector< ergo_real > &ycoords, const std::vector< ergo_real > &zcoords, int sparse_block_size_lowest, int first_factor, std::vector< int > &permutation, std::vector< int > &inversePermutation)
Definition: matrix_utilities.cc:205
bool operator()(int i, int j)
Definition: matrix_utilities.h:98
Header file with typedefs for matrix and vector types.
void write_matrix_in_matrix_market_format(Tmatrix const &A, std::vector< int > const &inversePermutation, std::string filename, std::string identifier, std::string method_and_basis)
Definition: matrix_utilities.h:476
matrix_utilities_CompareClass(RandomAccessIterator firstel)
Definition: matrix_utilities.h:96
void getMatrixPermutation(const BasisInfoStruct &basisInfo, int sparse_block_size, int factor1, int factor2, int factor3, std::vector< int > &permutation)
Definition: matrix_utilities.cc:188
RandomAccessIterator first
Definition: matrix_utilities.h:95
Symmetric matrix.
Definition: MatrixBase.h:49