for(j = 1; j < n; j++) a(0, j) /= a(0, 0);
    
	for(i = 1; i < n; i++)
	{
        for(j = 1; j <= i; j++)
            a(i, i) = a(i, i) - a((j - 1), i) * a((j - 1), i);
    
		a(i, i) = sqrt(a(i, i));
        if(i != (n - 1))
		{
			for(j = i + 1; j < n; j++)
			{
                for(k = 1; k <= i; k++)
					a(i, j) -= a((k - 1), i) * a((k - 1), j);
                
				a(i, j) /= a(i, i);
			}
		}
	}
    
	for(j = 0; j < m; j++)
	{
		d(0, j) = d(0, j) / a(0, 0);
        for(i = 1; i < n; i++)
		{
            for(k = 1; k <= i; k++)
				d(i, j) -= a((k - 1), i) * d((k - 1), j);
            
			d(i, j) /= a(i, i);
		}
	}
    
	for(j = 0; j < m; j++)
	{
        d((n - 1), j) = d((n - 1), j) / a((n - 1), (n - 1));
		
		for(k = n - 1; k >= 1; k --)
		{
		    for(i = k; i < n; i++)
		        d((k - 1), j) -= a((k - 1), i) * d(i, j);
        
            d((k - 1), j) /= a((k - 1), (k - 1));			
		}
	}
    
	return (1);	//運行成功，正常返回
}

//求解大型稀疏方程組的全選主元高斯-約當消去法
template <class _Ty>
int LE_SparseEuationTotalChoiceGaussJordan(matrix<_Ty> & a, valarray<_Ty> & b)
{ 
	int i, j, k, is;
    _Ty d, t;

    int n = a.GetColNum();		//方程組階數
    valarray<int> js(n);
	
    for(k = 0; k < n; k++)
	{ 
		d = 0.0;
    	for(i = k; i < n; i++)
			for(j = k; j < n; j++)
			{
				t = Abs(a(i, j));
		
				if(t > d)
				{ 
					d = t;
					js[k] = j; 
					is = i;
				}
			}
		
		if( FloatEqual(d, 0.0) )	//主元為0
		{					
			cout<<"fail 1"<<endl;
			return (0);
		}
		
		if(is != k)
		{ 
			for(j = k; j < n; j++)
			{
				t = a(k, j);
				a(k, j) = a(is, j);
				a(is, j) = t;
			}
			
			t = b[k];
			b[k] = b[is];
			b[is] = t;
		}
		
		if(js[k] != k)
			for(i = 0; i < n; i++)
			{
				t = a(i, k);
				a(i, k) = a(i, js[k]); 
				a(i, js[k]) = t;
			}
		
		t = a(k, k);
		
		for(j = k + 1; j < n; j++)
		{
			if(a(k, j) != 0.0)
				a(k, j) = a(k, j) / t;
		}
		
		b[k] = b[k] / t;
		for(j = k + 1; j < n; j++)
		{ 
			if(a(k, j) != 0.0)
			{
				for(i = 0; i < n; i++)
				{
					if((i != k) && (a(i, k) != 0.0))
					{ 
						a(i, j) = a(i, j) - a(i, k) * a(k, j);
					}
				}
			}
		}
		
		for(i = 0; i <  n; i++)
		{
			if((i != k) && (a(i, k) != 0.0))
				b[i] = b[i] - a(i, k) * b[k];
		}
	}
	
	for(k = n - 1; k >= 0; k--)
		if(k != js[k])
		{
			t = b[k];
			b[k] = b[js[k]]; 
			b[js[k]] = t;
		}
	   
	return (1);	//運行成功，正常返回
}

