src/math/LU.hpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: LU.hpp,v $

Copyright (c) 1993-2003 Ken Martin, Will Schroeder, Bill Lorensen 
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.

 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

 * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
   of any contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

 * Modified source versions must be plainly marked as such, and must not be
   misrepresented as being the original software.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

=========================================================================*/ 
#ifndef MATH_LU_HPP
#define MATH_LU_HPP

#include "utils/NumericConstant.hpp"

namespace OpenMD {

/**
 * Invert input square matrix A into matrix AI. 
 * @param A input square matrix
 * @param AI output square matrix
 * @return true if inverse is computed, otherwise return false
 * @note A is modified during the inversion
 */
template<class MatrixType>
bool invertMatrix(MatrixType& A, MatrixType& AI)
{
  typedef typename MatrixType::ElemType Real;
  if (A.getNRow() != A.getNCol() || A.getNRow() != AI.getNRow() || A.getNCol() != AI.getNCol()) {
    return false;
  }
  
  int size = A.getNRow();
  int *index=NULL, iScratch[10];
  Real *column=NULL, dScratch[10];

  // Check on allocation of working vectors
  //
  if ( size <= 10 ) {
    index = iScratch;
    column = dScratch;
  } else {
    index = new int[size];
    column = new Real[size];
  }

  bool retVal = invertMatrix(A, AI, size, index, column);

  if ( size > 10 ) {
    delete [] index;
    delete [] column;
  }
  
  return retVal;
}

/**
 * Invert input square matrix A into matrix AI (Thread safe versions). 
 * @param A input square matrix
 * @param AI output square matrix
 * @param size size of the matrix and temporary arrays
 * @param tmp1Size temporary array
 * @param tmp2Size temporary array
 * @return true if inverse is computed, otherwise return false
 * @note A is modified during the inversion.
 */

template<class MatrixType>
bool invertMatrix(MatrixType& A , MatrixType& AI, int size,
                          int *tmp1Size, typename MatrixType::ElemPoinerType tmp2Size)
{
  if (A.getNRow() != A.getNCol() || A.getNRow() != AI.getNRow() || A.getNCol() != AI.getNCol() || A.getNRow() != size) {
    return false;
  }
  
  int i, j;

  //
  // Factor matrix; then begin solving for inverse one column at a time.
  // Note: tmp1Size returned value is used later, tmp2Size is just working
  // memory whose values are not used in LUSolveLinearSystem
  //
  if ( LUFactorLinearSystem(A, tmp1Size, size, tmp2Size) == 0 ){
    return false;
  }
  
  for ( j=0; j < size; j++ ) {
    for ( i=0; i < size; i++ ) {
      tmp2Size[i] = 0.0;
    }
    tmp2Size[j] = 1.0;

    LUSolveLinearSystem(A,tmp1Size,tmp2Size,size);

    for ( i=0; i < size; i++ ) {
      AI(i, j) = tmp2Size[i];
    }
  }

  return true;
}

/**
 * Factor linear equations Ax = b using LU decompostion A = LU where L is
 * lower triangular matrix and U is upper triangular matrix. 
 * @param A input square matrix
 * @param index pivot indices
 * @param size size of the matrix and temporary arrays
 * @param tmpSize temporary array
 * @return true if inverse is computed, otherwise return false
 * @note A is modified during the inversion.
 */
template<class MatrixType>
int LUFactorLinearSystem(MatrixType& A, int *index, int size,
                                  typename MatrixType::ElemPoinerType tmpSize)
{
  typedef typename MatrixType::ElemType Real;
  int i, j, k;
  int maxI = 0;
  Real largest, temp1, temp2, sum;

  //
  // Loop over rows to get implicit scaling information
  //
  for ( i = 0; i < size; i++ ) {
    for ( largest = 0.0, j = 0; j < size; j++ ) {
      if ( (temp2 = fabs(A(i, j))) > largest ) {
        largest = temp2;
      }
    }

    if ( largest == 0.0 ) {
      //vtkGenericWarningMacro(<<"Unable to factor linear system");
      return 0;
    }
      tmpSize[i] = 1.0 / largest;
  }
  //
  // Loop over all columns using Crout's method
  //
  for ( j = 0; j < size; j++ ) {
    for (i = 0; i < j; i++) {
      sum = A(i, j);
      for ( k = 0; k < i; k++ ) {
        sum -= A(i, k) * A(k, j);
      }
      A(i, j) = sum;
    }
    //
    // Begin search for largest pivot element
    //
    for ( largest = 0.0, i = j; i < size; i++ ) {
      sum = A(i, j);
      for ( k = 0; k < j; k++ ) {
        sum -= A(i, k) * A(k, j);
      }
      A(i, j) = sum;

      if ( (temp1 = tmpSize[i]*fabs(sum)) >= largest ) {
        largest = temp1;
        maxI = i;
      }
    }
    //
    // Check for row interchange
    //
    if ( j != maxI ) {
      for ( k = 0; k < size; k++ ) {
        temp1 = A(maxI, k);
        A(maxI, k) = A(j, k);
        A(j, k) = temp1;
      }
      tmpSize[maxI] = tmpSize[j];
    }
    //
    // Divide by pivot element and perform elimination
    //
    index[j] = maxI;

    if ( fabs(A(j, j)) <= OpenMD::NumericConstant::epsilon ) {
      //vtkGenericWarningMacro(<<"Unable to factor linear system");
      return false;
      }

    if ( j != (size-1) ) {
      temp1 = 1.0 / A(j, j);
      for ( i = j + 1; i < size; i++ ) {
        A(i, j) *= temp1;
      }
    }
  }

  return 1;
}

/**
 * Solve linear equations Ax = b using LU decompostion A = LU where L is
 * lower triangular matrix and U is upper triangular matrix. 
 * @param A input square matrix
 * @param index pivot indices
 * @param x vector
 * @param size size of the matrix and temporary arrays
 * @return true if inverse is computed, otherwise return false
 * @note A=LU and index[] are generated from method LUFactorLinearSystem). 
 * Also, solution vector is written directly over input load vector.
 */
template<class MatrixType>
void LUSolveLinearSystem(MatrixType& A, int *index, 
                                  typename MatrixType::ElemPoinerType x, int size)
{
  typedef typename MatrixType::ElemType Real;
  int i, j, ii, idx;
  Real sum;
//
// Proceed with forward and backsubstitution for L and U
// matrices.  First, forward substitution.
//
  for ( ii = -1, i = 0; i < size; i++ ) {
    idx = index[i];
    sum = x[idx];
    x[idx] = x[i];

    if ( ii >= 0 ) {
      for ( j = ii; j <= (i-1); j++ ) {
        sum -= A(i, j)*x[j];
      }
    } else if (sum) {
      ii = i;
    }

    x[i] = sum;
  }
//
// Now, back substitution
//
  for ( i = size-1; i >= 0; i-- ) {
    sum = x[i];
    for ( j = i + 1; j < size; j++ ) {
      sum -= A(i, j)*x[j];
    }
    x[i] = sum / A(i, i);
  }
}

}

