ergo
matrix_utilities.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 MATRIX_UTILITIES_HEADER
39 #define MATRIX_UTILITIES_HEADER
40 
41 #include "matrix_typedefs.h"
42 #include "basisinfo.h"
43 
44 #if 0
45 
47 Perm* prepare_matrix_permutation(const BasisInfoStruct& basisInfo,
48  int sparse_block_size,
49  int factor1, int factor2, int factor3);
50 #else
51 
53  int sparse_block_size,
54  int factor1,
55  int factor2,
56  int factor3);
57 
58 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
59  int sparse_block_size,
60  int factor1,
61  int factor2,
62  int factor3,
63  std::vector<int> & permutation);
64 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
65  int sparse_block_size,
66  int factor1,
67  int factor2,
68  int factor3,
69  std::vector<int> & permutation,
70  std::vector<int> & inversePermutation);
71 void getMatrixPermutationOnlyFactor2(const std::vector<ergo_real> & xcoords,
72  const std::vector<ergo_real> & ycoords,
73  const std::vector<ergo_real> & zcoords,
74  int sparse_block_size_lowest,
75  int first_factor, // this factor may be different from 2, all other factors are always 2.
76  std::vector<int> & permutation,
77  std::vector<int> & inversePermutation);
79  int sparse_block_size_lowest,
80  int first_factor, // this factor may be different from 2, all other factors are always 2.
81  std::vector<int> & permutation,
82  std::vector<int> & inversePermutation);
83 
84 #endif
86 
88  symmMatrix & M,
89  ergo_real eps);
90 
92  std::vector<int> const & inversePermutationHML);
93 
94 void output_matrix(int n, const ergo_real* matrix, const char* matrixName);
95 
96 template<class Tmatrix>
97 ergo_real compute_maxabs_sparse(const Tmatrix & M)
98 {
99  return M.maxAbsValue();
100 
101 }
102 
103 template<typename RandomAccessIterator>
105  RandomAccessIterator first;
106  explicit matrix_utilities_CompareClass(RandomAccessIterator firstel)
107  : first(firstel){}
108  bool operator() (int i,int j) { return (*(first + i) < *(first + j));}
109 };
110 
111 template<typename Tmatrix>
112 void get_all_nonzeros( Tmatrix const & A,
113  std::vector<int> const & inversePermutation,
114  std::vector<int> & rowind,
115  std::vector<int> & colind,
116  std::vector<ergo_real> & values) {
117  rowind.resize(0);
118  colind.resize(0);
119  values.resize(0);
120  size_t nvalues = 0;
121  size_t nvalues_tmp = A.nvalues();
122  std::vector<int> rowind_tmp; rowind_tmp.reserve(nvalues_tmp);
123  std::vector<int> colind_tmp; colind_tmp.reserve(nvalues_tmp);
124  std::vector<ergo_real> values_tmp; values_tmp.reserve(nvalues_tmp);
125  A.get_all_values(rowind_tmp,
126  colind_tmp,
127  values_tmp,
128  inversePermutation,
129  inversePermutation);
130  // Count the number of nonzeros
131  for(size_t i = 0; i < nvalues_tmp; i++) {
132  nvalues += ( values_tmp[i] != 0 );
133  }
134  rowind.reserve(nvalues);
135  colind.reserve(nvalues);
136  values.reserve(nvalues);
137  // Extract all nonzeros
138  for(size_t i = 0; i < nvalues_tmp; i++) {
139  if ( values_tmp[i] != 0 ) {
140  rowind.push_back( rowind_tmp[i] );
141  colind.push_back( colind_tmp[i] );
142  values.push_back( values_tmp[i] );
143  }
144  }
145 } // end get_all_nonzeros(...)
146 
147 
148 template<typename Tmatrix>
150  std::vector<int> const & inversePermutation,
151  std::string filename,
152  std::string identifier,
153  std::string method_and_basis)
154 {
155 
156  // Get all matrix elements
157  size_t nvalues = 0;
158  std::vector<int> rowind;
159  std::vector<int> colind;
160  std::vector<ergo_real> values;
161  get_all_nonzeros( A, inversePermutation, rowind, colind, values);
162  nvalues = values.size();
163  // Now we have all matrix elements
164  // Open file stream
165  std::string mtx_filename = filename + ".mtx";
166  std::ofstream os(mtx_filename.c_str());
167 
168  time_t rawtime;
169  struct tm * timeinfo;
170  time ( &rawtime );
171  timeinfo = localtime ( &rawtime );
172 
173  std::string matrix_market_matrix_type = "general";
174  bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
175  if (matrixIsSymmetric)
176  matrix_market_matrix_type = "symmetric";
177  os << "%%MatrixMarket matrix coordinate real " << matrix_market_matrix_type << std::endl
178  << "%===============================================================================" << std::endl
179  << "% Generated by the Ergo quantum chemistry program version " << VERSION << " (www.ergoscf.org)" << std::endl
180  << "% Date : " << asctime (timeinfo) // newline added by asctime
181  << "% ID-string : " << identifier << std::endl
182  << "% Method : " << method_and_basis << std::endl
183  << "%" << std::endl
184  << "% MatrixMarket file format:" << std::endl
185  << "% +-----------------" << std::endl
186  << "% | % comments" << std::endl
187  << "% | nrows ncols nentries" << std::endl
188  << "% | i_1 j_1 A(i_1,j_1)" << std::endl
189  << "% | i_2 j_2 A(i_2,j_2)" << std::endl
190  << "% | ..." << std::endl
191  << "% | i_nentries j_nentries A(i_nentries,j_nentries) " << std::endl
192  << "% +----------------" << std::endl
193  << "% Note that indices are 1-based, i.e. A(1,1) is the first element." << std::endl
194  << "%" << std::endl
195  << "%===============================================================================" << std::endl;
196  os << A.get_nrows() << " " << A.get_ncols() << " " << nvalues << std::endl;
197  if (matrixIsSymmetric)
198  for(size_t i = 0; i < nvalues; i++) {
199  // Output lower triangle
200  if ( rowind[i] < colind[i] )
201  os << colind[i]+1 << " " << rowind[i]+1 << " " << std::setprecision(10) << (double)values[i] << std::endl;
202  else
203  os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << (double)values[i] << std::endl;
204  }
205  else
206  for(size_t i = 0; i < nvalues; i++) {
207  os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << (double)values[i] << std::endl;
208  }
209  os.close();
210 } // end write_matrix_in_matrix_market_format(...)
211 
212 
213 #endif
#define A
void fill_matrix_with_random_numbers(int n, symmMatrix &M)
Definition: matrix_utilities.cc:310
double ergo_real
Definition: realtype.h:69
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:362
#define VERSION
Definition: config.h:268
mat::SizesAndBlocks prepareMatrixSizesAndBlocks(int n_basis_functions, int sparse_block_size, int factor1, int factor2, int factor3)
Definition: matrix_utilities.cc:47
MatrixSymmetric< real, matri > symmMatrix
Definition: test_LanczosSeveralLargestEig.cc:69
void add_random_diag_perturbation(int n, symmMatrix &M, ergo_real eps)
Definition: matrix_utilities.cc:341
Code for setting up basis functions starting from shells.
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
ergo_real compute_maxabs_sparse(const Tmatrix &M)
Definition: matrix_utilities.h:97
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:112
Definition: matrix_utilities.h:104
void output_matrix(int n, const ergo_real *matrix, const char *matrixName)
Definition: matrix_utilities.cc:390
Definition: basisinfo.h:112
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:240
bool operator()(int i, int j)
Definition: matrix_utilities.h:108
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:149
matrix_utilities_CompareClass(RandomAccessIterator firstel)
Definition: matrix_utilities.h:106
void getMatrixPermutation(const BasisInfoStruct &basisInfo, int sparse_block_size, int factor1, int factor2, int factor3, std::vector< int > &permutation)
Definition: matrix_utilities.cc:223
RandomAccessIterator first
Definition: matrix_utilities.h:105