ergo
matrix_utilities.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
47Perm* 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
58void getMatrixPermutation(const BasisInfoStruct& basisInfo,
59 int sparse_block_size,
60 int factor1,
61 int factor2,
62 int factor3,
63 std::vector<int> & permutation);
64void 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);
71void 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
94void output_matrix(int n, const ergo_real* matrix, const char* matrixName);
95
96template<class Tmatrix>
98{
99 return M.maxAbsValue();
100
101}
102
103template<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
111template<typename Tmatrix>
112void 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
148template<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
Code for setting up basis functions starting from shells.
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
#define VERSION
Definition config.h:273
#define A
Header file with typedefs for matrix and vector types.
mat::SizesAndBlocks prepareMatrixSizesAndBlocks(int n_basis_functions, int sparse_block_size, int factor1, int factor2, int factor3)
Definition matrix_utilities.cc:47
void output_matrix(int n, const ergo_real *matrix, const char *matrixName)
Definition matrix_utilities.cc:390
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
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
void add_random_diag_perturbation(int n, symmMatrix &M, ergo_real eps)
Definition matrix_utilities.cc:341
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
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
ergo_real compute_maxabs_sparse(const Tmatrix &M)
Definition matrix_utilities.h:97
void getMatrixPermutation(const BasisInfoStruct &basisInfo, int sparse_block_size, int factor1, int factor2, int factor3, std::vector< int > &permutation)
Definition matrix_utilities.cc:223
void fill_matrix_with_random_numbers(int n, symmMatrix &M)
Definition matrix_utilities.cc:310
double ergo_real
Definition realtype.h:69
Definition basisinfo.h:112
Definition matrix_utilities.h:104
matrix_utilities_CompareClass(RandomAccessIterator firstel)
Definition matrix_utilities.h:106
RandomAccessIterator first
Definition matrix_utilities.h:105
bool operator()(int i, int j)
Definition matrix_utilities.h:108
MatrixSymmetric< real, matri > symmMatrix
Definition test_LanczosSeveralLargestEig.cc:69