#endif
Revision:	1808
Committed:	Mon Oct 22 20:42:10 2012 UTC (12 years, 6 months ago) by gezelter
File size:	9987 byte(s)
Log Message:	A bug fix in the electric field for the new electrostatic code. Also comment fixes for Doxygen
#	User	Rev	Content
1	tim	891	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	tim	891	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	tim	891	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	tim	891	*/
42
43			/*=========================================================================
44
45			Program: Visualization Toolkit
46			Module: $RCSfile: LU.hpp,v $
47
48			Copyright (c) 1993-2003 Ken Martin, Will Schroeder, Bill Lorensen
49			All rights reserved.
50
51			Redistribution and use in source and binary forms, with or without
52			modification, are permitted provided that the following conditions are met:
53
54			* Redistributions of source code must retain the above copyright notice,
55			this list of conditions and the following disclaimer.
56
57			* Redistributions in binary form must reproduce the above copyright notice,
58			this list of conditions and the following disclaimer in the documentation
59			and/or other materials provided with the distribution.
60
61			* Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
62			of any contributors may be used to endorse or promote products derived
63			from this software without specific prior written permission.
64
65			* Modified source versions must be plainly marked as such, and must not be
66			misrepresented as being the original software.
67
68			THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
69			AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
70			IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
71			ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
72			ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
73			DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
74			SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
75			CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
76			OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
77			OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
78
79			=========================================================================*/
80			#ifndef MATH_LU_HPP
81			#define MATH_LU_HPP
82
83			#include "utils/NumericConstant.hpp"
84
85	gezelter	1390	namespace OpenMD {
86	tim	891
87			/**
88			* Invert input square matrix A into matrix AI.
89			* @param A input square matrix
90			* @param AI output square matrix
91			* @return true if inverse is computed, otherwise return false
92			* @note A is modified during the inversion
93			*/
94			template<class MatrixType>
95			bool invertMatrix(MatrixType& A, MatrixType& AI)
96			{
97			typedef typename MatrixType::ElemType Real;
98			if (A.getNRow() != A.getNCol() \|\| A.getNRow() != AI.getNRow() \|\| A.getNCol() != AI.getNCol()) {
99			return false;
100			}
101
102			int size = A.getNRow();
103			int *index=NULL, iScratch[10];
104			Real *column=NULL, dScratch[10];
105
106			// Check on allocation of working vectors
107			//
108			if ( size <= 10 ) {
109			index = iScratch;
110			column = dScratch;
111			} else {
112			index = new int[size];
113			column = new Real[size];
114			}
115
116			bool retVal = invertMatrix(A, AI, size, index, column);
117
118			if ( size > 10 ) {
119			delete [] index;
120			delete [] column;
121			}
122
123			return retVal;
124			}
125
126			/**
127			* Invert input square matrix A into matrix AI (Thread safe versions).
128			* @param A input square matrix
129			* @param AI output square matrix
130			* @param size size of the matrix and temporary arrays
131			* @param tmp1Size temporary array
132			* @param tmp2Size temporary array
133			* @return true if inverse is computed, otherwise return false
134			* @note A is modified during the inversion.
135			*/
136
137			template<class MatrixType>
138			bool invertMatrix(MatrixType& A , MatrixType& AI, int size,
139			int *tmp1Size, typename MatrixType::ElemPoinerType tmp2Size)
140			{
141			if (A.getNRow() != A.getNCol() \|\| A.getNRow() != AI.getNRow() \|\| A.getNCol() != AI.getNCol() \|\| A.getNRow() != size) {
142			return false;
143			}
144
145			int i, j;
146
147			//
148			// Factor matrix; then begin solving for inverse one column at a time.
149			// Note: tmp1Size returned value is used later, tmp2Size is just working
150			// memory whose values are not used in LUSolveLinearSystem
151			//
152			if ( LUFactorLinearSystem(A, tmp1Size, size, tmp2Size) == 0 ){
153			return false;
154			}
155
156			for ( j=0; j < size; j++ ) {
157			for ( i=0; i < size; i++ ) {
158			tmp2Size[i] = 0.0;
159			}
160			tmp2Size[j] = 1.0;
161
162			LUSolveLinearSystem(A,tmp1Size,tmp2Size,size);
163
164			for ( i=0; i < size; i++ ) {
165			AI(i, j) = tmp2Size[i];
166			}
167			}
168
169			return true;
170			}
171
172			/**
173			* Factor linear equations Ax = b using LU decompostion A = LU where L is
174			* lower triangular matrix and U is upper triangular matrix.
175			* @param A input square matrix
176			* @param index pivot indices
177			* @param size size of the matrix and temporary arrays
178			* @param tmpSize temporary array
179			* @return true if inverse is computed, otherwise return false
180			* @note A is modified during the inversion.
181			*/
182			template<class MatrixType>
183			int LUFactorLinearSystem(MatrixType& A, int *index, int size,
184			typename MatrixType::ElemPoinerType tmpSize)
185			{
186			typedef typename MatrixType::ElemType Real;
187			int i, j, k;
188			int maxI = 0;
189			Real largest, temp1, temp2, sum;
190
191			//
192			// Loop over rows to get implicit scaling information
193			//
194			for ( i = 0; i < size; i++ ) {
195			for ( largest = 0.0, j = 0; j < size; j++ ) {
196			if ( (temp2 = fabs(A(i, j))) > largest ) {
197			largest = temp2;
198			}
199			}
200
201			if ( largest == 0.0 ) {
202			//vtkGenericWarningMacro(<<"Unable to factor linear system");
203			return 0;
204			}
205			tmpSize[i] = 1.0 / largest;
206			}
207			//
208			// Loop over all columns using Crout's method
209			//
210			for ( j = 0; j < size; j++ ) {
211			for (i = 0; i < j; i++) {
212			sum = A(i, j);
213			for ( k = 0; k < i; k++ ) {
214			sum -= A(i, k) * A(k, j);
215			}
216			A(i, j) = sum;
217			}
218			//
219			// Begin search for largest pivot element
220			//
221			for ( largest = 0.0, i = j; i < size; i++ ) {
222			sum = A(i, j);
223			for ( k = 0; k < j; k++ ) {
224			sum -= A(i, k) * A(k, j);
225			}
226			A(i, j) = sum;
227
228			if ( (temp1 = tmpSize[i]*fabs(sum)) >= largest ) {
229			largest = temp1;
230			maxI = i;
231			}
232			}
233			//
234			// Check for row interchange
235			//
236			if ( j != maxI ) {
237			for ( k = 0; k < size; k++ ) {
238			temp1 = A(maxI, k);
239			A(maxI, k) = A(j, k);
240			A(j, k) = temp1;
241			}
242			tmpSize[maxI] = tmpSize[j];
243			}
244			//
245			// Divide by pivot element and perform elimination
246			//
247			index[j] = maxI;
248
249	gezelter	1390	if ( fabs(A(j, j)) <= OpenMD::NumericConstant::epsilon ) {
250	tim	891	//vtkGenericWarningMacro(<<"Unable to factor linear system");
251			return false;
252			}
253
254			if ( j != (size-1) ) {
255			temp1 = 1.0 / A(j, j);
256			for ( i = j + 1; i < size; i++ ) {
257			A(i, j) *= temp1;
258			}
259			}
260			}
261
262			return 1;
263			}
264
265			/**
266			* Solve linear equations Ax = b using LU decompostion A = LU where L is
267			* lower triangular matrix and U is upper triangular matrix.
268			* @param A input square matrix
269			* @param index pivot indices
270	gezelter	1808	* @param x vector
271	tim	891	* @param size size of the matrix and temporary arrays
272			* @return true if inverse is computed, otherwise return false
273			* @note A=LU and index[] are generated from method LUFactorLinearSystem).
274			* Also, solution vector is written directly over input load vector.
275			*/
276			template<class MatrixType>
277			void LUSolveLinearSystem(MatrixType& A, int *index,
278			typename MatrixType::ElemPoinerType x, int size)
279			{
280			typedef typename MatrixType::ElemType Real;
281			int i, j, ii, idx;
282			Real sum;
283			//
284			// Proceed with forward and backsubstitution for L and U
285			// matrices. First, forward substitution.
286			//
287			for ( ii = -1, i = 0; i < size; i++ ) {
288			idx = index[i];
289			sum = x[idx];
290			x[idx] = x[i];
291
292			if ( ii >= 0 ) {
293			for ( j = ii; j <= (i-1); j++ ) {
294			sum -= A(i, j)*x[j];
295			}
296			} else if (sum) {
297			ii = i;
298			}
299
300			x[i] = sum;
301			}
302			//
303			// Now, back substitution
304			//
305			for ( i = size-1; i >= 0; i-- ) {
306			sum = x[i];
307			for ( j = i + 1; j < size; j++ ) {
308			sum -= A(i, j)*x[j];
309			}
310			x[i] = sum / A(i, i);
311			}
312			}
313
314			}
315
316			#endif