source: src/LinearAlgebra/MatrixContent.cpp@ cca9ef

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 cca9ef was cca9ef, checked in by Frederik Heber <heber@…>, 14 years ago

Renamed Matrix to RealSpaceMatrix.

  • class Matrix only represents 3x3 matrices, whereas we are now going to extend the class MatrixContent to contain arbritrary dimensions.
  • renamed class and file
  • Property mode set to 100644
File size: 9.3 KB
Line 
1/*
2 * MatrixContent.cpp
3 *
4 * Created on: Nov 14, 2010
5 * Author: heber
6 */
7
8
9// include config.h
10#ifdef HAVE_CONFIG_H
11#include <config.h>
12#endif
13
14#include "Helpers/MemDebug.hpp"
15
16#include "LinearAlgebra/RealSpaceMatrix.hpp"
17#include "Helpers/Assert.hpp"
18#include "Exceptions/NotInvertibleException.hpp"
19#include "Helpers/fast_functions.hpp"
20#include "Helpers/Assert.hpp"
21#include "LinearAlgebra/Vector.hpp"
22#include "LinearAlgebra/VectorContent.hpp"
23#include "LinearAlgebra/MatrixContent.hpp"
24
25#include <gsl/gsl_blas.h>
26#include <gsl/gsl_eigen.h>
27#include <gsl/gsl_matrix.h>
28#include <gsl/gsl_multimin.h>
29#include <gsl/gsl_vector.h>
30#include <cmath>
31#include <iostream>
32
33using namespace std;
34
35
36/** Constructor for class MatrixContent.
37 * \param rows number of rows
38 * \param columns number of columns
39 */
40MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
41 rows(_rows),
42 columns(_columns)
43{
44 content = gsl_matrix_calloc(rows, columns);
45}
46
47/** Constructor for class MatrixContent.
48 * \param rows number of rows
49 * \param columns number of columns
50 * \param *src array with components to initialize matrix with
51 */
52MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
53 rows(_rows),
54 columns(_columns)
55{
56 content = gsl_matrix_calloc(rows, columns);
57 set(0,0, src[0]);
58 set(1,0, src[1]);
59 set(2,0, src[2]);
60
61 set(0,1, src[3]);
62 set(1,1, src[4]);
63 set(2,1, src[5]);
64
65 set(0,2, src[6]);
66 set(1,2, src[7]);
67 set(2,2, src[8]);
68}
69
70/** Constructor for class MatrixContent.
71 * We embed the given gls_matrix pointer within this class and set it to NULL
72 * afterwards.
73 * \param *src source gsl_matrix vector to embed within this class
74 */
75MatrixContent::MatrixContent(gsl_matrix *&src) :
76 rows(src->size1),
77 columns(src->size2)
78{
79 content = gsl_matrix_alloc(src->size1, src->size2);
80 gsl_matrix_memcpy(content,src);
81// content = src;
82// src = NULL;
83}
84
85/** Copy constructor for class MatrixContent.
86 * \param &src reference to source MatrixContent
87 */
88MatrixContent::MatrixContent(const MatrixContent &src) :
89 rows(src.rows),
90 columns(src.columns)
91{
92 content = gsl_matrix_alloc(src.rows, src.columns);
93 gsl_matrix_memcpy(content,src.content);
94}
95
96/** Copy constructor for class MatrixContent.
97 * \param *src pointer to source MatrixContent
98 */
99MatrixContent::MatrixContent(const MatrixContent *src) :
100 rows(src->rows),
101 columns(src->columns)
102{
103 ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
104 content = gsl_matrix_alloc(src->rows, src->columns);
105 gsl_matrix_memcpy(content,src->content);
106}
107
108/** Destructor for class MatrixContent.
109 */
110MatrixContent::~MatrixContent()
111{
112 gsl_matrix_free(content);
113}
114
115/** Set matrix to identity.
116 */
117void MatrixContent::setIdentity()
118{
119 for(int i=rows;i--;){
120 for(int j=columns;j--;){
121 set(i,j,i==j);
122 }
123 }
124}
125
126/** Set all matrix components to zero.
127 */
128void MatrixContent::setZero()
129{
130 for(int i=rows;i--;){
131 for(int j=columns;j--;){
132 set(i,j,0.);
133 }
134 }
135}
136
137/** Copy operator for MatrixContent with self-assignment check.
138 * \param &src matrix to compare to
139 * \return reference to this
140 */
141MatrixContent &MatrixContent::operator=(const MatrixContent &src)
142{
143 if(&src!=this){
144 gsl_matrix_memcpy(content,src.content);
145 }
146 return *this;
147}
148
149/** Addition operator.
150 * \param &rhs matrix to add
151 * \return reference to this
152 */
153const MatrixContent &MatrixContent::operator+=(const MatrixContent &rhs)
154{
155 gsl_matrix_add(content, rhs.content);
156 return *this;
157}
158
159/** Subtraction operator.
160 * \param &rhs matrix to subtract
161 * \return reference to this
162 */
163const MatrixContent &MatrixContent::operator-=(const MatrixContent &rhs)
164 {
165 gsl_matrix_sub(content, rhs.content);
166 return *this;
167}
168
169/** Multiplication operator.
170 * Note that here matrix have to have same dimensions.
171 * \param &rhs matrix to multiply with
172 * \return reference to this
173 */
174const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
175{
176 ASSERT(rows == rhs.rows,
177 "MatrixContent::operator*=() - row dimension differ: "+toString(rows)+" != "+toString(rhs.rows)+".");
178 ASSERT(columns == rhs.columns,
179 "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.columns)+".");
180 (*this) = (*this)*rhs;
181 return *this;
182}
183
184/** Multiplication with copy operator.
185 * \param &rhs matrix to multiply with
186 * \return reference to newly allocated MatrixContent
187 */
188const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
189{
190 gsl_matrix *res = gsl_matrix_alloc(rows, columns);
191 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
192 // gsl_matrix is taken over by constructor, hence no free
193 MatrixContent tmp(res);
194 gsl_matrix_free(res);
195 return tmp;
196}
197
198/** Accessor for manipulating component (i,j).
199 * \param i row number
200 * \param j column number
201 * \return reference to component (i,j)
202 */
203double &MatrixContent::at(size_t i, size_t j)
204{
205 ASSERT((i>=0) && (i<rows),
206 "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
207 ASSERT((j>=0) && (j<columns),
208 "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
209 return *gsl_matrix_ptr (content, i, j);
210}
211
212/** Constant accessor for (value of) component (i,j).
213 * \param i row number
214 * \param j column number
215 * \return const component (i,j)
216 */
217const double MatrixContent::at(size_t i, size_t j) const
218{
219 ASSERT((i>=0) && (i<rows),
220 "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
221 ASSERT((j>=0) && (j<columns),
222 "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
223 return gsl_matrix_get(content, i, j);
224}
225
226/** Setter for component (i,j).
227 * \param i row numbr
228 * \param j column numnber
229 * \param value value to set componnt (i,j) to
230 */
231void MatrixContent::set(size_t i, size_t j, const double value)
232{
233 ASSERT((i>=0) && (i<rows),
234 "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
235 ASSERT((j>=0) && (j<columns),
236 "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
237 gsl_matrix_set(content,i,j,value);
238}
239
240/** Return transposed matrix.
241 * \return new matrix that is transposed of this.
242 */
243MatrixContent MatrixContent::transpose() const
244{
245 std::cout << "const MatrixContent::transpose()." << std::endl;
246 gsl_matrix *res = gsl_matrix_alloc(rows, columns);
247 gsl_matrix_transpose_memcpy(res, content);
248 MatrixContent newContent(res);
249 gsl_matrix_free(res);
250 return newContent;
251}
252
253/** Turn this matrix into its transposed.
254 * Note that this is only possible if rows == columns.
255 */
256MatrixContent &MatrixContent::transpose()
257{
258 std::cout << "MatrixContent::transpose()." << std::endl;
259 ASSERT( rows == columns,
260 "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
261 double tmp;
262 for (size_t i=0;i<rows;i++)
263 for (size_t j=i+1;j<rows;j++) {
264 tmp = at(j,i);
265 at(j,i) = at(i,j);
266 at(i,j) = tmp;
267 }
268 return *this;
269}
270
271/** Transform the matrix to its eigenbasis ans return resulting eigenvalues.
272 * \warn return vector has to be freed'd
273 * \return gsl_vector pointer to vector of eigenvalues
274 */
275gsl_vector* MatrixContent::transformToEigenbasis()
276{
277 ASSERT (rows == columns,
278 "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");
279 gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
280 gsl_vector *eval = gsl_vector_alloc(rows);
281 gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
282 gsl_eigen_symmv(content, eval, evec, T);
283 gsl_eigen_symmv_free(T);
284 gsl_matrix_memcpy(content, evec);
285 gsl_matrix_free(evec);
286 return eval;
287}
288
289/** Scalar multiplication operator.
290 * \param factor factor to scale with
291 */
292const MatrixContent &MatrixContent::operator*=(const double factor)
293{
294 gsl_matrix_scale(content, factor);
295 return *this;
296}
297
298/** Scalar multiplication and copy operator.
299 * \param factor factor to scale with
300 * \param &mat MatrixContent to scale
301 * \return copied and scaled MatrixContent
302 */
303const MatrixContent operator*(const double factor,const MatrixContent& mat)
304{
305 MatrixContent tmp = mat;
306 tmp*=factor;
307 return tmp;
308}
309
310/** Scalar multiplication and copy operator (with operands exchanged).
311 * \param &mat MatrixContent to scale
312 * \param factor factor to scale with
313 * \return copied and scaled MatrixContent
314 */
315const MatrixContent operator*(const MatrixContent &mat,const double factor)
316{
317 return factor*mat;
318}
319
320/** Equality operator.
321 * Note that we use numerical sensible checking, i.e. with threshold MYEPSILON.
322 * \param &rhs MatrixContent to checks against
323 */
324bool MatrixContent::operator==(const MatrixContent &rhs) const
325 {
326 if ((rows == rhs.rows) && (columns == rhs.columns)) {
327 for(int i=rows;i--;){
328 for(int j=columns;j--;){
329 if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
330 return false;
331 }
332 }
333 }
334 return true;
335 }
336 return false;
337}
338
339Vector operator*(const MatrixContent &mat,const Vector &vec)
340{
341 Vector res;
342 gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res.content->content);
343 return res;
344}
Note: See TracBrowser for help on using the repository browser.