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