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