3.cholesky分解法.cpp

cholesky算法

#include "iostream.h"
#include "math.h"
#define n 3

void main()
{
	int i , j , k ;
	double sum ;
	double a[n][n] , b[n]  ;
	double x[n], y[n] ;

	cout << "输入对称正定矩阵a[n][n]:" << endl;
	for (i = 0 ; i < n ; i++)
	{
		for(j = 0 ; j < n ; j++)
		{
			cin >> a[i][j] ;
		}
	}

	cout << "输入b[n]:" << endl ;
	for(i = 0 ; i < n ; i++)
	{
		cin >> b[i] ;
	}  //  输入部分

	
	for(k = 0 ; k < n ; k++)
	{
		a[k][k] = sqrt(a[k][k]) ;
		for(i = k + 1 ; i < n ; i++)
		{
			a[i][k] = a[i][k] / a[k][k] ;
		}
		for(j = k + 1 ; j < n ; j++)
		{
			for(i = j ; i < n ; i++)
			{
				a[i][j] = a[i][j] - a[i][k] * a[j][k] ;
			}
		}
	}    //   求L*sqrt(D)

	for(i = 0 ; i < n ; i++)
	{
		for(j = i ; j < n ; j++)
		{
			a[i][j] = a[j][i] ;
		}
	}     //    得到LU矩阵
   
	y[0] = b[0] ;
	for(i = 1 ; i < n ; i++)
	{
		sum = 0.0 ;
		for( j = 0 ; j <= i-1 ; j++)
		{
			sum = sum + a[i][j] * y[j] ;
		}
		y[i] = b[i] - sum ;
	}
	cout << "输出y1到yn:" << endl ;
	for(i = 0 ; i < n ; i++)
	{
		cout << y[i] << ' ' ;
	}  //  按LU分解法求y[i]并输出
	cout << endl;

	x[n-1] = y[n-1] / a[n-1][n-1] ;
	for (i = n-2 ; i >= 0 ; i--)
	{
		sum = 0.0 ;
		for( j = i + 1 ; j < n ; j++)
		{
			sum = sum + a[i][j] * x[j] ;
		}
		x[i] = (y[i] - sum) / a[i][i] ;
	}
	cout << "输出x1到xn:" << endl ;
	for( i = 0 ; i < n ; i++)
	{
		cout << x[i] << ' ' ;
	}  //  按x[i]并输出
	cout << endl ;
}



/*  实验矩阵
       5  2  -4
       2  1  -2
	   -4 -2  5
	   */