general.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 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

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1