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
00036 #ifndef MAT_MATRIXBASE
00037 #define MAT_MATRIXBASE
00038 #include <iostream>
00039 #include <fstream>
00040 #include <ios>
00041 #include "FileWritable.h"
00042 #include "matrix_proxy.h"
00043 #include "ValidPtr.h"
00044 #include "SizesAndBlocks.h"
00045 namespace mat {
00046 template<typename Treal, typename Tmatrix>
00047 class MatrixGeneral;
00048 template<typename Treal, typename Tmatrix>
00049 class MatrixSymmetric;
00050 template<typename Treal, typename Tmatrix>
00051 class MatrixTriangular;
00052 template<typename Treal, typename Tvector>
00053 class VectorGeneral;
00054 enum matrix_type {matrix_matr, matrix_symm, matrix_triang};
00055
00066 template<typename Treal, typename Tmatrix>
00067 class MatrixBase : public FileWritable {
00068 public:
00069 friend class MatrixGeneral<Treal, Tmatrix>;
00070 friend class MatrixSymmetric<Treal, Tmatrix>;
00071 friend class MatrixTriangular<Treal, Tmatrix>;
00072
00073
00074 inline void resetSizesAndBlocks(SizesAndBlocks const & newRows,
00075 SizesAndBlocks const & newCols) {
00076 matrixPtr.haveDataStructureSet(true);
00077 matrixPtr->resetRows(newRows);
00078 matrixPtr->resetCols(newCols);
00079 }
00080 inline void getRows(SizesAndBlocks & rowsCopy) const {
00081 matrixPtr->getRows(rowsCopy);
00082 }
00083 inline void getCols(SizesAndBlocks & colsCopy) const {
00084 matrixPtr->getCols(colsCopy);
00085 }
00086
00091 inline bool is_empty() const {
00092 return !matrixPtr.haveDataStructureGet();
00093 }
00094
00095 inline Treal trace() const {
00096 return matrixPtr->trace();
00097 }
00098
00099 inline void add_identity(Treal alpha) {
00100 matrixPtr->addIdentity(alpha);
00101 }
00102 inline MatrixBase<Treal, Tmatrix>& operator*=(Treal const alpha) {
00103 *matrixPtr *= alpha;
00104 return *this;
00105 }
00106
00107 inline bool operator==(int k) const {
00108 if (k == 0)
00109 return *matrixPtr == 0;
00110 else
00111 throw Failure("MatrixBase::operator== only implemented for k == 0");
00112 }
00113
00114
00115
00116 inline void clear() {
00117 if (is_empty())
00118
00119
00120 return;
00121 matrixPtr->clear();
00122 }
00123
00124 inline size_t memory_usage() const {
00125 return matrixPtr->memory_usage();
00126 }
00127
00128 inline void write_to_buffer_count(int& n_bytes) const {
00129 int ib_length = 3;
00130 int vb_length = 0;
00131 this->matrixPtr->write_to_buffer_count(ib_length, vb_length);
00132 n_bytes = vb_length * sizeof(Treal) + ib_length * sizeof(int);
00133 }
00134
00135 #if 1
00136 inline int get_nrows() const {
00137 return matrixPtr->nScalarsRows();
00138 }
00139 inline int get_ncols() const {
00140 return matrixPtr->nScalarsCols();
00141 }
00142 #endif
00143
00144 inline Tmatrix const & getMatrix() const {return *matrixPtr;}
00145 inline Tmatrix & getMatrix() {return *matrixPtr;}
00146
00148 inline Treal maxAbsValue() const {return matrixPtr->maxAbsValue();}
00149
00150 protected:
00151 ValidPtr<Tmatrix> matrixPtr;
00152
00153 MatrixBase():matrixPtr(new Tmatrix) {}
00154 MatrixBase(const MatrixBase<Treal, Tmatrix>& other)
00155 :FileWritable(other), matrixPtr(new Tmatrix) {
00156 matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet());
00157
00158
00159 *matrixPtr = other.matrixPtr.getConstRefForCopying();
00160 matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet());
00161 }
00162
00163 MatrixBase<Treal, Tmatrix>&
00164 operator=(const MatrixBase<Treal, Tmatrix>& other) {
00165 FileWritable::operator=(other);
00166 matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet());
00167
00168
00169 *matrixPtr = other.matrixPtr.getConstRefForCopying();
00170 matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet());
00171 return *this;
00172 }
00173
00174 MatrixBase<Treal, Tmatrix>&
00175 operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) {
00176 if (mt.A.matrixPtr.haveDataStructureGet()) {
00177 matrixPtr.haveDataStructureSet(true);
00178 }
00179 if (mt.tA)
00180 Tmatrix::transpose(*mt.A.matrixPtr, *this->matrixPtr);
00181 else
00182 *this->matrixPtr = *mt.A.matrixPtr;
00183 return *this;
00184
00185 }
00186
00187
00188 void write_to_buffer_base(void* buffer, const int n_bytes,
00189 const matrix_type mattype) const;
00190 void read_from_buffer_base(void* buffer, const int n_bytes,
00191 const matrix_type mattype);
00192
00193 void writeToFileBase(std::ofstream & file,
00194 matrix_type const mattype) const;
00195 void readFromFileBase(std::ifstream & file,
00196 matrix_type const mattype);
00197
00198 std::string obj_type_id() const {return "MatrixBase";}
00199 inline void inMemorySet(bool inMem) {
00200 matrixPtr.inMemorySet(inMem);
00201 }
00202
00203 static void getPermutedIndexes(std::vector<int> const & index,
00204 std::vector<int> const & permutation,
00205 std::vector<int> & newIndex) {
00206 newIndex.resize(index.size());
00207 for (unsigned int i = 0; i < index.size(); ++i)
00208 newIndex[i] = permutation[index[i]];
00209 }
00210
00211
00212 private:
00213
00214 };
00215
00216
00217 template<typename Treal, typename Tmatrix>
00218 void MatrixBase<Treal, Tmatrix>::
00219 writeToFileBase(std::ofstream & file,
00220 matrix_type const mattype) const {
00221 int type = (int)mattype;
00222 file.write((char*)&type,sizeof(int));
00223
00224 if (is_empty())
00225
00226
00227
00228 return;
00229 matrixPtr->writeToFile(file);
00230 }
00231
00232 template<typename Treal, typename Tmatrix>
00233 void MatrixBase<Treal, Tmatrix>::
00234 readFromFileBase(std::ifstream & file,
00235 matrix_type const mattype) {
00236 char type[sizeof(int)];
00237 file.read(type, sizeof(int));
00238 if (((int)*type) != mattype)
00239 throw Failure("MatrixBase<Treal, Tmatrix>::"
00240 "readFromFile(std::ifstream &, "
00241 "matrix_type const): Wrong matrix type");
00242 if (is_empty())
00243
00244 return;
00245 matrixPtr->readFromFile(file);
00246 }
00247
00248
00249
00250 template<typename Treal, typename Tmatrix>
00251 void MatrixBase<Treal, Tmatrix>::
00252 write_to_buffer_base(void* buffer, const int n_bytes,
00253 const matrix_type mattype) const {
00254 int ib_length = 3;
00255
00256 int vb_length = 0;
00257 this->matrixPtr->write_to_buffer_count(ib_length, vb_length);
00258 if (n_bytes >=
00259 (int)(vb_length * sizeof(Treal) + ib_length * sizeof(int))) {
00260 int* int_buf = (int*)buffer;
00261 int_buf[0] = mattype;
00262 int_buf[1] = ib_length;
00263 int_buf[2] = vb_length;
00264 Treal* value_buf = (Treal*)&(int_buf[ib_length]);
00265
00266 int ib_index = 0;
00267 int vb_index = 0;
00268 this->matrixPtr->write_to_buffer(&int_buf[3], ib_length - 3,
00269 value_buf, vb_length,
00270 ib_index, vb_index);
00271 }
00272 else {
00273 throw Failure("MatrixBase::write_to_buffer: Buffer is too small");
00274 }
00275 }
00276
00277 template<typename Treal, typename Tmatrix>
00278 void MatrixBase<Treal, Tmatrix>::
00279 read_from_buffer_base(void* buffer, const int n_bytes,
00280 const matrix_type mattype) {
00281 int* int_buf = (int*)buffer;
00282 if(int_buf[0] == mattype) {
00283 int ib_length = int_buf[1];
00284 int vb_length = int_buf[2];
00285 int ib_index = 0;
00286 int vb_index = 0;
00287 Treal* value_buf = (Treal*)&(int_buf[ib_length]);
00288 this->matrixPtr->read_from_buffer(&int_buf[3], ib_length - 3,
00289 value_buf, vb_length,
00290 ib_index, vb_index);
00291 }
00292 else {
00293 throw Failure("MatrixBase::read_from_buffer: Wrong matrix type");
00294 }
00295 }
00296
00297
00298 }
00299 #endif
00300
00301