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