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 #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
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
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 }
00200
00201 #endif