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
00029
00038 #ifndef MATRIX_UTILITIES_HEADER
00039 #define MATRIX_UTILITIES_HEADER
00040
00041 #include "matrix_typedefs.h"
00042 #include "basisinfo.h"
00043
00044 #if 0
00045
00047 Perm* prepare_matrix_permutation(const BasisInfoStruct& basisInfo,
00048 int sparse_block_size,
00049 int factor1, int factor2, int factor3);
00050 #else
00051
00052 mat::SizesAndBlocks prepareMatrixSizesAndBlocks(int n_basis_functions,
00053 int sparse_block_size,
00054 int factor1,
00055 int factor2,
00056 int factor3);
00057
00058 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
00059 int sparse_block_size,
00060 int factor1,
00061 int factor2,
00062 int factor3,
00063 std::vector<int> & permutation);
00064 void getMatrixPermutation(const BasisInfoStruct& basisInfo,
00065 int sparse_block_size,
00066 int factor1,
00067 int factor2,
00068 int factor3,
00069 std::vector<int> & permutation,
00070 std::vector<int> & inversePermutation);
00071 void getMatrixPermutationOnlyFactor2(const std::vector<ergo_real> & xcoords,
00072 const std::vector<ergo_real> & ycoords,
00073 const std::vector<ergo_real> & zcoords,
00074 int sparse_block_size_lowest,
00075 int first_factor,
00076 std::vector<int> & permutation,
00077 std::vector<int> & inversePermutation);
00078 void getMatrixPermutationOnlyFactor2(const BasisInfoStruct& basisInfo,
00079 int sparse_block_size_lowest,
00080 int first_factor,
00081 std::vector<int> & permutation,
00082 std::vector<int> & inversePermutation);
00083
00084 #endif
00085 void fill_matrix_with_random_numbers(int n, symmMatrix & M);
00086
00087 void add_random_diag_perturbation(int n,
00088 symmMatrix & M,
00089 ergo_real eps);
00090
00091 bool check_if_matrix_contains_strange_elements(const symmMatrix & M,
00092 std::vector<int> const & inversePermutationHML);
00093
00094 void output_matrix(int n, const ergo_real* matrix, const char* matrixName);
00095
00096 template<class Tmatrix>
00097 ergo_real compute_maxabs_sparse(const Tmatrix & M)
00098 {
00099 return M.maxAbsValue();
00100
00101 }
00102
00103 template<typename RandomAccessIterator>
00104 struct matrix_utilities_CompareClass {
00105 RandomAccessIterator first;
00106 explicit matrix_utilities_CompareClass(RandomAccessIterator firstel)
00107 : first(firstel){}
00108 bool operator() (int i,int j) { return (*(first + i) < *(first + j));}
00109 };
00110
00111 template<typename Tmatrix>
00112 void get_all_nonzeros( Tmatrix const & A,
00113 std::vector<int> const & inversePermutation,
00114 std::vector<int> & rowind,
00115 std::vector<int> & colind,
00116 std::vector<ergo_real> & values) {
00117 rowind.resize(0);
00118 colind.resize(0);
00119 values.resize(0);
00120 size_t nvalues = 0;
00121 size_t nvalues_tmp = A.nvalues();
00122 std::vector<int> rowind_tmp; rowind_tmp.reserve(nvalues_tmp);
00123 std::vector<int> colind_tmp; colind_tmp.reserve(nvalues_tmp);
00124 std::vector<ergo_real> values_tmp; values_tmp.reserve(nvalues_tmp);
00125 A.get_all_values(rowind_tmp,
00126 colind_tmp,
00127 values_tmp,
00128 inversePermutation,
00129 inversePermutation);
00130
00131 for(size_t i = 0; i < nvalues_tmp; i++) {
00132 nvalues += ( values_tmp[i] != 0 );
00133 }
00134 rowind.reserve(nvalues);
00135 colind.reserve(nvalues);
00136 values.reserve(nvalues);
00137
00138 for(size_t i = 0; i < nvalues_tmp; i++) {
00139 if ( values_tmp[i] != 0 ) {
00140 rowind.push_back( rowind_tmp[i] );
00141 colind.push_back( colind_tmp[i] );
00142 values.push_back( values_tmp[i] );
00143 }
00144 }
00145 }
00146
00147
00148 template<typename Tmatrix>
00149 void write_matrix_in_matrix_market_format( Tmatrix const & A,
00150 std::vector<int> const & inversePermutation,
00151 std::string filename,
00152 std::string identifier,
00153 std::string method_and_basis)
00154 {
00155
00156
00157 size_t nvalues = 0;
00158 std::vector<int> rowind;
00159 std::vector<int> colind;
00160 std::vector<ergo_real> values;
00161 get_all_nonzeros( A, inversePermutation, rowind, colind, values);
00162 nvalues = values.size();
00163
00164
00165 std::string mtx_filename = filename + ".mtx";
00166 std::ofstream os(mtx_filename.c_str());
00167
00168 time_t rawtime;
00169 struct tm * timeinfo;
00170 time ( &rawtime );
00171 timeinfo = localtime ( &rawtime );
00172
00173 std::string matrix_market_matrix_type = "general";
00174 bool matrixIsSymmetric = (A.obj_type_id() == "MatrixSymmetric");
00175 if (matrixIsSymmetric)
00176 matrix_market_matrix_type = "symmetric";
00177 os << "%%MatrixMarket matrix coordinate real " << matrix_market_matrix_type << std::endl
00178 << "%===============================================================================" << std::endl
00179 << "% Generated by the Ergo quantum chemistry program version " << VERSION << " (www.ergoscf.org)" << std::endl
00180 << "% Date : " << asctime (timeinfo)
00181 << "% ID-string : " << identifier << std::endl
00182 << "% Method : " << method_and_basis << std::endl
00183 << "%" << std::endl
00184 << "% MatrixMarket file format:" << std::endl
00185 << "% +-----------------" << std::endl
00186 << "% | % comments" << std::endl
00187 << "% | nrows ncols nentries" << std::endl
00188 << "% | i_1 j_1 A(i_1,j_1)" << std::endl
00189 << "% | i_2 j_2 A(i_2,j_2)" << std::endl
00190 << "% | ..." << std::endl
00191 << "% | i_nentries j_nentries A(i_nentries,j_nentries) " << std::endl
00192 << "% +----------------" << std::endl
00193 << "% Note that indices are 1-based, i.e. A(1,1) is the first element." << std::endl
00194 << "%" << std::endl
00195 << "%===============================================================================" << std::endl;
00196 os << A.get_nrows() << " " << A.get_ncols() << " " << nvalues << std::endl;
00197 if (matrixIsSymmetric)
00198 for(size_t i = 0; i < nvalues; i++) {
00199
00200 if ( rowind[i] < colind[i] )
00201 os << colind[i]+1 << " " << rowind[i]+1 << " " << std::setprecision(10) << (double)values[i] << std::endl;
00202 else
00203 os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << (double)values[i] << std::endl;
00204 }
00205 else
00206 for(size_t i = 0; i < nvalues; i++) {
00207 os << rowind[i]+1 << " " << colind[i]+1 << " " << std::setprecision(10) << (double)values[i] << std::endl;
00208 }
00209 os.close();
00210 }
00211
00212
00213 #endif