www.pudn.com > Kriging 算法实现 2D和3D地图等高线.zip > Numeric.h, change:2003-12-07,size:2484b

// Numeric.h: interface for the Numeric class. 
// 
////////////////////////////////////////////////////////////////////// 
 
#if !defined(AFX_NUMERIC_H__960148E2_5114_4465_AECE_6FC7433EC9CD__INCLUDED_) 
#define AFX_NUMERIC_H__960148E2_5114_4465_AECE_6FC7433EC9CD__INCLUDED_ 
 
#if _MSC_VER > 1000 
#pragma once 
#endif // _MSC_VER > 1000 
 
#include "Matrix.h" 
#include "BaseException.h" 
#include <vector> 
 
#define TINY_VALUE			1.5e-16 
 
class NumericException : public BaseException 
{ 
public: 
	NumericException() throw() : BaseException(_T("Numeric Exception")) {} 
	NumericException(std::string message) throw() : BaseException(message) {} 
}; 
 
template<class T> 
void LUDecompose(TMatrix<T>& A, std::vector<int>& Permutation, int& d) throw(NumericException) 
{ 
	int n = A.GetHeight(); 
	vector<T> vv(n); 
	Permutation.resize(n); 
 
	d=1; 
	 
	T amax; 
	for(int i=0; i<n; i++) { 
		amax = 0.0; 
		for(int j=0; j<n; j++) 
			if(fabs(A(i, j)) > amax) 
				amax = fabs(A(i, j)); 
 
		if(amax  TINY_VALUE) 
			throw NumericException(); 
 
		vv[i] = 1.0 / amax; 
	} 
 
	T sum, dum; 
	int imax; 
	for(int j=0; j<n; j++) { 
		for (i=0; i<j; i++) { 
			sum = A(i, j); 
			for(int k=0; k<i; k++) 
				sum -= A(i, k) * A(k, j); 
			A(i, j) = sum; 
		} 
		amax = 0.0; 
 
		for(i=j; i<n; i++) { 
			sum = A(i, j); 
			for(int k=0; k<j; k++) 
				sum -= A(i, k) * A(k, j); 
 
			A(i, j) = sum; 
			dum = vv[i] * fabs(sum); 
 
			if(dum >= amax) { 
				imax = i; 
				amax = dum; 
			} 
		} 
 
		if(j != imax) { 
			for(int k=0; k<n; k++) { 
				dum = A(imax, k); 
				A(imax, k) = A(j, k); 
				A(j, k) = dum; 
			} 
			d = -d; 
			vv[imax] = vv[j]; 
		} 
		Permutation[j] = imax; 
 
		if(fabs(A(j, j))  TINY_VALUE) 
			A(j, j) = TINY_VALUE; 
 
		if(j != n) { 
			dum = 1.0 / A(j, j); 
			for(i=j+1; i<n; i++) 
				A(i, j) *= dum; 
		} 
	}  
} 
 
template<class T> 
void LUBackSub(TMatrix<T>& A, std::vector<int>& Permutation, std::vector<T>& B) throw(NumericException) 
{ 
	int n = A.GetHeight(); 
	T sum; 
	int ii = 0; 
	int ll; 
	for(int i=0; i<n; i++) { 
		ll = Permutation[i]; 
		sum = B[ll]; 
		B[ll] = B[i]; 
		if(ii != 0) 
			for(int j=ii; j<i; j++) 
				sum -= A(i, j) * B[j]; 
		else if(sum != 0.0)\ 
			ii = i; 
		B[i] = sum; 
	} 
 
	for(i=n-1; i>=0; i--) { 
		sum = B[i]; 
		if(i n) { 
			for(int j=i+1; j<n; j++) 
				sum -= A(i, j) * B[j]; 
		} 
		B[i] = sum / A(i, i); 
	} 
}   
 
#endif // !defined(AFX_NUMERIC_H__960148E2_5114_4465_AECE_6FC7433EC9CD__INCLUDED_)