matrix_utilities.h

Go to the documentation of this file.
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

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1