PLU Decomposition

PLU分解的好处是，可以将Ax=b的矩阵，转换成Ly=b, Ux = y

的形式，当我们改变系数矩阵b时，此时由于矩阵L和U均是固定

的，所以总能高效的求出矩阵的解。

// LU.cpp : Defines the entry point for the console application.
//
/************************************************
*   Author:     JohnsonDu
*   From:       Institute of Computing Technology
*               University of Chinese Academy of Science
*   Time:       2014-10-7
*   Content:    PLU decomposition
*************************************************/

#include "stdafx.h"

#define MAXN 5005
#define eps 1e-9         // 精度
int n, m;
double mat[MAXN][MAXN];  // 输入矩阵
double matL[MAXN][MAXN]; // 矩阵L
double matU[MAXN][MAXN]; // 矩阵U
int matP[MAXN][MAXN];    // 矩阵P
int seq[MAXN];           // 记录行变换
//double vecB[MAXN];
//double vecY[MAXN];
//double vecX[MAXN];
//double matPb[MAXN];

void menu()
{
	printf("----------------PLU Factorization---------------\n");
	printf("|         Please follow the instruction        |\n");
	printf("|         to determine the LU decomposition    |\n");
	printf("|         PA = LU                              |\n");
	printf("------------------------------------------------\n\n");

}

void initLMatrix()
{
	memset(matU, 0, sizeof(matU));
	memset(matL, 0, sizeof(matL));
	memset(matP, 0, sizeof(matP));
}

void padLMatrix()
{
	for(int i = 0; i < n; i ++)
		matL[i][i] = 1.0;
}

inline double Abs(double x)
{
	return x < 0 ? -x : x;
}

void displayLU()
{
	// 输出矩阵L
	printf("\n----------------------\n");
	printf("Matrix L follows: \n");
	for(int i = 0; i < n; i ++)
	{
		for(int j = 0; j < (n < m ? n : m); j ++)
			printf("%.3f  ", matL[i][j]);
		printf("\n");
	}

	// 输出矩阵U
	printf("\nMatrix U follows: \n");
	for(int i = 0; i < (n < m ? n : m); i ++)
	{
		for(int j = 0; j < m; j ++)
			printf("%.3f  ", matU[i][j]);
		printf("\n");
	}

	// 输出矩阵P
	printf("\nMatrix P follows: \n");
	for(int i = 0; i < n; i ++)
	{
		for(int j = 0; j < n; j ++)
			printf("%d  ", matP[i][j]);
		printf("\n");
	}
	printf("----------------------\n");
}

/*
// 输出LU的过程及最终解
void displaySolution()
{
	// 输出矩阵Pb
	printf("\nMatrix Pb follows: \n");
	for(int i = 0; i < n; i ++)
	{
		printf("%.3f\n", matPb[i]);
	}

	// 输出向量y
	printf("\nVector Y follows: \n");
	for(int i = 0; i < n; i ++)
	{
		printf("%.3f\n", vecY[i]);
	}
	printf("\n");

	// 输出解向量x
	printf("\Vector X follows: \n");
	for(int i = 0; i < n; i ++)
	{
		printf("%.3f\n", vecX[i]);
	}
	printf("\n");
}
*/

// 交换元素
inline void swap(int &a, int &b)
{
	int t = a;
	a = b;
	b = t;
}

// 高斯消元部分
void gauss()
{
	int i;
	int col;
	int max_r;  

	col = 0;    //处理的当前列

	// 从第一行开始进行消元
	// k为处理的当前行
	for(int k = 0; k < n && col < min(n, m); k ++, col ++)
	{
		// 寻找当前col列的绝对值最大值
		max_r = k;
		for(i = k + 1; i < n; i ++)
			if(Abs(mat[i][col]) > Abs(mat[max_r][col]))
				max_r = i;

		// 进行行交换
		if(max_r != k)
		{
			for(int j = col; j < m; j ++)
				swap(mat[k][j], mat[max_r][j]);
			swap(seq[k], seq[max_r]);
			for(int j = 0; j < n; j ++)
				swap(matL[k][j], matL[max_r][j]);
		}

		// 当前主元为零， 继续
		if(Abs(mat[k][col]) < eps){
			continue;
		}

		// 消元部分，并获得L矩阵
		for(int i = k + 1; i < n; i ++)
		{

			double t = mat[i][col] / mat[k][col];
			matL[i][col] = t;
			for(int j = col; j < m; j ++)
				mat[i][j] -= t * mat[k][j];
		}

	}

	// 为矩阵U进行赋值
	for(int i = 0; i < n; i ++)
		for(int j = 0; j < m; j ++)
			matU[i][j] = mat[i][j];

	// 生成矩阵P
	for(int i = 0; i < n; i ++)  matP[i][seq[i]] = 1.0;

	// 为矩阵L添加对角线元素
	padLMatrix();
}

/*
// 计算Pb的值
void calcPb()
{
	for(int i = 0; i < n; i ++)
		matPb[i] = 0.0;
	//cout << "-----------" << endl;
	for(int i = 0; i < n; i ++)
	{
		double t = 0.0;
		for(int j = 0; j < n; j ++)
		{
			t = t + 1.0 * matP[i][j] * vecB[j];
			//cout << t << endl;
			//cout << matP[i][j] * vecB[j] << "---" << endl;
		}
		matPb[i] = t;
		//cout << matPb[i] << endl;
	}
}

// 计算Ly = Pb, y向量
void calcY()
{
	vecY[0] = matPb[0];
	for(int i = 1; i < n; i ++)
	{
		double t = 0.0;
		for(int j = 0; j < i; j ++)
			 t += vecY[j] * matL[i][j];
		vecY[i] = matPb[i] - t;
	}
}

// 计算Ux = y, y向量
void calcX()
{
	vecX[n-1] = vecY[n-1] / matU[n-1][n-1];
	for(int i = n-2; i >= 0; i --)
	{
		double t = 0.0;
		for(int j = n-1; j > i; j --)
			 t += vecX[j] * matU[i][j];
		vecX[i] = (vecY[i] - t) / matU[i][i];
	}
}
*/

int _tmain(int argc, _TCHAR* argv[])
{
	menu();
	while(true)
	{
		printf("Please input the matrix's dimension n & m: ");
		// 输入矩阵的行n和列m
		scanf("%d%d", &n, &m);
		printf("Please input the matrix: \n");

		// 输入矩阵
		for(int i = 0; i < n; i ++)
		{
			for(int j = 0; j < m; j ++)
				cin >> mat[i][j];
			seq[i] = i;
		}

		// 初始化为0
		initLMatrix();

		// 高斯消元
		gauss();

		// 输出P, L, U矩阵
		displayLU();
		system("pause");
		system("cls");
		menu();

/*
        //此处是输入b，求取x, y 和 pb
		while(true){
			printf("please input vector b(whose length equals to %d): \n", n);
			for(int i = 0; i < n; i ++) cin >> vecB[i];
			calcPb();
			calcY();
			calcX();
			displaySolution();
		}
*/
	}

	return 0;
}

其中stdafx.h的头文件：

// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//

#pragma once
#define  _CRT_SECURE_NO_WARNINGS

#include "targetver.h"

#include <stdio.h>
#include <tchar.h>
#include <iostream>
using namespace std;