lup.cpp

数值计算的例子和参考书籍

#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps = 1e-10;

const int MAXN = 1000;

bool lup_decomposition(double a[][MAXN], int n)
{
	int i, j, k, d, p[MAXN];

	for(i = 1; i <= n; i++) p[i] = i;
	int Psgn = 1;

	for(k = 1; k <= n; k++)
	{
		double max = 0;
		d = 0;
		for(i = k; i <= n; i++)
		{
			if(fabs(a[i][k]) > max)
			{
				max = fabs(a[i][k]);
				d = i;
			}
		}
		if(max < eps) return false;//singular matrix
		if(d != k)
		{
			Psgn = -Psgn;
			swap(p[k], p[d]);
			for(i = 1; i <= n; i++) swap(a[k][i], a[d][i]);
		}
		for(i = k + 1; i <= n; i++)
		{
			a[i][k] /= a[k][k];
			for(j = k + 1; j <= n; j++) a[i][j] -= a[i][k] * a[k][j];
		}
	}
//	for(i = 1; i <= n; i++)
//	{
//		for(j = 1; j <= n; j++)
//			printf("%.2lf ", a[i][j]);
//		putchar('\n');
//	}
	//double answer = 1.0;
	//for(i = 1; i <= n; i++) answer *= a[i][i];
	//return answer * Psgn;//det of matrix
	return true;
}

int main(void)
{
	int i, j;
	while(scanf("%d", &n) != EOF)
	{
		for(i = 1; i <= n; i++) for(j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
		lup_decomposition(a, n);
	}
	return 0;
}