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

#include

#include

using namespace std;