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 MAT_GENERAL 00029 #define MAT_GENERAL 00030 #include <cassert> 00031 namespace mat { 00032 00033 00034 00035 template<class Treal> 00036 static Treal maxdiff(const Treal* f1,const Treal* f2,int size) { 00037 Treal diff = 0; 00038 Treal tmpdiff; 00039 for(int i = 0; i < size * size; i++) { 00040 tmpdiff = template_blas_fabs(f1[i] - f2[i]); 00041 if (tmpdiff > 0) 00042 diff = (diff > tmpdiff ? diff : tmpdiff); 00043 } 00044 return diff; 00045 } 00046 00047 template<class Treal> 00048 static Treal maxdiff_tri(const Treal* f1,const Treal* f2,int size) { 00049 Treal diff = 0; 00050 Treal tmpdiff; 00051 for (int col = 0; col < size; col++) 00052 for (int row = 0; row < col + 1; row++) { 00053 tmpdiff = template_blas_fabs(f1[col * size + row] - f2[col * size + row]); 00054 diff = (diff > tmpdiff ? diff : tmpdiff); 00055 } 00056 return diff; 00057 } 00058 00059 00060 template<class Treal> 00061 static Treal frobdiff(const Treal* f1,const Treal* f2,int size) { 00062 Treal diff = 0; 00063 Treal tmp; 00064 for(int i = 0; i < size * size; i++) { 00065 tmp = f1[i] - f2[i]; 00066 diff += tmp * tmp; 00067 } 00068 return template_blas_sqrt(diff); 00069 } 00070 00071 #if 0 00072 template<class T> 00073 static void fileread(T *ptr,int size,FILE*) { 00074 std::cout<<"error reading file"<<std::endl; 00075 } 00076 template<> 00077 void fileread<double>(double *ptr,int size,FILE* file) { 00078 fread(ptr,sizeof(double),size*size,file); 00079 } 00080 template<> 00081 void fileread<float>(float *ptr,int size,FILE* file) { 00082 double* tmpptr=new double [size*size]; 00083 fread(tmpptr,sizeof(double),size*size,file); 00084 for (int i=0;i<size*size;i++) 00085 { 00086 ptr[i]=(float)tmpptr[i]; 00087 } 00088 delete[] tmpptr; 00089 } 00090 #else 00091 template<typename Treal, typename Trealonfile> 00092 static void fileread(Treal *ptr, int size, FILE* file) { 00093 if (sizeof(Trealonfile) == sizeof(Treal)) 00094 fread(ptr,sizeof(Treal),size,file); 00095 else { 00096 Trealonfile* tmpptr=new Trealonfile[size]; 00097 fread(tmpptr,sizeof(Trealonfile),size,file); 00098 for (int i = 0; i < size; i++) { 00099 ptr[i]=(Treal)tmpptr[i]; 00100 } 00101 delete[] tmpptr; 00102 } 00103 } 00104 #endif 00105 00106 template<typename Treal, typename Tmatrix> 00107 static void read_matrix(Tmatrix& A, 00108 char const * const matrixPath, 00109 int const size) { 00110 FILE* matrixfile=fopen(matrixPath,"rb"); 00111 if (!matrixfile) { 00112 throw Failure("read_matrix: Cannot open inputfile"); 00113 } 00114 Treal* matrixfull = new Treal [size*size]; 00115 fileread<Treal, double>(matrixfull, size*size, matrixfile); 00116 /* A must already have built data structure */ 00117 A.assign_from_full(matrixfull, size, size); 00118 delete[] matrixfull; 00119 return; 00120 } 00121 00122 template<typename Treal, typename Trealonfile, typename Tmatrix> 00123 static void read_sparse_matrix(Tmatrix& A, 00124 char const * const rowPath, 00125 char const * const colPath, 00126 char const * const valPath, 00127 int const nval) { 00128 FILE* rowfile=fopen(rowPath,"rb"); 00129 if (!rowfile) { 00130 throw Failure("read_matrix: Cannot open inputfile rowfile"); 00131 } 00132 FILE* colfile=fopen(colPath,"rb"); 00133 if (!colfile) { 00134 throw Failure("read_matrix: Cannot open inputfile colfile"); 00135 } 00136 FILE* valfile=fopen(valPath,"rb"); 00137 if (!valfile) { 00138 throw Failure("read_matrix: Cannot open inputfile valfile"); 00139 } 00140 int* row = new int[nval]; 00141 int* col = new int[nval]; 00142 Treal* val = new Treal[nval]; 00143 fileread<int, int>(row, nval, rowfile); 00144 fileread<int, int>(col, nval, colfile); 00145 fileread<Treal, Trealonfile>(val, nval, valfile); 00146 00147 /* A must already have built data structure */ 00148 A.assign_from_sparse(row, col, val, nval); 00149 #if 0 00150 Treal* compval = new Treal[nval]; 00151 A.get_values(row, col, compval, nval); 00152 Treal maxdiff = 0; 00153 Treal diff; 00154 for (int i = 0; i < nval; i++) { 00155 diff = template_blas_fabs(compval[i] - val[i]); 00156 maxdiff = diff > maxdiff ? diff : maxdiff; 00157 } 00158 std::cout<<"Maxdiff: "<<maxdiff<<std::endl; 00159 #endif 00160 delete[] row; 00161 delete[] col; 00162 delete[] val; 00163 return; 00164 } 00165 00166 template<typename Treal> 00167 static void read_xyz(Treal* x, Treal* y, Treal* z, 00168 char * atomsPath, 00169 int const natoms, 00170 int const size) { 00171 char* atomfile(atomsPath); 00172 std::ifstream input(atomfile); 00173 if (!input) { 00174 throw Failure("read_xyz: Cannot open inputfile"); 00175 } 00176 input >> std::setprecision(10); 00177 Treal* xtmp = new Treal[natoms]; 00178 Treal* ytmp = new Treal[natoms]; 00179 Treal* ztmp = new Treal[natoms]; 00180 int* atomstart = new int[natoms+1]; 00181 for(int i = 0 ; i < natoms ; i++) { 00182 input >> x[i]; 00183 input >> y[i]; 00184 input >> z[i]; 00185 input >> atomstart[i]; 00186 } 00187 atomstart[natoms] = size; 00188 for (int atom = 0; atom < natoms; atom++) 00189 for (int bf = atomstart[atom]; bf < atomstart[atom + 1]; bf++) { 00190 x[bf] = x[atom]; 00191 y[bf] = y[atom]; 00192 z[bf] = z[atom]; 00193 } 00194 delete[] xtmp; 00195 delete[] ytmp; 00196 delete[] ztmp; 00197 delete[] atomstart; 00198 } 00199 } /* end namespace mat */ 00200 00201 #endif