source: src/LinearAlgebra/MatrixContent.cpp@ 0d4424

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

Made gslmatrix obsolete.

  • Property mode set to 100644
File size: 13.7 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_linalg.h>
28#include <gsl/gsl_matrix.h>
29#include <gsl/gsl_multimin.h>
30#include <gsl/gsl_vector.h>
31#include <cmath>
32#include <iostream>
33
34using namespace std;
35
36
37/** Constructor for class MatrixContent.
38 * \param rows number of rows
39 * \param columns number of columns
40 */
41MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
42 rows(_rows),
43 columns(_columns)
44{
45 content = gsl_matrix_calloc(rows, columns);
46}
47
48/** Constructor for class MatrixContent.
49 * \param rows number of rows
50 * \param columns number of columns
51 * \param *src array with components to initialize matrix with
52 */
53MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
54 rows(_rows),
55 columns(_columns)
56{
57 content = gsl_matrix_calloc(rows, columns);
58 set(0,0, src[0]);
59 set(1,0, src[1]);
60 set(2,0, src[2]);
61
62 set(0,1, src[3]);
63 set(1,1, src[4]);
64 set(2,1, src[5]);
65
66 set(0,2, src[6]);
67 set(1,2, src[7]);
68 set(2,2, src[8]);
69}
70
71/** Constructor for class MatrixContent.
72 * We embed the given gls_matrix pointer within this class and set it to NULL
73 * afterwards.
74 * \param *src source gsl_matrix vector to embed within this class
75 */
76MatrixContent::MatrixContent(gsl_matrix *&src) :
77 rows(src->size1),
78 columns(src->size2)
79{
80 content = gsl_matrix_alloc(src->size1, src->size2);
81 gsl_matrix_memcpy(content,src);
82// content = src;
83// src = NULL;
84}
85
86/** Copy constructor for class MatrixContent.
87 * \param &src reference to source MatrixContent
88 */
89MatrixContent::MatrixContent(const MatrixContent &src) :
90 rows(src.rows),
91 columns(src.columns)
92{
93 content = gsl_matrix_alloc(src.rows, src.columns);
94 gsl_matrix_memcpy(content,src.content);
95}
96
97/** Copy constructor for class MatrixContent.
98 * \param *src pointer to source MatrixContent
99 */
100MatrixContent::MatrixContent(const MatrixContent *src) :
101 rows(src->rows),
102 columns(src->columns)
103{
104 ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
105 content = gsl_matrix_alloc(src->rows, src->columns);
106 gsl_matrix_memcpy(content,src->content);
107}
108
109/** Destructor for class MatrixContent.
110 */
111MatrixContent::~MatrixContent()
112{
113 gsl_matrix_free(content);
114}
115
116/** Set matrix to identity.
117 */
118void MatrixContent::setIdentity()
119{
120 for(int i=rows;i--;){
121 for(int j=columns;j--;){
122 set(i,j,i==j);
123 }
124 }
125}
126
127/** Set all matrix components to zero.
128 */
129void MatrixContent::setZero()
130{
131 for(int i=rows;i--;){
132 for(int j=columns;j--;){
133 set(i,j,0.);
134 }
135 }
136}
137
138/** Set all matrix components to a given value.
139 * \param _value value to set each component to
140 */
141void MatrixContent::setValue(double _value)
142{
143 for(int i=rows;i--;){
144 for(int j=columns;j--;){
145 set(i,j,_value);
146 }
147 }
148}
149
150/** Copy operator for MatrixContent with self-assignment check.
151 * \param &src matrix to compare to
152 * \return reference to this
153 */
154MatrixContent &MatrixContent::operator=(const MatrixContent &src)
155{
156 if(&src!=this){
157 gsl_matrix_memcpy(content,src.content);
158 }
159 return *this;
160}
161
162/** Addition operator.
163 * \param &rhs matrix to add
164 * \return reference to this
165 */
166const MatrixContent &MatrixContent::operator+=(const MatrixContent &rhs)
167{
168 gsl_matrix_add(content, rhs.content);
169 return *this;
170}
171
172/** Subtraction operator.
173 * \param &rhs matrix to subtract
174 * \return reference to this
175 */
176const MatrixContent &MatrixContent::operator-=(const MatrixContent &rhs)
177 {
178 gsl_matrix_sub(content, rhs.content);
179 return *this;
180}
181
182/** Multiplication operator.
183 * Note that here matrix have to have same dimensions.
184 * \param &rhs matrix to multiply with
185 * \return reference to this
186 */
187const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
188{
189 ASSERT(rows == rhs.rows,
190 "MatrixContent::operator*=() - row dimension differ: "+toString(rows)+" != "+toString(rhs.rows)+".");
191 ASSERT(columns == rhs.columns,
192 "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.columns)+".");
193 (*this) = (*this)*rhs;
194 return *this;
195}
196
197/** Multiplication with copy operator.
198 * \param &rhs matrix to multiply with
199 * \return reference to newly allocated MatrixContent
200 */
201const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
202{
203 gsl_matrix *res = gsl_matrix_alloc(rows, columns);
204 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
205 // gsl_matrix is taken over by constructor, hence no free
206 MatrixContent tmp(res);
207 gsl_matrix_free(res);
208 return tmp;
209}
210
211/* ========================== Accessing =============================== */
212
213/** Accessor for manipulating component (i,j).
214 * \param i row number
215 * \param j column number
216 * \return reference to component (i,j)
217 */
218double &MatrixContent::at(size_t i, size_t j)
219{
220 ASSERT((i>=0) && (i<rows),
221 "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
222 ASSERT((j>=0) && (j<columns),
223 "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
224 return *gsl_matrix_ptr (content, i, j);
225}
226
227/** Constant accessor for (value of) component (i,j).
228 * \param i row number
229 * \param j column number
230 * \return const component (i,j)
231 */
232const double MatrixContent::at(size_t i, size_t j) const
233{
234 ASSERT((i>=0) && (i<rows),
235 "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
236 ASSERT((j>=0) && (j<columns),
237 "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
238 return gsl_matrix_get(content, i, j);
239}
240
241/** These functions return a pointer to the \a m-th element of a matrix.
242 * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and a null pointer is returned.
243 * \param m index
244 * \return pointer to \a m-th element
245 */
246double *MatrixContent::Pointer(size_t m, size_t n)
247{
248 return gsl_matrix_ptr (content, m, n);
249};
250
251/** These functions return a constant pointer to the \a m-th element of a matrix.
252 * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and a null pointer is returned.
253 * \param m index
254 * \return const pointer to \a m-th element
255 */
256const double *MatrixContent::const_Pointer(size_t m, size_t n) const
257{
258 return gsl_matrix_const_ptr (content, m, n);
259};
260
261/* ========================== Initializing =============================== */
262
263/** Setter for component (i,j).
264 * \param i row numbr
265 * \param j column numnber
266 * \param value value to set componnt (i,j) to
267 */
268void MatrixContent::set(size_t i, size_t j, const double value)
269{
270 ASSERT((i>=0) && (i<rows),
271 "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
272 ASSERT((j>=0) && (j<columns),
273 "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
274 gsl_matrix_set(content,i,j,value);
275}
276
277/** This function sets the matrix from a double array.
278 * Creates a matrix view of the array and performs a memcopy.
279 * \param *x array of values (no dimension check is performed)
280 */
281void MatrixContent::setFromDoubleArray(double * x)
282{
283 gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
284 gsl_matrix_memcpy (content, &m.matrix);
285};
286
287/* ====================== Exchanging elements ============================ */
288/** This function exchanges the \a i-th and \a j-th row of the matrix in-place.
289 * \param i i-th row to swap with ...
290 * \param j ... j-th row to swap against
291 */
292bool MatrixContent::SwapRows(size_t i, size_t j)
293{
294 return (gsl_matrix_swap_rows (content, i, j) == GSL_SUCCESS);
295};
296
297/** This function exchanges the \a i-th and \a j-th column of the matrix in-place.
298 * \param i i-th column to swap with ...
299 * \param j ... j-th column to swap against
300 */
301bool MatrixContent::SwapColumns(size_t i, size_t j)
302{
303 return (gsl_matrix_swap_columns (content, i, j) == GSL_SUCCESS);
304};
305
306/** This function exchanges the \a i-th row and \a j-th column of the matrix in-place.
307 * The matrix must be square for this operation to be possible.
308 * \param i i-th row to swap with ...
309 * \param j ... j-th column to swap against
310 */
311bool MatrixContent::SwapRowColumn(size_t i, size_t j)
312{
313 assert (rows == columns && "The matrix must be square for swapping row against column to be possible.");
314 return (gsl_matrix_swap_rowcol (content, i, j) == GSL_SUCCESS);
315};
316
317/** Return transposed matrix.
318 * \return new matrix that is transposed of this.
319 */
320MatrixContent MatrixContent::transpose() const
321{
322 std::cout << "const MatrixContent::transpose()." << std::endl;
323 gsl_matrix *res = gsl_matrix_alloc(columns, rows); // column and row dimensions exchanged!
324 gsl_matrix_transpose_memcpy(res, content);
325 MatrixContent newContent(res);
326 gsl_matrix_free(res);
327 return newContent;
328}
329
330/** Turn this matrix into its transposed.
331 * Note that this is only possible if rows == columns.
332 */
333MatrixContent &MatrixContent::transpose()
334{
335 std::cout << "MatrixContent::transpose()." << std::endl;
336 ASSERT( rows == columns,
337 "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
338 double tmp;
339 for (size_t i=0;i<rows;i++)
340 for (size_t j=i+1;j<rows;j++) {
341 tmp = at(j,i);
342 at(j,i) = at(i,j);
343 at(i,j) = tmp;
344 }
345 return *this;
346}
347
348/** Transform the matrix to its eigenbasis ans return resulting eigenvalues.
349 * \warn return vector has to be freed'd
350 * \return gsl_vector pointer to vector of eigenvalues
351 */
352gsl_vector* MatrixContent::transformToEigenbasis()
353{
354 ASSERT (rows == columns,
355 "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");
356 gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
357 gsl_vector *eval = gsl_vector_alloc(rows);
358 gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
359 gsl_eigen_symmv(content, eval, evec, T);
360 gsl_eigen_symmv_free(T);
361 gsl_matrix_memcpy(content, evec);
362 gsl_matrix_free(evec);
363 return eval;
364}
365
366/* ============================ Properties ============================== */
367/** Checks whether matrix' elements are strictly null.
368 * \return true - is null, false - else
369 */
370bool MatrixContent::IsNull() const
371{
372 return gsl_matrix_isnull (content);
373};
374
375/** Checks whether matrix' elements are strictly positive.
376 * \return true - is positive, false - else
377 */
378bool MatrixContent::IsPositive() const
379{
380 return gsl_matrix_ispos (content);
381};
382
383/** Checks whether matrix' elements are strictly negative.
384 * \return true - is negative, false - else
385 */
386bool MatrixContent::IsNegative() const
387{
388 return gsl_matrix_isneg (content);
389};
390
391/** Checks whether matrix' elements are strictly non-negative.
392 * \return true - is non-negative, false - else
393 */
394bool MatrixContent::IsNonNegative() const
395{
396 return gsl_matrix_isnonneg (content);
397};
398
399/** This function performs a Cholesky decomposition to determine whether matrix is positive definite.
400 * We check whether GSL returns GSL_EDOM as error, indicating that decomposition failed due to matrix not being positive-definite.
401 * \return true - matrix is positive-definite, false - else
402 */
403bool MatrixContent::IsPositiveDefinite() const
404{
405 if (rows != columns) // only possible for square matrices.
406 return false;
407 else
408 return (gsl_linalg_cholesky_decomp (content) != GSL_EDOM);
409};
410
411
412/** Calculates the determinant of the matrix.
413 * if matrix is square, uses LU decomposition.
414 */
415double MatrixContent::Determinant() const
416{
417 int signum = 0;
418 assert (rows == columns && "Determinant can only be calculated for square matrices.");
419 gsl_permutation *p = gsl_permutation_alloc(rows);
420 gsl_linalg_LU_decomp(content, p, &signum);
421 gsl_permutation_free(p);
422 return gsl_linalg_LU_det(content, signum);
423};
424
425/* ============================= Operators =============================== */
426
427/** Scalar multiplication operator.
428 * \param factor factor to scale with
429 */
430const MatrixContent &MatrixContent::operator*=(const double factor)
431{
432 gsl_matrix_scale(content, factor);
433 return *this;
434}
435
436/** Scalar multiplication and copy operator.
437 * \param factor factor to scale with
438 * \param &mat MatrixContent to scale
439 * \return copied and scaled MatrixContent
440 */
441const MatrixContent operator*(const double factor,const MatrixContent& mat)
442{
443 MatrixContent tmp = mat;
444 tmp*=factor;
445 return tmp;
446}
447
448/** Scalar multiplication and copy operator (with operands exchanged).
449 * \param &mat MatrixContent to scale
450 * \param factor factor to scale with
451 * \return copied and scaled MatrixContent
452 */
453const MatrixContent operator*(const MatrixContent &mat,const double factor)
454{
455 return factor*mat;
456}
457
458/** Equality operator.
459 * Note that we use numerical sensible checking, i.e. with threshold MYEPSILON.
460 * \param &rhs MatrixContent to checks against
461 */
462bool MatrixContent::operator==(const MatrixContent &rhs) const
463 {
464 if ((rows == rhs.rows) && (columns == rhs.columns)) {
465 for(int i=rows;i--;){
466 for(int j=columns;j--;){
467 if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
468 return false;
469 }
470 }
471 }
472 return true;
473 }
474 return false;
475}
476
477Vector operator*(const MatrixContent &mat,const Vector &vec)
478{
479 Vector res;
480 gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res.content->content);
481 return res;
482}
Note: See TracBrowser for help on using the repository browser.