OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
LU.hpp
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
47 Program: Visualization Toolkit
48 Module: $RCSfile: LU.hpp,v $
49
50 Copyright (c) 1993-2003 Ken Martin, Will Schroeder, Bill Lorensen
51 All rights reserved.
52
53 Redistribution and use in source and binary forms, with or without
54 modification, are permitted provided that the following conditions are met:
55
56 * Redistributions of source code must retain the above copyright notice,
57 this list of conditions and the following disclaimer.
58
59 * Redistributions in binary form must reproduce the above copyright notice,
60 this list of conditions and the following disclaimer in the documentation
61 and/or other materials provided with the distribution.
62
63 * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
64 of any contributors may be used to endorse or promote products derived
65 from this software without specific prior written permission.
66
67 * Modified source versions must be plainly marked as such, and must not be
68 misrepresented as being the original software.
69
70 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
71 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
72 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
73 ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
74 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
75 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
76 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
77 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
78 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
79 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
80
81 =========================================================================*/
82#ifndef MATH_LU_HPP
83#define MATH_LU_HPP
84
85#include <iostream>
86#include <limits>
87
88namespace OpenMD {
89
90 /**
91 * Invert input square matrix A into matrix AI.
92 * @param A input square matrix
93 * @param AI output square matrix
94 * @return true if inverse is computed, otherwise return false
95 * @note A is modified during the inversion
96 */
97 template<class MatrixType>
98 bool invertMatrix(MatrixType& A, MatrixType& AI) {
99 using Real = typename MatrixType::ElemType;
100
101 if (A.getNRow() != A.getNCol() || A.getNRow() != AI.getNRow() ||
102 A.getNCol() != AI.getNCol()) {
103 return false;
104 }
105
106 int size = A.getNRow();
107 int *index = NULL, iScratch[10];
108 Real *column = NULL, dScratch[10];
109
110 // Check on allocation of working vectors
111 //
112 if (size <= 10) {
113 index = iScratch;
114 column = dScratch;
115 } else {
116 index = new int[size];
117 column = new Real[size];
118 }
119
120 bool retVal = invertMatrix(A, AI, size, index, column);
121
122 if (size > 10) {
123 delete[] index;
124 delete[] column;
125 }
126
127 return retVal;
128 }
129
130 /**
131 * Invert input square matrix A into matrix AI (Thread safe versions).
132 * @param A input square matrix
133 * @param AI output square matrix
134 * @param size size of the matrix and temporary arrays
135 * @param tmp1Size temporary array
136 * @param tmp2Size temporary array
137 * @return true if inverse is computed, otherwise return false
138 * @note A is modified during the inversion.
139 */
140 template<class MatrixType>
141 bool invertMatrix(MatrixType& A, MatrixType& AI, unsigned int size,
142 int* tmp1Size,
143 typename MatrixType::ElemPoinerType tmp2Size) {
144 if (A.getNRow() != A.getNCol() || A.getNRow() != AI.getNRow() ||
145 A.getNCol() != AI.getNCol() || A.getNRow() != size) {
146 return false;
147 }
148
149 unsigned int i, j;
150
151 //
152 // Factor matrix; then begin solving for inverse one column at a time.
153 // Note: tmp1Size returned value is used later, tmp2Size is just working
154 // memory whose values are not used in LUSolveLinearSystem
155 //
156 if (LUFactorLinearSystem(A, tmp1Size, size, tmp2Size) == 0) {
157 return false;
158 }
159
160 for (j = 0; j < size; j++) {
161 for (i = 0; i < size; i++) {
162 tmp2Size[i] = 0.0;
163 }
164 tmp2Size[j] = 1.0;
165
166 LUSolveLinearSystem(A, tmp1Size, tmp2Size, size);
167
168 for (i = 0; i < size; i++) {
169 AI(i, j) = tmp2Size[i];
170 }
171 }
172
173 return true;
174 }
175
176 /**
177 * Factor linear equations Ax = b using LU decompostion A = LU where L is
178 * lower triangular matrix and U is upper triangular matrix.
179 * @param A input square matrix
180 * @param index pivot indices
181 * @param size size of the matrix and temporary arrays
182 * @param tmpSize temporary array
183 * @return true if inverse is computed, otherwise return false
184 * @note A is modified during the inversion.
185 */
186 template<class MatrixType>
187 int LUFactorLinearSystem(MatrixType& A, int* index, int size,
188 typename MatrixType::ElemPoinerType tmpSize) {
189 using Real = typename MatrixType::ElemType;
190
191 int i, j, k;
192 int maxI = 0;
193 Real largest, temp1, temp2, sum;
194
195 //
196 // Loop over rows to get implicit scaling information
197 //
198 for (i = 0; i < size; i++) {
199 for (largest = 0.0, j = 0; j < size; j++) {
200 if ((temp2 = abs(A(i, j))) > largest) { largest = temp2; }
201 }
202
203 if (largest == 0.0) {
204 std::cerr << "Unable to factor linear system";
205 return 0;
206 }
207 tmpSize[i] = 1.0 / largest;
208 }
209 //
210 // Loop over all columns using Crout's method
211 //
212 for (j = 0; j < size; j++) {
213 for (i = 0; i < j; i++) {
214 sum = A(i, j);
215 for (k = 0; k < i; k++) {
216 sum -= A(i, k) * A(k, j);
217 }
218 A(i, j) = sum;
219 }
220 //
221 // Begin search for largest pivot element
222 //
223 for (largest = 0.0, i = j; i < size; i++) {
224 sum = A(i, j);
225 for (k = 0; k < j; k++) {
226 sum -= A(i, k) * A(k, j);
227 }
228 A(i, j) = sum;
229
230 if ((temp1 = tmpSize[i] * abs(sum)) >= largest) {
231 largest = temp1;
232 maxI = i;
233 }
234 }
235 //
236 // Check for row interchange
237 //
238 if (j != maxI) {
239 for (k = 0; k < size; k++) {
240 temp1 = A(maxI, k);
241 A(maxI, k) = A(j, k);
242 A(j, k) = temp1;
243 }
244 tmpSize[maxI] = tmpSize[j];
245 }
246 //
247 // Divide by pivot element and perform elimination
248 //
249 index[j] = maxI;
250
251 if (abs(A(j, j)) <= std::numeric_limits<Real>::epsilon()) {
252 std::cerr << "Unable to factor linear system";
253 return false;
254 }
255
256 if (j != (size - 1)) {
257 temp1 = 1.0 / A(j, j);
258 for (i = j + 1; i < size; i++) {
259 A(i, j) *= temp1;
260 }
261 }
262 }
263
264 return 1;
265 }
266
267 /**
268 * Solve linear equations Ax = b using LU decompostion A = LU where L is
269 * lower triangular matrix and U is upper triangular matrix.
270 * @param A input square matrix
271 * @param index pivot indices
272 * @param x vector
273 * @param size size of the matrix and temporary arrays
274 * @return true if inverse is computed, otherwise return false
275 * @note A=LU and index[] are generated from method LUFactorLinearSystem).
276 * Also, solution vector is written directly over input load vector.
277 */
278 template<class MatrixType>
279 void LUSolveLinearSystem(MatrixType& A, int* index,
280 typename MatrixType::ElemPoinerType x, int size) {
281 using Real = typename MatrixType::ElemType;
282
283 int i, j, ii, idx;
284 Real sum;
285 //
286 // Proceed with forward and backsubstitution for L and U
287 // matrices. First, forward substitution.
288 //
289 for (ii = -1, i = 0; i < size; i++) {
290 idx = index[i];
291 sum = x[idx];
292 x[idx] = x[i];
293
294 if (ii >= 0) {
295 for (j = ii; j <= (i - 1); j++) {
296 sum -= A(i, j) * x[j];
297 }
298 } else if (sum) {
299 ii = i;
300 }
301
302 x[i] = sum;
303 }
304 //
305 // Now, back substitution
306 //
307 for (i = size - 1; i >= 0; i--) {
308 sum = x[i];
309 for (j = i + 1; j < size; j++) {
310 sum -= A(i, j) * x[j];
311 }
312 x[i] = sum / A(i, i);
313 }
314 }
315} // namespace OpenMD
316
317#endif
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
void LUSolveLinearSystem(MatrixType &A, int *index, typename MatrixType::ElemPoinerType x, int size)
Solve linear equations Ax = b using LU decompostion A = LU where L is lower triangular matrix and U i...
Definition LU.hpp:279
int LUFactorLinearSystem(MatrixType &A, int *index, int size, typename MatrixType::ElemPoinerType tmpSize)
Factor linear equations Ax = b using LU decompostion A = LU where L is lower triangular matrix and U ...
Definition LU.hpp:187
bool invertMatrix(MatrixType &A, MatrixType &AI)
Invert input square matrix A into matrix AI.
Definition LU.hpp:98