28 #ifndef MATRIX_UTILITIES_HEADER
29 #define MATRIX_UTILITIES_HEADER
38 int sparse_block_size,
39 int factor1,
int factor2,
int factor3);
43 int sparse_block_size,
49 int sparse_block_size,
53 std::vector<int> & permutation);
55 int sparse_block_size,
59 std::vector<int> & permutation,
60 std::vector<int> & inversePermutation);
62 const std::vector<ergo_real> & ycoords,
63 const std::vector<ergo_real> & zcoords,
64 int sparse_block_size_lowest,
66 std::vector<int> & permutation,
67 std::vector<int> & inversePermutation);
69 int sparse_block_size_lowest,
71 std::vector<int> & permutation,
72 std::vector<int> & inversePermutation);
82 std::vector<int>
const & inversePermutationHML);
86 template<
class Tmatrix>
89 return M.maxAbsValue();
93 template<
typename RandomAccessIterator>
101 template<
typename Tmatrix>
103 std::vector<int>
const & inversePermutation,
104 std::vector<int> & rowind,
105 std::vector<int> & colind,
106 std::vector<ergo_real> & values) {
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,
121 for(
size_t i = 0; i < nvalues_tmp; i++) {
122 nvalues += ( values_tmp[i] != 0 );
124 rowind.reserve(nvalues);
125 colind.reserve(nvalues);
126 values.reserve(nvalues);
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] );
139 template<
typename Tmatrix>
142 std::vector<int>
const & inversePermutation,
147 std::string m_name = name +
".m";
148 std::ofstream os(m_name.c_str());
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++) {
160 size_t number_of_stored_zeros = 0;
166 std::vector<int> rowind;
167 std::vector<int> colind;
168 std::vector<ergo_real> values;
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);
175 bool matrixIsSymmetric = (A.obj_type_id() ==
"MatrixSymmetric");
176 if (matrixIsSymmetric) {
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] );
199 nvalues = values.size();
201 for(
size_t i = 0; i < nvalues; i++) {
204 minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
205 maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
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;
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 );
223 std::vector<size_t> index(nvalues);
224 for (
size_t ind = 0; ind < index.size(); ++ind ) {
230 compareDist( distances.begin() );
231 std::sort ( index.begin(), index.end(), compareDist );
234 ergo_real minDistance = *std::min_element( distances.begin(), distances.end() );
235 ergo_real maxDistance = *std::max_element( distances.begin(), distances.end() );
237 ergo_real rbox_length = (maxDistance - minDistance) / resolution_r;
240 ergo_real maxMagLog10 = std::log10(maxAbsValue);
241 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
243 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
245 os <<
"A = [ " << std::endl;
247 size_t start_ind = 0;
249 for (
int rbox = 0; rbox < resolution_r; rbox++ ) {
252 size_t end_ind = start_ind;
253 while ( end_ind < nvalues && distances[index[end_ind]] < r_upp )
258 compareMagnitude( values.begin() );
259 std::sort ( index.begin() + start_ind, index.begin() + end_ind, compareMagnitude );
262 size_t ind_m = start_ind;
265 while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_low )
267 size_t skipped_small = ind_m - start_ind;
271 << std::pow(10,m_low) <<
" "
275 for (
int mbox = 0; mbox < resolution_m; mbox++ ) {
278 while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_upp ) {
286 << std::pow(10,m_low) <<
" "
287 << std::pow(10,m_upp) <<
" "
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;
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;
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;
343 template<
typename Tmatrix>
348 std::string m_name = name +
".m";
349 std::ofstream os(m_name.c_str());
351 size_t number_of_stored_zeros = 0;
357 std::vector<int> rowind;
358 std::vector<int> colind;
359 std::vector<ergo_real> values;
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,
373 for(
size_t i = 0; i < nvalues_tmp; i++) {
374 nvalues += ( values_tmp[i] != 0 );
377 bool matrixIsSymmetric = (A.obj_type_id() ==
"MatrixSymmetric");
378 if (matrixIsSymmetric) {
380 rowind.reserve(nvalues*2);
381 colind.reserve(nvalues*2);
382 values.reserve(nvalues*2);
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] );
396 nvalues = values.size();
399 rowind.reserve(nvalues);
400 colind.reserve(nvalues);
401 values.reserve(nvalues);
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] );
410 assert( nvalues == values.size() );
413 for(
size_t i = 0; i < nvalues; i++) {
416 minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
417 maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
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;
426 std::vector<size_t> index(nvalues);
427 for (
size_t ind = 0; ind < index.size(); ++ind ) {
432 ergo_real maxMagLog10 = std::log10(maxAbsValue);
433 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
435 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
437 os <<
"A = [ " << std::endl;
440 compareMagnitude( values.begin() );
441 std::sort ( index.begin(), index.end(), compareMagnitude );
447 while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_low )
449 size_t skipped_small = ind_m;
451 << std::pow(10,m_low) <<
" "
455 for (
int mbox = 0; mbox < resolution_m; mbox++ ) {
458 while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_upp ) {
464 os << std::pow(10,m_low) <<
" "
465 << std::pow(10,m_upp) <<
" "
470 os <<
"];" << std::endl;
475 template<
typename Tmatrix>
477 std::vector<int>
const & inversePermutation,
478 std::string filename,
479 std::string identifier,
480 std::string method_and_basis)
485 std::vector<int> rowind;
486 std::vector<int> colind;
487 std::vector<ergo_real> values;
489 nvalues = values.size();
492 std::string mtx_filename = filename +
".mtx";
493 std::ofstream os(mtx_filename.c_str());
496 struct tm * timeinfo;
498 timeinfo = localtime ( &rawtime );
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)
508 <<
"% ID-string : " << identifier << std::endl
509 <<
"% Method : " << method_and_basis << 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
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++) {
527 if ( rowind[i] < colind[i] )
528 os << colind[i]+1 <<
" " << rowind[i]+1 <<
" " << std::setprecision(10) << values[i] << std::endl;
530 os << rowind[i]+1 <<
" " << colind[i]+1 <<
" " << std::setprecision(10) << values[i] << std::endl;
533 for(
size_t i = 0; i < nvalues; i++) {
534 os << rowind[i]+1 <<
" " << colind[i]+1 <<
" " << std::setprecision(10) << values[i] << std::endl;
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