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. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. 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.
 */
 
/*=========================================================================

  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 oopse {

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