source: src/vector.cpp@ 5e588b5

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 5e588b5 was b5bf84, checked in by Tillmann Crueger <crueger@…>, 15 years ago

Added a method that allows to scale all components of one vector with all components of another vector.

  • Property mode set to 100644
File size: 13.2 KB
Line 
1/** \file vector.cpp
2 *
3 * Function implementations for the class vector.
4 *
5 */
6
7#include "Helpers/MemDebug.hpp"
8
9#include "vector.hpp"
10#include "Matrix.hpp"
11#include "verbose.hpp"
12#include "World.hpp"
13#include "Helpers/Assert.hpp"
14#include "Helpers/fast_functions.hpp"
15#include "Exceptions/MathException.hpp"
16
17#include <iostream>
18#include <gsl/gsl_blas.h>
19
20
21using namespace std;
22
23
24/************************************ Functions for class vector ************************************/
25
26/** Constructor of class vector.
27 */
28Vector::Vector()
29{
30 content = gsl_vector_calloc (NDIM);
31};
32
33/**
34 * Copy constructor
35 */
36
37Vector::Vector(const Vector& src)
38{
39 content = gsl_vector_alloc(NDIM);
40 gsl_vector_memcpy(content, src.content);
41}
42
43/** Constructor of class vector.
44 */
45Vector::Vector(const double x1, const double x2, const double x3)
46{
47 content = gsl_vector_alloc(NDIM);
48 gsl_vector_set(content,0,x1);
49 gsl_vector_set(content,1,x2);
50 gsl_vector_set(content,2,x3);
51};
52
53Vector::Vector(gsl_vector *_content) :
54 content(_content)
55{}
56
57/**
58 * Assignment operator
59 */
60Vector& Vector::operator=(const Vector& src){
61 // check for self assignment
62 if(&src!=this){
63 gsl_vector_memcpy(content, src.content);
64 }
65 return *this;
66}
67
68/** Desctructor of class vector.
69 */
70Vector::~Vector() {
71 gsl_vector_free(content);
72};
73
74/** Calculates square of distance between this and another vector.
75 * \param *y array to second vector
76 * \return \f$| x - y |^2\f$
77 */
78double Vector::DistanceSquared(const Vector &y) const
79{
80 double res = 0.;
81 for (int i=NDIM;i--;)
82 res += (at(i)-y[i])*(at(i)-y[i]);
83 return (res);
84};
85
86/** Calculates distance between this and another vector.
87 * \param *y array to second vector
88 * \return \f$| x - y |\f$
89 */
90double Vector::distance(const Vector &y) const
91{
92 return (sqrt(DistanceSquared(y)));
93};
94
95Vector Vector::getClosestPoint(const Vector &point) const{
96 // the closest point to a single point space is always the single point itself
97 return *this;
98}
99
100/** Calculates scalar product between this and another vector.
101 * \param *y array to second vector
102 * \return \f$\langle x, y \rangle\f$
103 */
104double Vector::ScalarProduct(const Vector &y) const
105{
106 double res = 0.;
107 gsl_blas_ddot(content, y.content, &res);
108 return (res);
109};
110
111
112/** Calculates VectorProduct between this and another vector.
113 * -# returns the Product in place of vector from which it was initiated
114 * -# ATTENTION: Only three dim.
115 * \param *y array to vector with which to calculate crossproduct
116 * \return \f$ x \times y \f&
117 */
118void Vector::VectorProduct(const Vector &y)
119{
120 Vector tmp;
121 for(int i=NDIM;i--;)
122 tmp[i] = at((i+1)%NDIM)*y[(i+2)%NDIM] - at((i+2)%NDIM)*y[(i+1)%NDIM];
123 (*this) = tmp;
124};
125
126
127/** projects this vector onto plane defined by \a *y.
128 * \param *y normal vector of plane
129 * \return \f$\langle x, y \rangle\f$
130 */
131void Vector::ProjectOntoPlane(const Vector &y)
132{
133 Vector tmp;
134 tmp = y;
135 tmp.Normalize();
136 tmp.Scale(ScalarProduct(tmp));
137 *this -= tmp;
138};
139
140/** Calculates the minimum distance of this vector to the plane.
141 * \sa Vector::GetDistanceVectorToPlane()
142 * \param *out output stream for debugging
143 * \param *PlaneNormal normal of plane
144 * \param *PlaneOffset offset of plane
145 * \return distance to plane
146 */
147double Vector::DistanceToSpace(const Space &space) const
148{
149 return space.distance(*this);
150};
151
152/** Calculates the projection of a vector onto another \a *y.
153 * \param *y array to second vector
154 */
155void Vector::ProjectIt(const Vector &y)
156{
157 (*this) += (-ScalarProduct(y))*y;
158};
159
160/** Calculates the projection of a vector onto another \a *y.
161 * \param *y array to second vector
162 * \return Vector
163 */
164Vector Vector::Projection(const Vector &y) const
165{
166 Vector helper = y;
167 helper.Scale((ScalarProduct(y)/y.NormSquared()));
168
169 return helper;
170};
171
172/** Calculates norm of this vector.
173 * \return \f$|x|\f$
174 */
175double Vector::Norm() const
176{
177 return (sqrt(NormSquared()));
178};
179
180/** Calculates squared norm of this vector.
181 * \return \f$|x|^2\f$
182 */
183double Vector::NormSquared() const
184{
185 return (ScalarProduct(*this));
186};
187
188/** Normalizes this vector.
189 */
190void Vector::Normalize()
191{
192 double factor = Norm();
193 (*this) *= 1/factor;
194};
195
196/** Zeros all components of this vector.
197 */
198void Vector::Zero()
199{
200 at(0)=at(1)=at(2)=0;
201};
202
203/** Zeros all components of this vector.
204 */
205void Vector::One(const double one)
206{
207 at(0)=at(1)=at(2)=one;
208};
209
210/** Checks whether vector has all components zero.
211 * @return true - vector is zero, false - vector is not
212 */
213bool Vector::IsZero() const
214{
215 return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < MYEPSILON);
216};
217
218/** Checks whether vector has length of 1.
219 * @return true - vector is normalized, false - vector is not
220 */
221bool Vector::IsOne() const
222{
223 return (fabs(Norm() - 1.) < MYEPSILON);
224};
225
226/** Checks whether vector is normal to \a *normal.
227 * @return true - vector is normalized, false - vector is not
228 */
229bool Vector::IsNormalTo(const Vector &normal) const
230{
231 if (ScalarProduct(normal) < MYEPSILON)
232 return true;
233 else
234 return false;
235};
236
237/** Checks whether vector is normal to \a *normal.
238 * @return true - vector is normalized, false - vector is not
239 */
240bool Vector::IsEqualTo(const Vector &a) const
241{
242 bool status = true;
243 for (int i=0;i<NDIM;i++) {
244 if (fabs(at(i) - a[i]) > MYEPSILON)
245 status = false;
246 }
247 return status;
248};
249
250/** Calculates the angle between this and another vector.
251 * \param *y array to second vector
252 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
253 */
254double Vector::Angle(const Vector &y) const
255{
256 double norm1 = Norm(), norm2 = y.Norm();
257 double angle = -1;
258 if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
259 angle = this->ScalarProduct(y)/norm1/norm2;
260 // -1-MYEPSILON occured due to numerical imprecision, catch ...
261 //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
262 if (angle < -1)
263 angle = -1;
264 if (angle > 1)
265 angle = 1;
266 return acos(angle);
267};
268
269
270double& Vector::operator[](size_t i){
271 ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
272 return *gsl_vector_ptr (content, i);
273}
274
275const double& Vector::operator[](size_t i) const{
276 ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
277 return *gsl_vector_ptr (content, i);
278}
279
280double& Vector::at(size_t i){
281 return (*this)[i];
282}
283
284const double& Vector::at(size_t i) const{
285 return (*this)[i];
286}
287
288gsl_vector* Vector::get(){
289 return content;
290}
291
292/** Compares vector \a to vector \a b component-wise.
293 * \param a base vector
294 * \param b vector components to add
295 * \return a == b
296 */
297bool Vector::operator==(const Vector& b) const
298{
299 return IsEqualTo(b);
300};
301
302bool Vector::operator!=(const Vector& b) const
303{
304 return !IsEqualTo(b);
305}
306
307/** Sums vector \a to this lhs component-wise.
308 * \param a base vector
309 * \param b vector components to add
310 * \return lhs + a
311 */
312const Vector& Vector::operator+=(const Vector& b)
313{
314 this->AddVector(b);
315 return *this;
316};
317
318/** Subtracts vector \a from this lhs component-wise.
319 * \param a base vector
320 * \param b vector components to add
321 * \return lhs - a
322 */
323const Vector& Vector::operator-=(const Vector& b)
324{
325 this->SubtractVector(b);
326 return *this;
327};
328
329/** factor each component of \a a times a double \a m.
330 * \param a base vector
331 * \param m factor
332 * \return lhs.x[i] * m
333 */
334const Vector& operator*=(Vector& a, const double m)
335{
336 a.Scale(m);
337 return a;
338};
339
340/** Sums two vectors \a and \b component-wise.
341 * \param a first vector
342 * \param b second vector
343 * \return a + b
344 */
345Vector const Vector::operator+(const Vector& b) const
346{
347 Vector x = *this;
348 x.AddVector(b);
349 return x;
350};
351
352/** Subtracts vector \a from \b component-wise.
353 * \param a first vector
354 * \param b second vector
355 * \return a - b
356 */
357Vector const Vector::operator-(const Vector& b) const
358{
359 Vector x = *this;
360 x.SubtractVector(b);
361 return x;
362};
363
364Vector &Vector::operator*=(const Matrix &mat){
365 (*this) = mat*(*this);
366 return *this;
367}
368
369Vector operator*(const Matrix &mat,const Vector &vec){
370 gsl_vector *res = gsl_vector_calloc(NDIM);
371 gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content, 0.0, res);
372 return Vector(res);
373}
374
375
376/** Factors given vector \a a times \a m.
377 * \param a vector
378 * \param m factor
379 * \return m * a
380 */
381Vector const operator*(const Vector& a, const double m)
382{
383 Vector x(a);
384 x.Scale(m);
385 return x;
386};
387
388/** Factors given vector \a a times \a m.
389 * \param m factor
390 * \param a vector
391 * \return m * a
392 */
393Vector const operator*(const double m, const Vector& a )
394{
395 Vector x(a);
396 x.Scale(m);
397 return x;
398};
399
400ostream& operator<<(ostream& ost, const Vector& m)
401{
402 ost << "(";
403 for (int i=0;i<NDIM;i++) {
404 ost << m[i];
405 if (i != 2)
406 ost << ",";
407 }
408 ost << ")";
409 return ost;
410};
411
412
413void Vector::ScaleAll(const double *factor)
414{
415 for (int i=NDIM;i--;)
416 at(i) *= factor[i];
417};
418
419void Vector::ScaleAll(const Vector &factor){
420 gsl_vector_mul(content, factor.content);
421}
422
423
424void Vector::Scale(const double factor)
425{
426 gsl_vector_scale(content,factor);
427};
428
429std::pair<Vector,Vector> Vector::partition(const Vector &rhs) const{
430 double factor = ScalarProduct(rhs)/rhs.NormSquared();
431 Vector res= factor * rhs;
432 return make_pair(res,(*this)-res);
433}
434
435std::pair<pointset,Vector> Vector::partition(const pointset &points) const{
436 Vector helper = *this;
437 pointset res;
438 for(pointset::const_iterator iter=points.begin();iter!=points.end();++iter){
439 pair<Vector,Vector> currPart = helper.partition(*iter);
440 res.push_back(currPart.first);
441 helper = currPart.second;
442 }
443 return make_pair(res,helper);
444}
445
446/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
447 * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
448 * \param *x1 first vector
449 * \param *x2 second vector
450 * \param *x3 third vector
451 * \param *factors three-component vector with the factor for each given vector
452 */
453void Vector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
454{
455 (*this) = (factors[0]*x1) +
456 (factors[1]*x2) +
457 (factors[2]*x3);
458};
459
460/** Calculates orthonormal vector to one given vectors.
461 * Just subtracts the projection onto the given vector from this vector.
462 * The removed part of the vector is Vector::Projection()
463 * \param *x1 vector
464 * \return true - success, false - vector is zero
465 */
466bool Vector::MakeNormalTo(const Vector &y1)
467{
468 bool result = false;
469 double factor = y1.ScalarProduct(*this)/y1.NormSquared();
470 Vector x1 = factor * y1;
471 SubtractVector(x1);
472 for (int i=NDIM;i--;)
473 result = result || (fabs(at(i)) > MYEPSILON);
474
475 return result;
476};
477
478/** Creates this vector as one of the possible orthonormal ones to the given one.
479 * Just scan how many components of given *vector are unequal to zero and
480 * try to get the skp of both to be zero accordingly.
481 * \param *vector given vector
482 * \return true - success, false - failure (null vector given)
483 */
484bool Vector::GetOneNormalVector(const Vector &GivenVector)
485{
486 int Components[NDIM]; // contains indices of non-zero components
487 int Last = 0; // count the number of non-zero entries in vector
488 int j; // loop variables
489 double norm;
490
491 for (j=NDIM;j--;)
492 Components[j] = -1;
493
494 // in two component-systems we need to find the one position that is zero
495 int zeroPos = -1;
496 // find two components != 0
497 for (j=0;j<NDIM;j++){
498 if (fabs(GivenVector[j]) > MYEPSILON)
499 Components[Last++] = j;
500 else
501 // this our zero Position
502 zeroPos = j;
503 }
504
505 switch(Last) {
506 case 3: // threecomponent system
507 // the position of the zero is arbitrary in three component systems
508 zeroPos = Components[2];
509 case 2: // two component system
510 norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
511 at(zeroPos) = 0.;
512 // in skp both remaining parts shall become zero but with opposite sign and third is zero
513 at(Components[1]) = -1./GivenVector[Components[1]] / norm;
514 at(Components[0]) = 1./GivenVector[Components[0]] / norm;
515 return true;
516 break;
517 case 1: // one component system
518 // set sole non-zero component to 0, and one of the other zero component pendants to 1
519 at((Components[0]+2)%NDIM) = 0.;
520 at((Components[0]+1)%NDIM) = 1.;
521 at(Components[0]) = 0.;
522 return true;
523 break;
524 default:
525 return false;
526 }
527};
528
529/** Adds vector \a *y componentwise.
530 * \param *y vector
531 */
532void Vector::AddVector(const Vector &y)
533{
534 gsl_vector_add(content, y.content);
535}
536
537/** Adds vector \a *y componentwise.
538 * \param *y vector
539 */
540void Vector::SubtractVector(const Vector &y)
541{
542 gsl_vector_sub(content, y.content);
543}
544
545/**
546 * Checks whether this vector is within the parallelepiped defined by the given three vectors and
547 * their offset.
548 *
549 * @param offest for the origin of the parallelepiped
550 * @param three vectors forming the matrix that defines the shape of the parallelpiped
551 */
552bool Vector::IsInParallelepiped(const Vector &offset, const Matrix& _parallelepiped) const
553{
554 Matrix parallelepiped = _parallelepiped.invert();
555 Vector a = parallelepiped * ((*this)-offset);
556 bool isInside = true;
557
558 for (int i=NDIM;i--;)
559 isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
560
561 return isInside;
562}
563
564
565// some comonly used vectors
566const Vector zeroVec(0,0,0);
567const Vector e1(1,0,0);
568const Vector e2(0,1,0);
569const Vector e3(0,0,1);
Note: See TracBrowser for help on using the repository browser.