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_VECTORGENERAL
00037 #define MAT_VECTORGENERAL
00038 #include <iostream>
00039 #include <fstream>
00040 #include <ios>
00041 #include "FileWritable.h"
00042 #include "matrix_proxy.h"
00043 #include "ValidPtr.h"
00044 namespace mat {
00045 template<typename Treal, typename Tvector>
00046 class VectorGeneral : public FileWritable {
00047 public:
00048
00049 inline void resetSizesAndBlocks(SizesAndBlocks const & newRows) {
00050 vectorPtr.haveDataStructureSet(true);
00051 vectorPtr->resetRows(newRows);
00052 }
00053
00054 inline bool is_empty() const {
00055 return !vectorPtr.haveDataStructureGet();
00056 }
00057
00058
00059 VectorGeneral():vectorPtr(new Tvector) {}
00060 explicit VectorGeneral(const VectorGeneral<Treal, Tvector>& other)
00061 :FileWritable(other), vectorPtr(new Tvector) {
00062 if (other.vectorPtr.haveDataStructureGet()) {
00063 vectorPtr.haveDataStructureSet(true);
00064 }
00065 *vectorPtr = *other.vectorPtr;
00066 }
00067
00068 inline void assign_from_full
00069 (std::vector<Treal> const & fullVector,
00070 SizesAndBlocks const & newRows) {
00071 resetSizesAndBlocks(newRows);
00072 this->vectorPtr->assignFromFull(fullVector);
00073 }
00074 inline void fullvector(std::vector<Treal> & fullVector) const {
00075 this->vectorPtr->fullVector(fullVector);
00076 }
00077 VectorGeneral<Treal, Tvector>&
00078 operator=(const VectorGeneral<Treal, Tvector>& other) {
00079 if (other.vectorPtr.haveDataStructureGet()) {
00080 vectorPtr.haveDataStructureSet(true);
00081 }
00082 *this->vectorPtr = *other.vectorPtr;
00083 return *this;
00084 }
00085
00086 inline void clear() {
00087 if (is_empty())
00088
00089
00090 return;
00091 vectorPtr->clear();
00092 }
00093
00094 inline void rand() {
00095 vectorPtr->randomNormalized();
00096 }
00097
00098
00099
00101 template<typename Tmatrix>
00102 VectorGeneral<Treal, Tvector>& operator=
00103 (const XYZ<Treal,
00104 MatrixGeneral<Treal, Tmatrix>,
00105 VectorGeneral<Treal, Tvector> >& smv);
00106
00108 template<typename Tmatrix>
00109 VectorGeneral<Treal, Tvector>& operator+=
00110 (const XYZ<Treal,
00111 MatrixGeneral<Treal, Tmatrix>,
00112 VectorGeneral<Treal, Tvector> >& smv);
00114 template<typename Tmatrix>
00115 VectorGeneral<Treal, Tvector>& operator=
00116 (const XYZpUV<Treal,
00117 MatrixGeneral<Treal, Tmatrix>,
00118 VectorGeneral<Treal, Tvector>,
00119 Treal,
00120 VectorGeneral<Treal, Tvector> >& smvpsv);
00121
00123 template<typename Tmatrix>
00124 VectorGeneral<Treal, Tvector>& operator=
00125 (const XY<MatrixGeneral<Treal, Tmatrix>,
00126 VectorGeneral<Treal, Tvector> >& mv) {
00127 Treal ONE = 1.0;
00128 return this->operator=(XYZ<Treal, MatrixGeneral<Treal, Tmatrix>,
00129 VectorGeneral<Treal, Tvector> >(ONE, mv.A, mv.B,
00130 false, mv.tA, mv.tB));
00131 }
00132
00133
00135 template<typename Tmatrix>
00136 VectorGeneral<Treal, Tvector>& operator=
00137 (const XYZ<Treal,
00138 MatrixSymmetric<Treal, Tmatrix>,
00139 VectorGeneral<Treal, Tvector> >& smv);
00141 template<typename Tmatrix>
00142 VectorGeneral<Treal, Tvector>& operator+=
00143 (const XYZ<Treal,
00144 MatrixSymmetric<Treal, Tmatrix>,
00145 VectorGeneral<Treal, Tvector> >& smv);
00147 template<typename Tmatrix>
00148 VectorGeneral<Treal, Tvector>& operator=
00149 (const XYZpUV<Treal,
00150 MatrixSymmetric<Treal, Tmatrix>,
00151 VectorGeneral<Treal, Tvector>,
00152 Treal,
00153 VectorGeneral<Treal, Tvector> >& smvpsv);
00154
00155
00157 template<typename Tmatrix>
00158 VectorGeneral<Treal, Tvector>& operator=
00159 (const XY<MatrixTriangular<Treal, Tmatrix>,
00160 VectorGeneral<Treal, Tvector> >& mv);
00161
00162
00163
00164 inline Treal eucl() const {
00165 return vectorPtr->eucl();
00166 }
00167
00168 inline VectorGeneral<Treal, Tvector>&
00169 operator*=(Treal const alpha) {
00170 *vectorPtr *= alpha;
00171 return *this;
00172 }
00173
00174 inline VectorGeneral<Treal, Tvector>&
00175 operator=(int const k) {
00176 *vectorPtr = k;
00177 return *this;
00178 }
00179
00181 VectorGeneral<Treal, Tvector>& operator+=
00182 (const XY<Treal, VectorGeneral<Treal, Tvector> >& sv);
00183
00184
00185 inline Tvector const & getVector() const {return *vectorPtr;}
00186
00187 std::string obj_type_id() const {return "VectorGeneral";}
00188 protected:
00189 ValidPtr<Tvector> vectorPtr;
00190
00191 inline void writeToFileProt(std::ofstream & file) const {
00192 if (is_empty())
00193
00194 return;
00195 vectorPtr->writeToFile(file);
00196 }
00197 inline void readFromFileProt(std::ifstream & file) {
00198 if (is_empty())
00199
00200 return;
00201 vectorPtr->readFromFile(file);
00202 }
00203
00204 inline void inMemorySet(bool inMem) {
00205 vectorPtr.inMemorySet(inMem);
00206 }
00207
00208 private:
00209
00210 };
00211
00212
00213
00214
00216 template<typename Treal, typename Tvector>
00217 template<typename Tmatrix>
00218 VectorGeneral<Treal, Tvector>&
00219 VectorGeneral<Treal, Tvector>::operator=
00220 (const XYZ<Treal,
00221 MatrixGeneral<Treal, Tmatrix>,
00222 VectorGeneral<Treal, Tvector> >& smv) {
00223 assert(!smv.tC);
00224 vectorPtr.haveDataStructureSet(true);
00225 if ( this == &smv.C ) {
00226
00227 VectorGeneral<Treal, Tvector> tmp(smv.C);
00228 Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
00229 *tmp.vectorPtr, 0, *this->vectorPtr);
00230 }
00231 else
00232 Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
00233 *smv.C.vectorPtr, 0, *this->vectorPtr);
00234 return *this;
00235 }
00236
00238 template<typename Treal, typename Tvector>
00239 template<typename Tmatrix>
00240 VectorGeneral<Treal, Tvector>&
00241 VectorGeneral<Treal, Tvector>::operator+=
00242 (const XYZ<Treal,
00243 MatrixGeneral<Treal, Tmatrix>,
00244 VectorGeneral<Treal, Tvector> >& smv) {
00245 assert(!smv.tC);
00246 assert(this != &smv.C);
00247 Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
00248 *smv.C.vectorPtr, 1, *this->vectorPtr);
00249 return *this;
00250 }
00251
00252
00254 template<typename Treal, typename Tvector>
00255 template<typename Tmatrix>
00256 VectorGeneral<Treal, Tvector>&
00257 VectorGeneral<Treal, Tvector>::operator=
00258 (const XYZpUV<Treal,
00259 MatrixGeneral<Treal, Tmatrix>,
00260 VectorGeneral<Treal, Tvector>,
00261 Treal,
00262 VectorGeneral<Treal, Tvector> >& smvpsv) {
00263 assert(!smvpsv.tC && !smvpsv.tE);
00264 assert(this != &smvpsv.C);
00265 if (this == &smvpsv.E)
00266 Tvector::gemv(smvpsv.tB, smvpsv.A, smvpsv.B.getMatrix(),
00267 *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
00268 else
00269 throw Failure("VectorGeneral<Treal, Tvector>::operator="
00270 "(const XYZpUV<Treal, "
00271 "MatrixGeneral<Treal, Tmatrix>, "
00272 "VectorGeneral<Treal, Tvector>, "
00273 "Treal, "
00274 "VectorGeneral<Treal, Tvector> >&) : "
00275 "y = alpha * op(A) * x + beta * z "
00276 "not supported for z != y");
00277 return *this;
00278 }
00279
00280
00281
00283 template<typename Treal, typename Tvector>
00284 template<typename Tmatrix>
00285 VectorGeneral<Treal, Tvector>&
00286 VectorGeneral<Treal, Tvector>::operator=
00287 (const XYZ<Treal,
00288 MatrixSymmetric<Treal, Tmatrix>,
00289 VectorGeneral<Treal, Tvector> >& smv) {
00290 assert(!smv.tC);
00291 assert(this != &smv.C);
00292 vectorPtr.haveDataStructureSet(true);
00293 Tvector::symv('U', smv.A, smv.B.getMatrix(),
00294 *smv.C.vectorPtr, 0, *this->vectorPtr);
00295 return *this;
00296 }
00297
00298
00300 template<typename Treal, typename Tvector>
00301 template<typename Tmatrix>
00302 VectorGeneral<Treal, Tvector>&
00303 VectorGeneral<Treal, Tvector>::operator+=
00304 (const XYZ<Treal,
00305 MatrixSymmetric<Treal, Tmatrix>,
00306 VectorGeneral<Treal, Tvector> >& smv) {
00307 assert(!smv.tC);
00308 assert(this != &smv.C);
00309 Tvector::symv('U', smv.A, smv.B.getMatrix(),
00310 *smv.C.vectorPtr, 1, *this->vectorPtr);
00311 return *this;
00312 }
00313
00315 template<typename Treal, typename Tvector>
00316 template<typename Tmatrix>
00317 VectorGeneral<Treal, Tvector>&
00318 VectorGeneral<Treal, Tvector>::operator=
00319 (const XYZpUV<Treal,
00320 MatrixSymmetric<Treal, Tmatrix>,
00321 VectorGeneral<Treal, Tvector>,
00322 Treal,
00323 VectorGeneral<Treal, Tvector> >& smvpsv) {
00324 assert(!smvpsv.tC && !smvpsv.tE);
00325 assert(this != &smvpsv.C);
00326 if (this == &smvpsv.E)
00327 Tvector::symv('U', smvpsv.A, smvpsv.B.getMatrix(),
00328 *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
00329 else
00330 throw Failure("VectorGeneral<Treal, Tvector>::operator="
00331 "(const XYZpUV<Treal, "
00332 "MatrixSymmetric<Treal, Tmatrix>, "
00333 "VectorGeneral<Treal, Tvector>, "
00334 "Treal, "
00335 "VectorGeneral<Treal, Tvector> >&) : "
00336 "y = alpha * A * x + beta * z "
00337 "not supported for z != y");
00338 return *this;
00339 }
00340
00341
00344 template<typename Treal, typename Tvector>
00345 template<typename Tmatrix>
00346 VectorGeneral<Treal, Tvector>&
00347 VectorGeneral<Treal, Tvector>::operator=
00348 (const XY<MatrixTriangular<Treal, Tmatrix>,
00349 VectorGeneral<Treal, Tvector> >& mv) {
00350 assert(!mv.tB);
00351 if (this != &mv.B)
00352 throw Failure("y = A * x not supported for y != x ");
00353 Tvector::trmv('U', mv.tA,
00354 mv.A.getMatrix(),
00355 *this->vectorPtr);
00356 return *this;
00357 }
00358
00359
00360
00362 template<typename Treal, typename Tvector>
00363 VectorGeneral<Treal, Tvector>&
00364 VectorGeneral<Treal, Tvector>::operator+=
00365 (const XY<Treal, VectorGeneral<Treal, Tvector> >& sv) {
00366 assert(!sv.tB);
00367 assert(this != &sv.B);
00368 Tvector::axpy(sv.A, *sv.B.vectorPtr, *this->vectorPtr);
00369 return *this;
00370 }
00371
00372
00373
00374
00378 template<typename Treal, typename Tvector>
00379 Treal operator*(Xtrans<VectorGeneral<Treal, Tvector> > const & xT,
00380 VectorGeneral<Treal, Tvector> const & y) {
00381 if (xT.tA == false)
00382 throw Failure("operator*("
00383 "Xtrans<VectorGeneral<Treal, Tvector> > const &,"
00384 " VectorGeneral<Treal, Tvector> const &): "
00385 "Dimension mismatch in vector operation");
00386 return Tvector::dot(xT.A.getVector(), y.getVector());
00387 }
00388
00389
00390
00391 }
00392 #endif
00393