//求解托伯利茲方程組的列文遜方法
template <class _Ty>
int LE_ToeplitzEuationLevinson(valarray<_Ty>& t,  valarray<_Ty>& b, valarray<_Ty>& x)
{
	int i, j, k;
    _Ty a, beta, q, c, h;
	bool yn;

    int n = t.size();			//方程組階數
	valarray<_Ty> y(n);
	valarray<_Ty> s(n);

	a = t[0];
	yn = FloatEqual(a, 0.0);
    if(yn)
	{
		cout << "fail 1" << endl;
		return(-1);
	}
    
	y[0] = 1.0;
	x[0] = b[0] / a;

    for(k = 1; k < n; k++)
	{
		beta = 0.0;
		q = 0.0;
        
		for(j = 0; j < k; j++)
		{
			beta = beta + y[j] * t[j + 1];
            q = q + x[j] * t[k - j];
		}
        
		yn = FloatEqual(a, 0.0);
		if(yn)		
		{
			cout << "fail 2" << endl;
			return(-1);
		}
        
		c = -beta / a; 
		s[0] = c * y[k - 1]; 
		y[k] = y[k - 1];
        
		if(k != 1)
			for(i = 1; i < k; i++)
				s[i] = y[i - 1] + c * y[k - i - 1];
			
		a = a + c * beta;
        
		yn = FloatEqual(a, 0.0);
		if(yn)
		{
			cout << "fail 3" << endl;
			return (-1);
		}

        h = (b[k] - q) / a;
        
		for(i = 0; i < k ; i++)
		{
			x[i] = x[i] + h * s[i];
			y[i] = s[i];
		}
        
		x[k] = h * y[k];
	}
  
	return (1);	//運行成功，正常返回
}

//高斯-賽德爾迭代
template <class _Ty>
int LE_GaussSeidelIteration(matrix<_Ty>& a, valarray<_Ty>& b, 
											valarray<_Ty>& x, _Ty eps)
{
	int i, j;
    _Ty p, t, s, q;

	int n = a.GetColNum();		//方程組階數
    
	for(i = 0; i < n; i++)
	{	
		p = 0.0;
		x[i] = 0.0;
        
		for(j = 0; j < n; j++)
			if(i != j)
				p = p + Abs(a(i, j));
        
		if(p >= Abs(a(i, i)) )	//主對角線不占絕對優勢
		{
			cout<<"fail 1"<<endl;
			return(-1);
		}
	}
    
	p = eps + 1.0;
    while(p > eps  || FloatEqual(p, eps))
	{
		p = 0.0;
        for(i = 0; i < n; i++)
		{
			t = x[i];
			s = 0.0;
            for(j = 0; j < n; j++)
				if(j != i) 
					s = s + a(i, j) * x[j];
			
			x[i] = (b[i] - s) / a(i, i);
            q = Abs(x[i] - t) / (1.0 + Abs(x[i]));
            if(q > p) p = q;
		}
	}

	return (1);	//運行成功，正常返回
}
    
//求解對稱正定方程組的共軛梯度法
template <class _Ty>
int LE_SymmetryRegularEuationConjugateGradient(matrix<_Ty> & a, 
						valarray<_Ty> & b, _Ty eps, valarray<_Ty> & x)
{
	int i, k;
    _Ty  alpha, beta, d, e;
	
	if(MatrixSymmetryRegular(a,1) != 2)	//判別矩陣a是否對稱正定
		return(-1);						//矩陣a不對稱正定
    
	int n = a.GetColNum();				//方程組階數
    valarray<_Ty> r(n);
    
	matrix<double> q(n, 1);
    matrix<double> p(n, 1);
	matrix<double> s(n, 1);
    matrix<double> xx(n, 1);
    
	for(i = 0; i < n; i++) 
		xx(i, 0) = x[i];

    for(i = 0; i < n; i++)
	{
		xx(i,0) = 0.0;
		p(i, 0) = b[i];
		r[i] = b[i]; 
	}
    
	i = 0;

    while(i < n)
	{		
		MatrixMultiply(s, a, p);
        d = 0.0;
		e = 0.0;
        
		for(k = 0; k < n; k++)
		{
			d = d + p(k, 0) * b[k];
			e = e + p(k, 0) * s(k,0); 
		}
        
		alpha = d / e;
        
		for(k = 0; k < n; k++)
			xx(k,0) = xx(k, 0) + alpha * p(k, 0);
        		
		MatrixMultiply(q, a, xx);
        
		d = 0.0;
        
		for(k = 0; k < n; k++)
		{
			r[k] = b[k] - q(k, 0);
			d = d + r[k] * s(k, 0);
		}
        
		beta = d / e;
		d = 0.0;
        
		for(k = 0; k < n; k++) d = d + r[k] * r[k];
        
		d = sqrt(d);
        
		if(d < eps) 
		{
			for(i = 0; i < n; i++) 
				x[i] = xx(i, 0);

		}

        for(k = 0; k < n; k++) 
 			p(k, 0) = r[k] - beta * p(k, 0);
        
		i++;
	}

	for(i = 0; i < n; i++) 
		x[i] = xx(i, 0);

	return (1);	//運行成功，正常返回
}

