OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
SquareMatrix.hpp
Go to the documentation of this file.
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45/**
46 * @file SquareMatrix.hpp
47 * @author Teng Lin
48 * @date 10/11/2004
49 * @version 1.0
50 */
51#ifndef MATH_SQUAREMATRIX_HPP
52#define MATH_SQUAREMATRIX_HPP
53
54#include "math/LU.hpp"
55#include "math/RectMatrix.hpp"
56
57namespace OpenMD {
58
59 /**
60 * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
61 * @brief A square matrix class
62 * \tparam Real the element type
63 * \tparam Dim the dimension of the square matrix
64 */
65 template<typename Real, int Dim>
66 class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
67 public:
68 using ElemType = Real;
69 using ElemPoinerType = Real*;
70
71 /** default constructor */
73 for (unsigned int i = 0; i < Dim; i++)
74 for (unsigned int j = 0; j < Dim; j++)
75 this->data_[i][j] = 0.0;
76 }
77
78 /** Constructs and initializes every element of this matrix to a scalar */
79 SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s) {}
80
81 /** Constructs and initializes from an array */
82 SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array) {}
83
84 /** copy constructor */
86 RectMatrix<Real, Dim, Dim>(m) {}
87
88 /** copy assignment operator */
93
94 /** Returns an identity matrix*/
95
98
99 for (unsigned int i = 0; i < Dim; i++)
100 for (unsigned int j = 0; j < Dim; j++)
101 if (i == j)
102 m(i, j) = 1.0;
103 else
104 m(i, j) = 0.0;
105
106 return m;
107 }
108
109 /**
110 * Returns the inverse of this matrix.
111 */
113 // LU sacrifices the original matrix, so create a copy:
114 SquareMatrix<Real, Dim> tmp(this);
116
117 // Use LU-factorization to invert the matrix:
118 invertMatrix(tmp, m);
119 delete tmp;
120 return m;
121 }
122
123 /**
124 * Returns the determinant of this matrix.
125 */
126 Real determinant() const {
127 Real det;
128 // Base case : if matrix contains single element
129 if (Dim == 1) return this->data_[0][0];
130
131 int sign = 1; // To store sign multiplier
132
133 // Iterate for each element of first row
134 for (int f = 0; f < Dim; f++) {
135 // Getting Cofactor of A[0][f]
136 SquareMatrix<Real, Dim - 1> temp = cofactor(this, 0, f);
137 det += sign * this->data_[0][f] * temp.determinant();
138
139 // terms are to be added with alternate sign
140 sign = -sign;
141 }
142 return det;
143 }
144
145 SquareMatrix<Real, Dim - 1> cofactor(int p, int q) {
146 SquareMatrix<Real, Dim - 1> m;
147
148 int i = 0, j = 0;
149 // Looping for each element of the matrix
150 for (unsigned int row = 0; row < Dim; row++) {
151 for (unsigned int col = 0; col < Dim; col++) {
152 // Copying into temporary matrix only those element
153 // which are not in given row and column
154 if (row != p && col != q) {
155 m(i, j++) = this->data_[row][col];
156
157 // Row is filled, so increase row index and
158 // reset col index
159 if (j == Dim - 1) {
160 j = 0;
161 i++;
162 }
163 }
164 }
165 }
166 return m;
167 }
168
169 /** Returns the trace of this matrix. */
170 Real trace() const {
171 Real tmp = 0;
172
173 for (unsigned int i = 0; i < Dim; i++)
174 tmp += this->data_[i][i];
175
176 return tmp;
177 }
178
179 /** Tests if this matrix is symmetrix. */
180 bool isSymmetric() const {
181 for (unsigned int i = 0; i < Dim - 1; i++)
182 for (unsigned int j = i; j < Dim; j++)
183 if (fabs(this->data_[i][j] - this->data_[j][i]) > epsilon)
184 return false;
185
186 return true;
187 }
188
189 /** Tests if this matrix is orthogonal. */
192
193 tmp = *this * transpose();
194
195 return tmp.isDiagonal();
196 }
197
198 /** Tests if this matrix is diagonal. */
199 bool isDiagonal() const {
200 for (unsigned int i = 0; i < Dim; i++)
201 for (unsigned int j = 0; j < Dim; j++)
202 if (i != j && fabs(this->data_[i][j]) > epsilon) return false;
203
204 return true;
205 }
206
207 /**
208 * Returns a column vector that contains the elements from the
209 * diagonal of m in the order R(0) = m(0,0), R(1) = m(1,1), and so
210 * on.
211 */
213 Vector<Real, Dim> result;
214 for (unsigned int i = 0; i < Dim; i++) {
215 result(i) = this->data_[i][i];
216 }
217 return result;
218 }
219
220 /** Tests if this matrix is the unit matrix. */
221 bool isUnitMatrix() const {
222 if (!isDiagonal()) return false;
223
224 for (unsigned int i = 0; i < Dim; i++)
225 if (fabs(this->data_[i][i] - 1) > epsilon) return false;
226
227 return true;
228 }
229
230 /** Return the transpose of this matrix */
233
234 for (unsigned int i = 0; i < Dim; i++)
235 for (unsigned int j = 0; j < Dim; j++)
236 result(j, i) = this->data_[i][j];
237
238 return result;
239 }
240
241 /** @todo need implementation */
242 void diagonalize() {
243 // jacobi(m, eigenValues, ortMat);
244 }
245
246 /**
247 * Jacobi iteration routines for computing eigenvalues/eigenvectors of
248 * real symmetric matrix
249 *
250 * @return true if success, otherwise return false
251 * @param a symmetric matrix whose eigenvectors are to be computed. On
252 * return, the matrix is overwritten
253 * @param d will contain the eigenvalues of the matrix On return of this
254 * function
255 * @param v the columns of this matrix will contain the eigenvectors. The
256 * eigenvectors are normalized and mutually orthogonal.
257 */
258
261 }; // end SquareMatrix
262
263 /*=========================================================================
264
265 Program: Visualization Toolkit
266 Module: $RCSfile: SquareMatrix.hpp,v $
267
268 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
269 All rights reserved.
270 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
271
272 This software is distributed WITHOUT ANY WARRANTY; without even
273 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
274 PURPOSE. See the above copyright notice for more information.
275
276 =========================================================================*/
277
278#define VTK_ROTATE(a, i, j, k, l) \
279 g = a(i, j); \
280 h = a(k, l); \
281 a(i, j) = g - s * (h + g * tau); \
282 a(k, l) = h + s * (g - h * tau)
283
284#define VTK_MAX_ROTATIONS 20
285
286 // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
287 // real symmetric matrix. Square nxn matrix a; size of matrix in n;
288 // output eigenvalues in w; and output eigenvectors in v. Resulting
289 // eigenvalues/vectors are sorted in decreasing order; eigenvectors are
290 // normalized.
291 template<typename Real, int Dim>
295 const int n = Dim;
296 int i, j, k, iq, ip, numPos;
297 Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
298 Real bspace[4], zspace[4];
299 Real* b = bspace;
300 Real* z = zspace;
301
302 // only allocate memory if the matrix is large
303 if (n > 4) {
304 b = new Real[n];
305 z = new Real[n];
306 }
307
308 // initialize
309 for (ip = 0; ip < n; ip++) {
310 for (iq = 0; iq < n; iq++) {
311 v(ip, iq) = 0.0;
312 }
313 v(ip, ip) = 1.0;
314 }
315 for (ip = 0; ip < n; ip++) {
316 b[ip] = w[ip] = a(ip, ip);
317 z[ip] = 0.0;
318 }
319
320 // begin rotation sequence
321 for (i = 0; i < VTK_MAX_ROTATIONS; i++) {
322 sm = 0.0;
323 for (ip = 0; ip < n - 1; ip++) {
324 for (iq = ip + 1; iq < n; iq++) {
325 sm += fabs(a(ip, iq));
326 }
327 }
328 if (sm == 0.0) { break; }
329
330 if (i < 3) { // first 3 sweeps
331 tresh = 0.2 * sm / (n * n);
332 } else {
333 tresh = 0.0;
334 }
335
336 for (ip = 0; ip < n - 1; ip++) {
337 for (iq = ip + 1; iq < n; iq++) {
338 g = 100.0 * fabs(a(ip, iq));
339
340 // after 4 sweeps
341 if (i > 3 && (fabs(w[ip]) + g) == fabs(w[ip]) &&
342 (fabs(w[iq]) + g) == fabs(w[iq])) {
343 a(ip, iq) = 0.0;
344 } else if (fabs(a(ip, iq)) > tresh) {
345 h = w[iq] - w[ip];
346 if ((fabs(h) + g) == fabs(h)) {
347 t = (a(ip, iq)) / h;
348 } else {
349 theta = 0.5 * h / (a(ip, iq));
350 t = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta));
351 if (theta < 0.0) { t = -t; }
352 }
353 c = 1.0 / sqrt(1 + t * t);
354 s = t * c;
355 tau = s / (1.0 + c);
356 h = t * a(ip, iq);
357 z[ip] -= h;
358 z[iq] += h;
359 w[ip] -= h;
360 w[iq] += h;
361 a(ip, iq) = 0.0;
362
363 // ip already shifted left by 1 unit
364 for (j = 0; j <= ip - 1; j++) {
365 VTK_ROTATE(a, j, ip, j, iq);
366 }
367 // ip and iq already shifted left by 1 unit
368 for (j = ip + 1; j <= iq - 1; j++) {
369 VTK_ROTATE(a, ip, j, j, iq);
370 }
371 // iq already shifted left by 1 unit
372 for (j = iq + 1; j < n; j++) {
373 VTK_ROTATE(a, ip, j, iq, j);
374 }
375 for (j = 0; j < n; j++) {
376 VTK_ROTATE(v, j, ip, j, iq);
377 }
378 }
379 }
380 }
381
382 for (ip = 0; ip < n; ip++) {
383 b[ip] += z[ip];
384 w[ip] = b[ip];
385 z[ip] = 0.0;
386 }
387 }
388
389 //// this is NEVER called
390 if (i >= VTK_MAX_ROTATIONS) {
391 std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions"
392 << std::endl;
393 if (n > 4) {
394 delete[] b;
395 delete[] z;
396 }
397 return 0;
398 }
399
400 // sort eigenfunctions these changes do not affect accuracy
401 for (j = 0; j < n - 1; j++) { // boundary incorrect
402 k = j;
403 tmp = w[k];
404 for (i = j + 1; i < n; i++) { // boundary incorrect, shifted already
405 if (w[i] >= tmp) { // why exchage if same?
406 k = i;
407 tmp = w[k];
408 }
409 }
410 if (k != j) {
411 w[k] = w[j];
412 w[j] = tmp;
413 for (i = 0; i < n; i++) {
414 tmp = v(i, j);
415 v(i, j) = v(i, k);
416 v(i, k) = tmp;
417 }
418 }
419 }
420 // insure eigenvector consistency (i.e., Jacobi can compute vectors that
421 // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
422 // reek havoc in hyperstreamline/other stuff. We will select the most
423 // positive eigenvector.
424 int ceil_half_n = (n >> 1) + (n & 1);
425 for (j = 0; j < n; j++) {
426 for (numPos = 0, i = 0; i < n; i++) {
427 if (v(i, j) >= 0.0) { numPos++; }
428 }
429 // if ( numPos < ceil(RealType(n)/RealType(2.0)) )
430 if (numPos < ceil_half_n) {
431 for (i = 0; i < n; i++) {
432 v(i, j) *= -1.0;
433 }
434 }
435 }
436
437 if (n > 4) {
438 delete[] b;
439 delete[] z;
440 }
441 return 1;
442 }
443
445} // namespace OpenMD
446
447#endif // MATH_SQUAREMATRIX_HPP
rectangular matrix class
RectMatrix< Real, Row, Col > & operator=(const RectMatrix< Real, Row, Col > &m)
copy assignment operator
A square matrix class.
SquareMatrix()
default constructor
Vector< Real, Dim > diagonals() const
Returns a column vector that contains the elements from the diagonal of m in the order R(0) = m(0,...
SquareMatrix< Real, Dim > & operator=(const RectMatrix< Real, Dim, Dim > &m)
copy assignment operator
Real trace() const
Returns the trace of this matrix.
bool isOrthogonal()
Tests if this matrix is orthogonal.
static SquareMatrix< Real, Dim > identity()
Returns an identity matrix.
Real determinant() const
Returns the determinant of this matrix.
bool isSymmetric() const
Tests if this matrix is symmetrix.
bool isDiagonal() const
Tests if this matrix is diagonal.
SquareMatrix(const RectMatrix< Real, Dim, Dim > &m)
copy constructor
SquareMatrix< Real, Dim > inverse()
Returns the inverse of this matrix.
SquareMatrix(Real s)
Constructs and initializes every element of this matrix to a scalar.
SquareMatrix(Real *array)
Constructs and initializes from an array.
bool isUnitMatrix() const
Tests if this matrix is the unit matrix.
SquareMatrix< Real, Dim > transpose() const
Return the transpose of this matrix.
static int jacobi(SquareMatrix< Real, Dim > &a, Vector< Real, Dim > &d, SquareMatrix< Real, Dim > &v)
Jacobi iteration routines for computing eigenvalues/eigenvectors of real symmetric matrix.
Fix length vector class.
Definition Vector.hpp:78
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
bool invertMatrix(MatrixType &A, MatrixType &AI)
Invert input square matrix A into matrix AI.
Definition LU.hpp:98