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