ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 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 // Count the number of nonzeros 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 // Extract all nonzeros 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 } // end get_all_nonzeros(...) 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 // Get xyz coords for all indices 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 // Get all matrix elements 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 // Also include lower triangle 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 } // end for 00177 } // end if 00178 else { 00179 rowind = rowind_tmp; 00180 colind = colind_tmp; 00181 values = values_tmp; 00182 } // end else 00183 00184 nvalues = values.size(); 00185 // Take absolute value 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 // Get distances for all matrix elements 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 // Index vector 00208 std::vector<size_t> index(nvalues); 00209 for ( size_t ind = 0; ind < index.size(); ++ind ) { 00210 index[ind] = ind; 00211 } 00212 00213 // Sort based on distance 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 // Min and max distances 00219 ergo_real minDistance = *std::min_element( distances.begin(), distances.end() ); 00220 ergo_real maxDistance = *std::max_element( distances.begin(), distances.end() ); 00221 // Size of box in r direction 00222 ergo_real rbox_length = (maxDistance - minDistance) / resolution_r; 00223 00224 // Get max absolute value of A 00225 ergo_real maxMagLog10 = std::log10(maxAbsValue); 00226 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20; 00227 // Size of box in m direction 00228 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m; 00229 00230 os << "A = [ " << std::endl; 00231 // Loop over r boxes 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 // Find end index 00237 size_t end_ind = start_ind; 00238 while ( end_ind < nvalues && distances[index[end_ind]] < r_upp ) 00239 end_ind++; 00240 // Now we have the bounds for box in r-direction 00241 // Sort based on magnitude 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 // Loop over m boxes 00246 ergo_real m_low = minMagLog10; 00247 size_t ind_m = start_ind; 00248 00249 // Skip very small values 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 // Now we have r_low r_upp m_low m_upp count 00268 // Write to stream 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 } // end output_distance_vs_magnitude(...) 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 // Get all matrix elements 00341 size_t nvalues = 0; 00342 std::vector<int> rowind; 00343 std::vector<int> colind; 00344 std::vector<ergo_real> values; 00345 { 00346 // Get all nonzeros 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 // Count the number of nonzeros 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 // Also include lower triangle 00365 rowind.reserve(nvalues*2); 00366 colind.reserve(nvalues*2); 00367 values.reserve(nvalues*2); 00368 // Extract all nonzeros 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 } // end if 00383 else { 00384 rowind.reserve(nvalues); 00385 colind.reserve(nvalues); 00386 values.reserve(nvalues); 00387 // Extract all nonzeros 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 } // end else 00397 // Take absolute value 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 // Index vector 00411 std::vector<size_t> index(nvalues); 00412 for ( size_t ind = 0; ind < index.size(); ++ind ) { 00413 index[ind] = ind; 00414 } 00415 00416 // Get min and max absolute value of A 00417 ergo_real maxMagLog10 = std::log10(maxAbsValue); 00418 ergo_real minMagLog10 = std::log10(minAbsValue) > -20 ? std::log10(minAbsValue) : -20; 00419 // Size of box in m direction 00420 ergo_real mbox_length = (maxMagLog10 - minMagLog10) / resolution_m; 00421 00422 os << "A = [ " << std::endl; 00423 // Sort based on magnitude 00424 matrix_utilities_CompareClass<typename std::vector<ergo_real>::const_iterator> 00425 compareMagnitude( values.begin() ); 00426 std::sort ( index.begin(), index.end(), compareMagnitude ); 00427 // Loop over m boxes 00428 ergo_real m_low = minMagLog10; 00429 size_t ind_m = 0; 00430 00431 // Skip very small values 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 // Now we have m_low m_upp count 00448 // Write to stream 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 } // end output_magnitude_histogram 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 // Get all matrix elements 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 // Now we have all matrix elements 00476 // Open file stream 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) // newline added by asctime 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 // Output lower triangle 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 } // end write_matrix_in_matrix_market_format(...) 00523 00524 00525 #endif