source: src/gslmatrix.cpp@ 88104f

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 88104f was 0a4f7f, checked in by Tillmann Crueger <crueger@…>, 15 years ago

Made data internal data-structure of vector class private

  • Replaced occurences of access to internals with operator
  • moved Vector-class into LinAlg-Module
  • Reworked Vector to allow clean modularization
  • Added Plane class to describe arbitrary planes in 3d space
  • Property mode set to 100644
File size: 7.2 KB
Line 
1/*
2 * gslmatrix.cpp
3 *
4 * Created on: Jan 8, 2010
5 * Author: heber
6 */
7
8using namespace std;
9
10#include "gslmatrix.hpp"
11#include "helpers.hpp"
12#include "Helpers/fast_functions.hpp"
13
14#include <cassert>
15#include <gsl/gsl_linalg.h>
16
17/** Constructor of class GSLMatrix.
18 * Allocates GSL structures
19 * \param m dimension of matrix
20 */
21GSLMatrix::GSLMatrix(size_t m, size_t n) : rows(m), columns(n)
22{
23 matrix = gsl_matrix_calloc(rows, columns);
24};
25
26/** Copy constructor of class GSLMatrix.
27 * Allocates GSL structures and copies components from \a *src.
28 * \param *src source matrix
29 */
30GSLMatrix::GSLMatrix(const GSLMatrix * const src) : rows(src->rows), columns(src->columns)
31{
32 matrix = gsl_matrix_alloc(rows, columns);
33 gsl_matrix_memcpy (matrix, src->matrix);
34};
35
36/** Destructor of class GSLMatrix.
37 * Frees GSL structures
38 */
39GSLMatrix::~GSLMatrix()
40{
41 gsl_matrix_free(matrix);
42 rows = 0;
43 columns = 0;
44};
45
46/** Assignment operator.
47 * \param &rhs right hand side
48 * \return object itself
49 */
50GSLMatrix& GSLMatrix::operator=(const GSLMatrix& rhs)
51{
52 if (this == &rhs) // not necessary here, but identity assignment check is generally a good idea
53 return *this;
54 assert(rows == rhs.rows && columns == rhs.columns && "Number of rows and columns do not match!");
55
56 gsl_matrix_memcpy (matrix, rhs.matrix);
57
58 return *this;
59};
60
61/* ============================ Accessing =============================== */
62/** This function sets the matrix from a double array.
63 * Creates a matrix view of the array and performs a memcopy.
64 * \param *x array of values (no dimension check is performed)
65 */
66void GSLMatrix::SetFromDoubleArray(double * x)
67{
68 gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
69 gsl_matrix_memcpy (matrix, &m.matrix);
70};
71
72/** This function returns the i-th element of a matrix.
73 * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and 0 is returned.
74 * \param m row index
75 * \param n colum index
76 * \return (m,n)-th element of matrix
77 */
78double GSLMatrix::Get(size_t m, size_t n)
79{
80 return gsl_matrix_get (matrix, m, n);
81};
82
83/** This function sets the value of the \a m -th element of a matrix to \a x.
84 * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked.
85 * \param m row index
86 * \param m column index
87 * \param x value to set element (m,n) to
88 */
89void GSLMatrix::Set(size_t m, size_t n, double x)
90{
91 gsl_matrix_set (matrix, m, n, x);
92};
93
94/** These functions return a pointer to the \a m-th element of a matrix.
95 * 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.
96 * \param m index
97 * \return pointer to \a m-th element
98 */
99double *GSLMatrix::Pointer(size_t m, size_t n)
100{
101 return gsl_matrix_ptr (matrix, m, n);
102};
103
104/** These functions return a constant pointer to the \a m-th element of a matrix.
105 * 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.
106 * \param m index
107 * \return const pointer to \a m-th element
108 */
109const double *GSLMatrix::const_Pointer(size_t m, size_t n)
110{
111 return gsl_matrix_const_ptr (matrix, m, n);
112};
113
114/* ========================== Initializing =============================== */
115/** This function sets all the elements of the matrix to the value \a x.
116 * \param *x
117 */
118void GSLMatrix::SetAll(double x)
119{
120 gsl_matrix_set_all (matrix, x);
121};
122
123/** This function sets all the elements of the matrix to zero.
124 */
125void GSLMatrix::SetZero()
126{
127 gsl_matrix_set_zero (matrix);
128};
129
130/** This function sets the elements of the matrix to the corresponding elements of the identity matrix, \f$m(i,j) = \delta(i,j)\f$, i.e. a unit diagonal with all off-diagonal elements zero.
131 * This applies to both square and rectangular matrices.
132 */
133void GSLMatrix::SetIdentity()
134{
135 gsl_matrix_set_identity (matrix);
136};
137
138/* ====================== Exchanging elements ============================ */
139/** This function exchanges the \a i-th and \a j-th row of the matrix in-place.
140 * \param i i-th row to swap with ...
141 * \param j ... j-th row to swap against
142 */
143bool GSLMatrix::SwapRows(size_t i, size_t j)
144{
145 return (gsl_matrix_swap_rows (matrix, i, j) == GSL_SUCCESS);
146};
147
148/** This function exchanges the \a i-th and \a j-th column of the matrix in-place.
149 * \param i i-th column to swap with ...
150 * \param j ... j-th column to swap against
151 */
152bool GSLMatrix::SwapColumns(size_t i, size_t j)
153{
154 return (gsl_matrix_swap_columns (matrix, i, j) == GSL_SUCCESS);
155};
156
157/** This function exchanges the \a i-th row and \a j-th column of the matrix in-place.
158 * The matrix must be square for this operation to be possible.
159 * \param i i-th row to swap with ...
160 * \param j ... j-th column to swap against
161 */
162bool GSLMatrix::SwapRowColumn(size_t i, size_t j)
163{
164 assert (rows == columns && "The matrix must be square for swapping row against column to be possible.");
165 return (gsl_matrix_swap_rowcol (matrix, i, j) == GSL_SUCCESS);
166};
167
168/** This function transposes the matrix.
169 * Note that the function is extended to the non-square case, where the matrix is re-allocated and copied.
170 */
171bool GSLMatrix::Transpose()
172{
173 if (rows == columns)// if square, use GSL
174 return (gsl_matrix_transpose (matrix) == GSL_SUCCESS);
175 else { // otherwise we have to copy a bit
176 gsl_matrix *dest = gsl_matrix_alloc(columns, rows);
177 for (size_t i=0;i<rows; i++)
178 for (size_t j=0;j<columns;j++) {
179 gsl_matrix_set(dest, j,i, gsl_matrix_get(matrix, i,j) );
180 }
181 gsl_matrix_free(matrix);
182 matrix = dest;
183 flip(rows, columns);
184 return true;
185 }
186};
187
188/* ============================ Properties ============================== */
189/** Checks whether matrix' elements are strictly null.
190 * \return true - is null, false - else
191 */
192bool GSLMatrix::IsNull()
193{
194 return gsl_matrix_isnull (matrix);
195};
196
197/** Checks whether matrix' elements are strictly positive.
198 * \return true - is positive, false - else
199 */
200bool GSLMatrix::IsPositive()
201{
202 return gsl_matrix_ispos (matrix);
203};
204
205/** Checks whether matrix' elements are strictly negative.
206 * \return true - is negative, false - else
207 */
208bool GSLMatrix::IsNegative()
209{
210 return gsl_matrix_isneg (matrix);
211};
212
213/** Checks whether matrix' elements are strictly non-negative.
214 * \return true - is non-negative, false - else
215 */
216bool GSLMatrix::IsNonNegative()
217{
218 return gsl_matrix_isnonneg (matrix);
219};
220
221/** This function performs a Cholesky decomposition to determine whether matrix is positive definite.
222 * We check whether GSL returns GSL_EDOM as error, indicating that decomposition failed due to matrix not being positive-definite.
223 * \return true - matrix is positive-definite, false - else
224 */
225bool GSLMatrix::IsPositiveDefinite()
226{
227 if (rows != columns) // only possible for square matrices.
228 return false;
229 else
230 return (gsl_linalg_cholesky_decomp (matrix) != GSL_EDOM);
231};
232
233
234/** Calculates the determinant of the matrix.
235 * if matrix is square, uses LU decomposition.
236 */
237double GSLMatrix::Determinant()
238{
239 int signum = 0;
240 assert (rows == columns && "Determinant can only be calculated for square matrices.");
241 gsl_permutation *p = gsl_permutation_alloc(rows);
242 gsl_linalg_LU_decomp(matrix, p, &signum);
243 gsl_permutation_free(p);
244 return gsl_linalg_LU_det(matrix, signum);
245};
246
Note: See TracBrowser for help on using the repository browser.