//求解線性最小二乘問題的豪斯荷爾德變換法
template <class _Ty>
int LE_LinearLeastSquareHouseholder(matrix<_Ty>& a, valarray<_Ty>& b, matrix<_Ty>& q)
{
	int i,j;
    _Ty d;

    int n = a.GetColNum();		//系數矩陣a的行數，m>=n
	int m = a.GetRowNum();		//系數矩陣a的列數，n<=m
    
	valarray<_Ty> c(n);

    i = MatrixQR(a, q);			//一般m×n階的實矩陣QR分解函數

    if(i == 0)
	{
		return (-1);
	}

	for(i = 0; i < n; i++)
	{
		d = 0.0;

        for(j = 0; j < m; j ++)
				d = d + q(j ,i) * b[j];
        
		c[i] = d;
	}
    
	b[n - 1] = c[n - 1] / a((n - 1), (n - 1));
	
    for(i = n - 2; i >= 0; i--)
	{
		d = 0.0;

        for(j = i + 1; j < n; j++)
				d = d + a(i,j) * b[j];
        
		b[i] = (c[i] - d) / a(i, i);
	}
	return (1);	//運行成功，正常返回
}

//求解線性最小二乘問題的廣義逆法
template <class _Ty>
int LE_LinearLeastSquareGeneralizedInverse(matrix<_Ty>& a, 
				valarray<_Ty>& b, valarray<_Ty>& x, matrix<_Ty>& aa, 
								_Ty eps, matrix<_Ty>& u, matrix<_Ty>& v)
{
	int i, j, ii;
       
	int n = a.GetColNum();	//系數矩陣a的列數
	int m = a.GetRowNum();	//系數矩陣a的行數
    
	int ka = (n > m) ? n : m;	
	
	ka++;

	ii = GeneralizedInversionSingularValue(a, aa, eps, u, v);
    
	if(ii < 0)  return(-1);
    
	for(i = 0; i < n; i ++)
	{
		x[i] = 0.0;
        for(j = 0; j < m; j ++)
		{		
			x[i] += aa(i, j) * b[j];			
		}
	}
    
	return (1);	//運行成功，正常返回
}

//病態方程組的求解
template <class _Ty>
int LE_IllConditionedEquation(matrix<_Ty> & a, valarray<_Ty> & b, _Ty eps, valarray<_Ty> & x)
{
	int i(60), j, kk;
    _Ty q, qq;

    int n = a.GetColNum();			//方程組階數

	matrix<_Ty> p(n, n);
	matrix<_Ty> ee(n, 1);
	matrix<_Ty> xx(n, 1);
	valarray<double> rr(n);
	    
	for(int k = 0; k < n; k++)
		for(j = 0; j < n; j++)
			p(k, j) = a(k, j);
	
	for(k = 0; k < n; k++) x[k] = b[k];
    kk = LE_TotalChoiceGauss(p, x);
	
	for(int w = 0; w < n; w++) xx(w, 0) = x[w];  
	
	if(kk == 0)
	{ 
		for( w = 0; w < n; w++) x[w] = xx(w, 0);
	    
		return 0;
	}
    
	q = 1.0 + eps;
    while(q > eps || FloatEqual(q, eps))
	{
		if(i == 0)
		{ 
			for( w = 0; w < n; w++) x[w] = xx(w, 0);
			
			return 0;
		}
        i--;        
		
		MatrixMultiply(ee, a, xx);
		for( w = 0; w < n; w++) x[w] = xx(w, 0);        
			
		for(k = 0; k < n; k++) rr[k] = b[k] - ee(k, 0);

        for(k = 0; k < n; k++)
			for(j = 0; j < n; j++) p(k, j) = a(k, j);
			
		kk = LE_TotalChoiceGauss(p, rr);	

        if(kk == 0)
		{ 
			for( w = 0; w < n; w++) x[w] = xx(w, 0);
	
			return 0;
		}
        
		q = 0.0;
        for(k = 0; k < n; k++)
		{
			qq = Abs(rr[k]) / (1.0 + Abs(x[k] + rr[k]));
            if(qq > q) q = qq;
		}
        
		for(k = 0; k < n; k++) xx(k, 0) = xx(k, 0) + rr[k];
		for(int w = 0; w < n; w++)  x[w] = xx(w, 0);	    
	}
		
	for( w = 0; w < n; w++) x[w] = xx(w, 0);
	
	return (1);	//運行成功，正常返回
}

#endif //_LINEAREQUATION_INL