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);
71 template<
class Tmatrix>
74 return M.maxAbsValue();
78 template<
typename RandomAccessIterator>
86 template<
typename Tmatrix>
88 std::vector<int>
const & inversePermutation,
89 std::vector<int> & rowind,
90 std::vector<int> & colind,
91 std::vector<ergo_real> & values) {
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,
106 for(
size_t i = 0; i < nvalues_tmp; i++) {
107 nvalues += ( values_tmp[i] != 0 );
109 rowind.reserve(nvalues);
110 colind.reserve(nvalues);
111 values.reserve(nvalues);
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] );
124 template<
typename Tmatrix>
127 std::vector<int>
const & inversePermutation,
132 std::string m_name = name +
".m";
133 std::ofstream os(m_name.c_str());
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++) {
145 size_t number_of_stored_zeros = 0;
151 std::vector<int> rowind;
152 std::vector<int> colind;
153 std::vector<ergo_real> values;
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);
160 bool matrixIsSymmetric = (A.obj_type_id() ==
"MatrixSymmetric");
161 if (matrixIsSymmetric) {
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] );
184 nvalues = values.size();
186 for(
size_t i = 0; i < nvalues; i++) {
189 minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
190 maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
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;
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 );
208 std::vector<size_t> index(nvalues);
209 for (
size_t ind = 0; ind < index.size(); ++ind ) {
215 compareDist( distances.begin() );
216 std::sort ( index.begin(), index.end(), compareDist );
219 ergo_real minDistance = *std::min_element( distances.begin(), distances.end() );
220 ergo_real maxDistance = *std::max_element( distances.begin(), distances.end() );
222 ergo_real rbox_length = (maxDistance - minDistance) / resolution_r;
225 ergo_real maxMagLog10 = std::log10(maxAbsValue);
226 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
228 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
230 os <<
"A = [ " << std::endl;
232 size_t start_ind = 0;
234 for (
int rbox = 0; rbox < resolution_r; rbox++ ) {
237 size_t end_ind = start_ind;
238 while ( end_ind < nvalues && distances[index[end_ind]] < r_upp )
243 compareMagnitude( values.begin() );
244 std::sort ( index.begin() + start_ind, index.begin() + end_ind, compareMagnitude );
247 size_t ind_m = start_ind;
250 while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_low )
252 size_t skipped_small = ind_m - start_ind;
256 << std::pow(10,m_low) <<
" "
260 for (
int mbox = 0; mbox < resolution_m; mbox++ ) {
263 while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_upp ) {
271 << std::pow(10,m_low) <<
" "
272 << std::pow(10,m_upp) <<
" "
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;
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;
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;
328 template<
typename Tmatrix>
333 std::string m_name = name +
".m";
334 std::ofstream os(m_name.c_str());
336 size_t number_of_stored_zeros = 0;
342 std::vector<int> rowind;
343 std::vector<int> colind;
344 std::vector<ergo_real> values;
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,
358 for(
size_t i = 0; i < nvalues_tmp; i++) {
359 nvalues += ( values_tmp[i] != 0 );
362 bool matrixIsSymmetric = (A.obj_type_id() ==
"MatrixSymmetric");
363 if (matrixIsSymmetric) {
365 rowind.reserve(nvalues*2);
366 colind.reserve(nvalues*2);
367 values.reserve(nvalues*2);
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] );
381 nvalues = values.size();
384 rowind.reserve(nvalues);
385 colind.reserve(nvalues);
386 values.reserve(nvalues);
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] );
395 assert( nvalues == values.size() );
398 for(
size_t i = 0; i < nvalues; i++) {
401 minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
402 maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
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;
411 std::vector<size_t> index(nvalues);
412 for (
size_t ind = 0; ind < index.size(); ++ind ) {
417 ergo_real maxMagLog10 = std::log10(maxAbsValue);
418 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
420 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
422 os <<
"A = [ " << std::endl;
425 compareMagnitude( values.begin() );
426 std::sort ( index.begin(), index.end(), compareMagnitude );
432 while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_low )
434 size_t skipped_small = ind_m;
436 << std::pow(10,m_low) <<
" "
440 for (
int mbox = 0; mbox < resolution_m; mbox++ ) {
443 while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_upp ) {
449 os << std::pow(10,m_low) <<
" "
450 << std::pow(10,m_upp) <<
" "
455 os <<
"];" << std::endl;
460 template<
typename Tmatrix>
462 std::vector<int>
const & inversePermutation,
463 std::string filename,
464 std::string identifier,
465 std::string method_and_basis)
470 std::vector<int> rowind;
471 std::vector<int> colind;
472 std::vector<ergo_real> values;
474 nvalues = values.size();
477 std::string mtx_filename = filename +
".mtx";
478 std::ofstream os(mtx_filename.c_str());
481 struct tm * timeinfo;
483 timeinfo = localtime ( &rawtime );
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)
493 <<
"% ID-string : " << identifier << std::endl
494 <<
"% Method : " << method_and_basis << 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
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++) {
512 if ( rowind[i] < colind[i] )
513 os << colind[i]+1 <<
" " << rowind[i]+1 <<
" " << std::setprecision(10) << values[i] << std::endl;
515 os << rowind[i]+1 <<
" " << colind[i]+1 <<
" " << std::setprecision(10) << values[i] << std::endl;
518 for(
size_t i = 0; i < nvalues; i++) {
519 os << rowind[i]+1 <<
" " << colind[i]+1 <<
" " << std::setprecision(10) << values[i] << std::endl;