00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef MATRIX_UTILITIES_HEADER
00029 #define MATRIX_UTILITIES_HEADER
00030
00031 #include "matrix_typedefs.h"
00032 #include "basisinfo.h"
00033
00034 #if 0
00035
00037 Perm* prepare_matrix_permutation(const BasisInfoStruct& basisInfo,
00038 int sparse_block_size,
00039 int factor1, int factor2, int factor3);
00040 #else
00041
00042 mat::SizesAndBlocks prepareMatrixSizesAndBlocks(int n_basis_functions,
00043 int sparse_block_size,
00044 int factor1,
00045 int factor2,
00046 int factor3);
00047
00048 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
00049 int sparse_block_size,
00050 int factor1,
00051 int factor2,
00052 int factor3,
00053 std::vector<int> & permutation);
00054 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
00055 int sparse_block_size,
00056 int factor1,
00057 int factor2,
00058 int factor3,
00059 std::vector<int> & permutation,
00060 std::vector<int> & inversePermutation);
00061
00062 #endif
00063 void fill_matrix_with_random_numbers(int n, symmMatrix & M);
00064
00065 void add_random_diag_perturbation(int n,
00066 symmMatrix & M,
00067 ergo_real eps);
00068
00069 void output_matrix(int n, const ergo_real* matrix, const char* matrixName);
00070
00071 template<class Tmatrix>
00072 ergo_real compute_maxabs_sparse(const Tmatrix & M)
00073 {
00074 return M.maxAbsValue();
00075
00076 }
00077
00078 template<typename RandomAccessIterator>
00079 struct matrix_utilities_CompareClass {
00080 RandomAccessIterator first;
00081 explicit matrix_utilities_CompareClass(RandomAccessIterator firstel)
00082 : first(firstel){}
00083 bool operator() (int i,int j) { return (*(first + i) < *(first + j));}
00084 };
00085
00086 template<typename Tmatrix>
00087 void get_all_nonzeros( Tmatrix const & A,
00088 std::vector<int> const & inversePermutation,
00089 std::vector<int> & rowind,
00090 std::vector<int> & colind,
00091 std::vector<ergo_real> & values) {
00092 rowind.resize(0);
00093 colind.resize(0);
00094 values.resize(0);
00095 size_t nvalues = 0;
00096 size_t nvalues_tmp = A.nvalues();
00097 std::vector<int> rowind_tmp; rowind_tmp.reserve(nvalues_tmp);
00098 std::vector<int> colind_tmp; colind_tmp.reserve(nvalues_tmp);
00099 std::vector<ergo_real> values_tmp; values_tmp.reserve(nvalues_tmp);
00100 A.get_all_values(rowind_tmp,
00101 colind_tmp,
00102 values_tmp,
00103 inversePermutation,
00104 inversePermutation);
00105
00106 for(size_t i = 0; i < nvalues_tmp; i++) {
00107 nvalues += ( values_tmp[i] != 0 );
00108 }
00109 rowind.reserve(nvalues);
00110 colind.reserve(nvalues);
00111 values.reserve(nvalues);
00112
00113 for(size_t i = 0; i < nvalues_tmp; i++) {
00114 if ( values_tmp[i] != 0 ) {
00115 rowind.push_back( rowind_tmp[i] );
00116 colind.push_back( colind_tmp[i] );
00117 values.push_back( values_tmp[i] );
00118 }
00119 }
00120 }
00121
00122
00123
00124 template<typename Tmatrix>
00125 void output_distance_vs_magnitude( BasisInfoStruct const & basisInfo,
00126 Tmatrix const & A,
00127 std::vector<int> const & inversePermutation,
00128 std::string name,
00129 int resolution_r,
00130 int resolution_m
00131 ) {
00132 std::string m_name = name + ".m";
00133 std::ofstream os(m_name.c_str());
00134
00135 int n = basisInfo.noOfBasisFuncs;
00136 std::vector<ergo_real> x(n);
00137 std::vector<ergo_real> y(n);
00138 std::vector<ergo_real> z(n);
00139 for(int i = 0; i < n; i++) {
00140 x[i] = basisInfo.basisFuncList[i].centerCoords[0];
00141 y[i] = basisInfo.basisFuncList[i].centerCoords[1];
00142 z[i] = basisInfo.basisFuncList[i].centerCoords[2];
00143 }
00144
00145 size_t number_of_stored_zeros = 0;
00146 ergo_real minAbsValue = 1e22;
00147 ergo_real maxAbsValue = 0;
00148
00149
00150 size_t nvalues = 0;
00151 std::vector<int> rowind;
00152 std::vector<int> colind;
00153 std::vector<ergo_real> values;
00154 {
00155 std::vector<int> rowind_tmp;
00156 std::vector<int> colind_tmp;
00157 std::vector<ergo_real> values_tmp;
00158 get_all_nonzeros( A, inversePermutation, rowind_tmp, colind_tmp, values_tmp);
00159
00160 bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
00161 if (matrixIsSymmetric) {
00162
00163 size_t nvalues_tmp = values_tmp.size();
00164 rowind.reserve(nvalues_tmp*2);
00165 colind.reserve(nvalues_tmp*2);
00166 values.reserve(nvalues_tmp*2);
00167 for(size_t i = 0; i < nvalues_tmp; i++) {
00168 rowind.push_back( rowind_tmp[i] );
00169 colind.push_back( colind_tmp[i] );
00170 values.push_back( values_tmp[i] );
00171 if ( rowind_tmp[i] != colind_tmp[i] ) {
00172 rowind.push_back( colind_tmp[i] );
00173 colind.push_back( rowind_tmp[i] );
00174 values.push_back( values_tmp[i] );
00175 }
00176 }
00177 }
00178 else {
00179 rowind = rowind_tmp;
00180 colind = colind_tmp;
00181 values = values_tmp;
00182 }
00183
00184 nvalues = values.size();
00185
00186 for(size_t i = 0; i < nvalues; i++) {
00187 ergo_real fabsVal = fabs( values[i] );
00188 values[i] = fabsVal;
00189 minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
00190 maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
00191 }
00192 }
00193
00194 os << "%% Run for example like this: matlab -nosplash -nodesktop -r " << name << std::endl;
00195 os << "number_of_stored_zeros = " << number_of_stored_zeros << ";" << std::endl;
00196 os << "number_of_stored_nonzeros = " << nvalues << ";" << std::endl;
00197
00198
00199 std::vector<ergo_real> distances(nvalues);
00200 for(size_t i = 0; i < nvalues; i++) {
00201 ergo_real diff_x = x[ rowind[i] ] - x[ colind[i] ];
00202 ergo_real diff_y = y[ rowind[i] ] - y[ colind[i] ];
00203 ergo_real diff_z = z[ rowind[i] ] - z[ colind[i] ];
00204 distances[i] = std::sqrt( diff_x * diff_x + diff_y * diff_y + diff_z * diff_z );
00205 }
00206
00207
00208 std::vector<size_t> index(nvalues);
00209 for ( size_t ind = 0; ind < index.size(); ++ind ) {
00210 index[ind] = ind;
00211 }
00212
00213
00214 matrix_utilities_CompareClass<typename std::vector<ergo_real>::const_iterator>
00215 compareDist( distances.begin() );
00216 std::sort ( index.begin(), index.end(), compareDist );
00217
00218
00219 ergo_real minDistance = *std::min_element( distances.begin(), distances.end() );
00220 ergo_real maxDistance = *std::max_element( distances.begin(), distances.end() );
00221
00222 ergo_real rbox_length = (maxDistance - minDistance) / resolution_r;
00223
00224
00225 ergo_real maxMagLog10 = std::log10(maxAbsValue);
00226 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
00227
00228 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
00229
00230 os << "A = [ " << std::endl;
00231
00232 size_t start_ind = 0;
00233 ergo_real r_low = minDistance;
00234 for ( int rbox = 0; rbox < resolution_r; rbox++ ) {
00235 ergo_real r_upp = r_low + rbox_length;
00236
00237 size_t end_ind = start_ind;
00238 while ( end_ind < nvalues && distances[index[end_ind]] < r_upp )
00239 end_ind++;
00240
00241
00242 matrix_utilities_CompareClass<typename std::vector<ergo_real>::const_iterator>
00243 compareMagnitude( values.begin() );
00244 std::sort ( index.begin() + start_ind, index.begin() + end_ind, compareMagnitude );
00245
00246 ergo_real m_low = minMagLog10;
00247 size_t ind_m = start_ind;
00248
00249
00250 while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_low )
00251 ind_m++;
00252 size_t skipped_small = ind_m - start_ind;
00253 os << r_low << " "
00254 << r_upp << " "
00255 << 0 << " "
00256 << std::pow(10,m_low) << " "
00257 << skipped_small
00258 << std::endl;
00259
00260 for ( int mbox = 0; mbox < resolution_m; mbox++ ) {
00261 ergo_real m_upp = m_low + mbox_length;
00262 size_t count = 0;
00263 while ( ind_m < end_ind && std::log10( values[index[ind_m]] ) < m_upp ) {
00264 ind_m++;
00265 count++;
00266 }
00267
00268
00269 os << r_low << " "
00270 << r_upp << " "
00271 << std::pow(10,m_low) << " "
00272 << std::pow(10,m_upp) << " "
00273 << count
00274 << std::endl;
00275 m_low = m_upp;
00276 }
00277
00278 r_low = r_upp;
00279 start_ind = end_ind;
00280 }
00281 os << "];" << std::endl;
00282 os << "B=[];" << std::endl;
00283 os << "for ind = 1 : size(A,1)" << std::endl;
00284 os << " if (A(ind,3) ~= 0)" << std::endl;
00285 os << " B = [B; A(ind,:)];" << std::endl;
00286 os << " end" << std::endl;
00287 os << "end" << std::endl;
00288 os << "%col = jet(101);" << std::endl;
00289 os << "col = gray(101);col=col(end:-1:1,:);" << std::endl;
00290 os << "maxCount = max(B(:,5));" << std::endl;
00291 os << "ax = [0 30 1e-12 1e3]" << std::endl;
00292
00293 os << "fighandle = figure;" << std::endl;
00294 os << "for ind = 1 : size(B,1)" << std::endl;
00295 os << " rmin = B(ind, 1); rmax = B(ind, 2);" << std::endl;
00296 os << " mmin = B(ind, 3); mmax = B(ind, 4);" << std::endl;
00297 os << " colind = round( 1+100 * B(ind,5) / maxCount);" << std::endl;
00298 os << " fill([rmin rmin rmax rmax rmin], [mmin mmax mmax mmin mmin], col(colind,:), 'EdgeColor', col(colind,:) )" << std::endl;
00299 os << " hold on" << std::endl;
00300 os << "end" << std::endl;
00301 os << "set(gca,'YScale','log')" << std::endl;
00302 os << "axis(ax)" << std::endl;
00303 os << "set(gca,'FontSize',16)" << std::endl;
00304 os << "xlabel('Distance')" << std::endl;
00305 os << "ylabel('Magnitude')" << std::endl;
00306 os << "print( fighandle, '-depsc2', '" << name << "')" << std::endl;
00307
00308 os << "fighandle = figure;" << std::endl;
00309 os << "for ind = 1 : size(B,1)" << std::endl;
00310 os << " if (B(ind,5) ~= 0)" << std::endl;
00311 os << " rmin = B(ind, 1); rmax = B(ind, 2);" << std::endl;
00312 os << " mmin = B(ind, 3); mmax = B(ind, 4);" << std::endl;
00313 os << " msize = 3+1*ceil(20 * B(ind,5) / maxCount);" << std::endl;
00314 os << " plot((rmin+rmax)/2,(mmin+mmax)/2,'k.','MarkerSize',msize)" << std::endl;
00315 os << " hold on" << std::endl;
00316 os << " end" << std::endl;
00317 os << "end" << std::endl;
00318 os << "set(gca,'YScale','log')" << std::endl;
00319 os << "axis(ax)" << std::endl;
00320 os << "set(gca,'FontSize',16)" << std::endl;
00321 os << "xlabel('Distance')" << std::endl;
00322 os << "ylabel('Magnitude')" << std::endl;
00323 os << "print( fighandle, '-depsc2', '" << name << "_dots')" << std::endl;
00324 os << "exit(0);" << std::endl;
00325 os.close();
00326 }
00327
00328 template<typename Tmatrix>
00329 void output_magnitude_histogram( Tmatrix const & A,
00330 std::string name,
00331 int resolution_m
00332 ) {
00333 std::string m_name = name + ".m";
00334 std::ofstream os(m_name.c_str());
00335
00336 size_t number_of_stored_zeros = 0;
00337 ergo_real minAbsValue = 1e22;
00338 ergo_real maxAbsValue = 0;
00339
00340
00341 size_t nvalues = 0;
00342 std::vector<int> rowind;
00343 std::vector<int> colind;
00344 std::vector<ergo_real> values;
00345 {
00346
00347 rowind.resize(0);
00348 colind.resize(0);
00349 values.resize(0);
00350 size_t nvalues_tmp = A.nvalues();
00351 std::vector<int> rowind_tmp; rowind_tmp.reserve(nvalues_tmp);
00352 std::vector<int> colind_tmp; colind_tmp.reserve(nvalues_tmp);
00353 std::vector<ergo_real> values_tmp; values_tmp.reserve(nvalues_tmp);
00354 A.get_all_values(rowind_tmp,
00355 colind_tmp,
00356 values_tmp);
00357
00358 for(size_t i = 0; i < nvalues_tmp; i++) {
00359 nvalues += ( values_tmp[i] != 0 );
00360 }
00361
00362 bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
00363 if (matrixIsSymmetric) {
00364
00365 rowind.reserve(nvalues*2);
00366 colind.reserve(nvalues*2);
00367 values.reserve(nvalues*2);
00368
00369 for(size_t i = 0; i < nvalues_tmp; i++) {
00370 if ( values_tmp[i] != 0 ) {
00371 rowind.push_back( rowind_tmp[i] );
00372 colind.push_back( colind_tmp[i] );
00373 values.push_back( values_tmp[i] );
00374 if ( rowind_tmp[i] != colind_tmp[i] ) {
00375 rowind.push_back( colind_tmp[i] );
00376 colind.push_back( rowind_tmp[i] );
00377 values.push_back( values_tmp[i] );
00378 }
00379 }
00380 }
00381 nvalues = values.size();
00382 }
00383 else {
00384 rowind.reserve(nvalues);
00385 colind.reserve(nvalues);
00386 values.reserve(nvalues);
00387
00388 for(size_t i = 0; i < nvalues_tmp; i++) {
00389 if ( values_tmp[i] != 0 ) {
00390 rowind.push_back( rowind_tmp[i] );
00391 colind.push_back( colind_tmp[i] );
00392 values.push_back( values_tmp[i] );
00393 }
00394 }
00395 assert( nvalues == values.size() );
00396 }
00397
00398 for(size_t i = 0; i < nvalues; i++) {
00399 ergo_real fabsVal = fabs( values[i] );
00400 values[i] = fabsVal;
00401 minAbsValue = fabsVal < minAbsValue ? fabsVal : minAbsValue;
00402 maxAbsValue = fabsVal > maxAbsValue ? fabsVal : maxAbsValue;
00403 }
00404 }
00405
00406 os << "%% Run for example like this: matlab -nosplash -nodesktop -r " << name << std::endl;
00407 os << "number_of_stored_zeros = " << number_of_stored_zeros << ";" << std::endl;
00408 os << "number_of_stored_nonzeros = " << nvalues << ";" << std::endl;
00409
00410
00411 std::vector<size_t> index(nvalues);
00412 for ( size_t ind = 0; ind < index.size(); ++ind ) {
00413 index[ind] = ind;
00414 }
00415
00416
00417 ergo_real maxMagLog10 = std::log10(maxAbsValue);
00418 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20;
00419
00420 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m;
00421
00422 os << "A = [ " << std::endl;
00423
00424 matrix_utilities_CompareClass<typename std::vector<ergo_real>::const_iterator>
00425 compareMagnitude( values.begin() );
00426 std::sort ( index.begin(), index.end(), compareMagnitude );
00427
00428 ergo_real m_low = minMagLog10;
00429 size_t ind_m = 0;
00430
00431
00432 while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_low )
00433 ind_m++;
00434 size_t skipped_small = ind_m;
00435 os << 0 << " "
00436 << std::pow(10,m_low) << " "
00437 << skipped_small
00438 << std::endl;
00439
00440 for ( int mbox = 0; mbox < resolution_m; mbox++ ) {
00441 ergo_real m_upp = m_low + mbox_length;
00442 size_t count = 0;
00443 while ( ind_m < nvalues && std::log10( values[index[ind_m]] ) < m_upp ) {
00444 ind_m++;
00445 count++;
00446 }
00447
00448
00449 os << std::pow(10,m_low) << " "
00450 << std::pow(10,m_upp) << " "
00451 << count
00452 << std::endl;
00453 m_low = m_upp;
00454 }
00455 os << "];" << std::endl;
00456
00457 os.close();
00458 }
00459
00460 template<typename Tmatrix>
00461 void write_matrix_in_matrix_market_format( Tmatrix const & A,
00462 std::vector<int> const & inversePermutation,
00463 std::string filename,
00464 std::string identifier,
00465 std::string method_and_basis)
00466 {
00467
00468
00469 size_t nvalues = 0;
00470 std::vector<int> rowind;
00471 std::vector<int> colind;
00472 std::vector<ergo_real> values;
00473 get_all_nonzeros( A, inversePermutation, rowind, colind, values);
00474 nvalues = values.size();
00475
00476
00477 std::string mtx_filename = filename + ".mtx";
00478 std::ofstream os(mtx_filename.c_str());
00479
00480 time_t rawtime;
00481 struct tm * timeinfo;
00482 time ( &rawtime );
00483 timeinfo = localtime ( &rawtime );
00484
00485 std::string matrix_market_matrix_type = "general";
00486 bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
00487 if (matrixIsSymmetric)
00488 matrix_market_matrix_type = "symmetric";
00489 os << "%%MatrixMarket matrix coordinate real " << matrix_market_matrix_type << std::endl
00490 << "%===============================================================================" << std::endl
00491 << "% Generated by the Ergo quantum chemistry program version " << VERSION << " (www.ergoscf.org)" << std::endl
00492 << "% Date : " << asctime (timeinfo)
00493 << "% ID-string : " << identifier << std::endl
00494 << "% Method : " << method_and_basis << std::endl
00495 << "%" << std::endl
00496 << "% MatrixMarket file format:" << std::endl
00497 << "% +-----------------" << std::endl
00498 << "% | % comments" << std::endl
00499 << "% | nrows ncols nentries" << std::endl
00500 << "% | i_1 j_1 A(i_1,j_1)" << std::endl
00501 << "% | i_2 j_2 A(i_1,j_1)" << std::endl
00502 << "% | ..." << std::endl
00503 << "% | i_nentries j_nentries A(i_nentries,j_nentries) " << std::endl
00504 << "% +----------------" << std::endl
00505 << "% Note that indices are 1-based, i.e. A(1,1) is the first element." << std::endl
00506 << "%" << std::endl
00507 << "%===============================================================================" << std::endl;
00508 os << A.get_nrows() << " " << A.get_ncols() << " " << nvalues << std::endl;
00509 if (matrixIsSymmetric)
00510 for(size_t i = 0; i < nvalues; i++) {
00511
00512 if ( rowind[i] < colind[i] )
00513 os << colind[i]+1 << " " << rowind[i]+1 << " " << std::setprecision(10) << values[i] << std::endl;
00514 else
00515 os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << values[i] << std::endl;
00516 }
00517 else
00518 for(size_t i = 0; i < nvalues; i++) {
00519 os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << values[i] << std::endl;
00520 }
00521 os.close();
00522 }
00523
00524
00525 #endif