【算法设计与分析基础】高斯消元

参考《算法导论》第七部分算法问题选编中的第28章矩阵运算。

#include<iostream>
using namespace std;/*** 一些矩阵例子*//*
double Value[9][9]={{0,0,0,-85,0,59,0,0,72},{97,-137,122,81,97,55,0,132,0},{0,-69,145,70,82,0,66,0,88},{101,0,-103,76,0,0,-58,-143,0},{0,-100,0,-171,-93,0,89,0,0},{0,0,59,154,0,0,0,-50,0},{0,111,0,-161,-80,0,85,-173,0},{112,0,250,111,-167,208,96,0,116},{62,0,86,-111,-90,-228,132,-61,81}
},
b[9]={1,1,1,1,1,1,1,1,1};
*/
double Value[9][9]={{0,0,0,-85,0,59,0,0,72},{1,-137,122,81,97,55,0,132,0},{0,-69,145,70,82,0,66,0,88},{101,0,0,76,0,0,-58,-143,0},{0,-100,0,-171,-93,0,89,0,0},{0,0,59,154,0,0,0,-50,0},{0,111,0,-161,-80,0,85,0,0},{112,0,250,0,-167,208,96,0,116},{62,0,86,-111,0,0,132,-61,81}
},
b[9]={1,1,1,1,1,1,1,1,1};double Example[4][4]={{2,3,1,5},{6,13,5,19},{2,19,10,23},{4,10,11,31},
};
double Examplee[3][3]={{1,2,0},{3,4,4},{5,6,3},
},
Ex[3]={3,7,8};double Exampleee[4][4]={{2,0,2,0.6},{3,3,4,-2},{5,5,4,2},{-1,-2,3.4,-1},
},
Ex1[4]={1,1,1,1};class Gauss{
private:double **A;//原矩阵double *x;//线性方程组的解double *y;//采用正向替换，对Ly=Pb求解ydouble *b;//线性方程组最右边一列的常数double **L;//LU分解的L数组double **U;//LU分解的U数组double **l;//L的逆矩阵double **u;//U的逆矩阵double **a;//A的逆矩阵 = l * udouble **p;//根据数组P求出置换矩阵int *P;//行数排列int num;
public:Gauss(int n);//构造函数，构造一个 n x n 的矩阵~Gauss();//析构函数void setA(int x,int y,double v);//给数组A的[x][y]赋值为vvoid setB(int x,double v);//给数组B的[x]赋值为vvoid Add(int r1,int r2);//第r1行加上r2行void LU();// 计算 L 和 U 矩阵void LUP();// 计算出 L 和 U 矩阵 以及 排列行数Pvoid solve_p();//计算 置换矩阵p 并 修改原始矩阵的顺序void LUP_SOLVE();//通过 L矩阵 和 U矩阵 解决线性方程组问题void U_();//计算 U 的逆矩阵void L_();//计算 L 的逆矩阵void A_();//计算 A 的逆矩阵void A_Next();//计算 原始矩阵 的逆矩阵 即 A的逆矩阵乘以置换矩阵void displayB();//输出 线性方程组 最右边一列常数void displayL();//输出 L 矩阵void displayU();//输出 U 矩阵void displayDiagonal();//输出对角线void displayx();//输出线性方程组的解void displayA();//输出原矩阵void displayy();//输出 y 矩阵void displayu();//输出 U 的逆矩阵void displayl();//输出 L 的逆矩阵void displaya();//输出 A 的逆矩阵void displayP();//输出 排列行数Pvoid displayp();//输出 置换矩阵p
};Gauss::Gauss(int n){//构造函数，构造一个 n x n 的矩阵int i;num=n;A=new double*[n];for(i=-1;++i<n;){A[i]=new double[n];memset(A[i],0,n*sizeof(A[i]));}x=new double[n];b=new double[n];}Gauss::~Gauss(){//析构函数/*  int i;for(i=-1;++i<num;){delete[]A[i];delete[]L[i];delete[]U[i];}delete[]y;delete[]A;delete[]x;delete[]b;delete[]U;*/}void Gauss::setA(int x,int y,double v){//给数组A的[x][y]赋值为vif( (x<0 || x>=num) || (y<0 || y>=num)  )return;A[x][y]=v;}void Gauss::setB(int x,double v){//给数组B的[x]赋值为vb[x]=v;}void Gauss::Add(int r1,int r2){//r1行加上r2行int i=-1;for(;++i<num;){A[r1][i]+=A[r2][i];}//forb[r1]+=b[r2];}void Gauss::LU(){// 计算 L 和 U 矩阵//让L 和 U 成为 n x n 的矩阵L=new double*[num];U=new double*[num];int i=-1,j,k;for(;++i<num;){L[i]=new double[num];U[i]=new double[num];//L 和 U 全部初始化为0memset(L[i],0,num*sizeof(double));memset(U[i],0,num*sizeof(double));//L的对角线全是1L[i][i]=1;}for(k=-1;++k<num;){U[k][k]=A[k][k];for(i=k;++i<num;){L[i][k]=A[i][k]/U[k][k];U[k][i]=A[k][i];}for(i=k;++i<num;){for(j=k;++j<num;){A[i][j]=A[i][j]-L[i][k] * U[k][j];}}}  }void Gauss::LUP(){//计算出 L 和 U 矩阵 以及 排列行数PP=new int[num];double **_A=new double*[num];int i=-1,j,k,p,_k;double temp;/*** 给P赋值为0 到 num-1*/for(;++i<num;){_A[i]=new double[num];P[i]=i;for(j=-1;++j<num;){_A[i][j]=A[i][j];}}//forfor(k=-1;++k<num;){p=0;for(i=k;i<num;i++){if(abs(_A[i][k])>p){p=abs(_A[i][k]);_k=i;}}if(p==0){cout<<"singular matrix"<<endl;exit(0);}temp=P[k];P[k]=P[_k];P[_k]=temp;for(i=-1;++i<num;){temp=_A[k][i];_A[k][i]=_A[_k][i];_A[_k][i]=temp;}for(i=k;++i<num;){_A[i][k]/=_A[k][k];for(j=k;++j<num;){_A[i][j]-=_A[i][k]*_A[k][j];}}}}void Gauss::solve_p(){//计算 置换矩阵p 并 修改原始矩阵的顺序int i,j;double **k=new double*[num];p=new double*[num];for(i=-1;++i<num;){p[i]=new double[num];k[i]=new double[num];memset(p[i],0,sizeof(double)*num);for(j=-1;++j<num;){k[i][j]=A[i][j];}}for(i=-1;++i<num;){p[i][P[i]]=1;}for(i=-1;++i<num;){for(j=-1;++j<num;){A[i][j]=k[P[i]][j];}}}void Gauss::LUP_SOLVE(){//通过 L矩阵 和 U矩阵 解决线性方程组问题x=new double[num];y=new double[num];memset(x,0,sizeof(double)*num);memset(y,0,sizeof(double)*num);int i=-1,j;double sum;for(;++i<num;){sum=0;for(j=-1;++j<i;){sum+=L[i][j]*y[j];}y[i]=b[i]-sum;}for(i=num;--i>=0;){sum=0;for(j=num;--j>i;){sum+=U[i][j]*x[j];}x[i]=(y[i]-sum)/U[i][j];}}void Gauss::U_(){//计算 U 的逆矩阵int i,j,k;double sum;u=new double*[num];for(i=-1;++i<num;){u[i]=new double[num];memset(u[i],0,sizeof(double)*num);}for(i=-1;++i<num;){u[i][i]=1/U[i][i];for(k=i;--k>=0;){sum=0;for(j=k;++j<=i;){sum+=U[k][j]*u[j][i];u[k][i]=-sum/U[k][k];}}}}void Gauss::L_(){//计算 L 的逆矩阵int i,j,k;l=new double*[num];for(i=-1;++i<num;){l[i]=new double[num];memset(l[i],0,sizeof(double)*num);}for(i=-1;++i<num;){l[i][i]=1;//对角线依次赋值为1for(k=i;++k<num;){for(j=i-1;++j<=k-1;){l[k][i]=l[k][i]-L[k][j]*l[j][i];}}}}void Gauss::A_(){//计算 A 的逆矩阵int i,j,k;double sum;a=new double*[num];for(i=-1;++i<num;){a[i]=new double[num];memset(a[i],0,sizeof(double)*num);}for(i=-1;++i<num;){for(j=-1;++j<num;){sum=0;for(k=-1;++k<num;){sum+=u[i][k]*l[k][j];}a[i][j]=sum;}}}void Gauss::A_Next(){//计算 原始矩阵 的逆矩阵 即 A的逆矩阵乘以置换矩阵int i,j,k;double sum;double **aa=new double*[num];for(i=-1;++i<num;){aa[i]=new double[num];for(j=-1;++j<num;){aa[i][j]=a[i][j];}}for(i=-1;++i<num;){for(j=-1;++j<num;){for(j=-1;++j<num;){sum=0;for(k=-1;++k<num;){sum+=aa[i][k]*p[k][j];}a[i][j]=sum;}}}}void Gauss::displayB(){//输出 线性方程组 最右边一列常数int i=-1;for(;++i<num-1;){cout<<b[i]<<' ';}cout<<b[i]<<endl;}void Gauss::displayL(){//输出 L 矩阵int i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<L[i][j]<<' ';}cout<<L[i][j]<<endl;}cout<<endl;}void Gauss::displayU(){//输出 U 矩阵int i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<U[i][j]<<' ';}cout<<U[i][j]<<endl;}cout<<endl;}void Gauss::displayDiagonal(){//输出对角线int i=-1;for(;++i<num-1;){cout<<A[i][i]<<" ";}//forcout<<A[i][i]<<endl;}void Gauss::displayx(){//输出线性方程组的解int i=-1;for(;++i<num;){cout<<'x'<<i+1<<" = "<<x[i]<<endl;}//for}void Gauss::displayA(){//输出原矩阵int i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<A[i][j]<<',';}cout<<A[i][j]<<endl;}cout<<endl;}void Gauss::displayy(){//输出 y 矩阵int i=-1;for(;++i<num;){cout<<'y'<<i+1<<" = "<<y[i]<<endl;}//for}void Gauss::displayu(){//输出 U 的逆矩阵int i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<u[i][j]<<' ';}cout<<u[i][j]<<endl;}cout<<endl;}void Gauss::displayl(){//输出 L 的逆矩阵int i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<l[i][j]<<' ';}cout<<l[i][j]<<endl;}cout<<endl;}void Gauss::displaya(){//输出 A 的逆矩阵int i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<a[i][j]<<' ';}cout<<a[i][j]<<endl;}cout<<endl;}void Gauss::displayP(){//输出 排列行数Pint i=-1;for(;++i<num;){cout<<'p'<<i+1<<" = "<<P[i]<<endl;}//for}void Gauss::displayp(){//输出 置换矩阵pint i=-1,j;for(;++i<num;){for(j=-1;++j<num-1;){cout<<p[i][j]<<' ';}cout<<p[i][j]<<endl;}cout<<endl;}int main(){Gauss G(9);int i=-1,j;/*** 初始化线性方程组的矩阵*/for(;++i<9;){G.setB(i,b[i]);for(j=-1;++j<9;){G.setA(i,j,Value[i][j]);}}
/*  Gauss G(3);int i=-1,j;for(;++i<3;){G.setB(i,Ex[i]);for(j=-1;++j<3;){G.setA(i,j,Examplee[i][j]);}}*/
/*  Gauss G(4);int i=-1,j;for(;++i<4;){G.setB(i,Ex1[i]);for(j=-1;++j<4;){G.setA(i,j,Exampleee[i][j]);}}*//*** 一旦主元选取为0，就会有除以0的问题，主元处于矩阵U的对角线上。** 由于调用LU分解的时候有 除以0的现象* 所以把矩阵中为 0 的数字 去除* 方法是 这一行 加上 另一行*//*    G.Add(0,1);G.Add(0,3);G.Add(1,0);G.Add(2,1);G.Add(3,2);G.Add(4,3);G.Add(5,4);G.Add(6,5);G.Add(7,6);G.Add(8,7);*//*** 上面是靠初始行变换来创造调用 LU分解算法 的条件，这样的话无法求逆矩阵* 这里是计算置换矩阵P来创造LU分解算法的条件*/cout<<"原始矩阵"<<endl;G.displayA();G.LUP();//计算重排列行数Pcout<<"输出排列行数"<<endl;G.displayP();cout<<"输出置换矩阵P"<<endl;G.solve_p();//计算置换矩阵pG.displayp();cout<<endl;cout<<"b矩阵"<<endl;G.displayB();cout<<endl;cout<<"A矩阵"<<endl;G.displayA();G.LU();//计算LU矩阵G.LUP_SOLVE();//通过L矩阵和U矩阵计算线性方程组的解cout<<"L矩阵"<<endl;G.displayL();cout<<"U矩阵"<<endl;G.displayU();cout<<"y矩阵"<<endl;G.displayy();cout<<"线性方程组的解"<<endl;G.displayx();G.U_();//计算U的逆矩阵cout<<"U矩阵的逆矩阵"<<endl;G.displayu();G.L_();//计算L的逆矩阵cout<<"L矩阵的逆矩阵"<<endl;G.displayl();G.A_();//计算A的逆矩阵cout<<"A矩阵的逆矩阵"<<endl;G.displaya();G.A_Next();//计算原始矩阵的逆矩阵 = A的逆矩阵乘以置换矩阵pcout<<"原始矩阵的逆矩阵"<<endl;G.displaya();cout<<endl;return 0;
}

高斯消元LU分解PPT：http://download.csdn.net/detail/u013580497/8892729

【算法设计与分析基础】高斯消元相关推荐

算法设计与分析基础-笔记-上
算法设计与分析基础绪论什么是算法一系列解决问题的明确指令,对于符合一定规范的输入,能够在有限的时间内获得要求的输出. 例子:最大公约数:俩个不全为0 的非负整数 m m m和 n n n的最大公 ...
计算机算法设计与分析读后感,算法设计与分析基础经典读后感有感
<算法设计与分析基础>是一本由Anany levitin著作,清华大学出版社出版的胶版纸图书,本书定价:49.00元,页数:409,特精心从网络上整理的一些读者的读后感,希望对大家能有帮助 ...
第一章算法设计与分析基础知识
系列文章目录第一章算法设计与分析基础知识第二章算法的分治策略第三章算法的动态规划第四章算法的贪心法 -- @[TOC](这里写目录标题) # 一级目录 ## 二级目录 ### 三级目录 ...
算法设计与分析基础知识
一.算法设计基础算法是(algorithm)是对特定问题求解步骤的一种描述,是指令的有限序列. 算法的五个特性: 输入:一个算法可以有零个或多个输入. 输出:一个算法有一个输出或多个输出. 有穷性( ...
算法设计与分析基础知识点
前言:全文参考徐承志老师的PPT 适合期末复习,查缺补漏,有缺漏或错误欢迎指正,后面的第九章内容之后会继续补充. 目录一.算法基础概念二.算法分析基础 1.概念 2.算法设计的一般过程 3.时间复 ...
算法设计与分析基础第三版
课后题答案一.算法级基础知识 1.算法的基本概念解决问题的确定方法和有限步骤称为算法,对于计算机科学来说,算法指的是对特定问题的求解步骤的一种描述,是若干条指令的有穷序列.并有以下特性:输入.输出 ...
算法设计与分析基础第一章谜题
习题1.1 10.b 欧几里得游戏一开始,板上写有两个不相等的正整数,两个玩家交替写数字,每一次,当前玩家都必须在板上写出任意两个板上数字的差,而且这两个数字必须是新的,也就是说,不能与板上任何一个 ...
算法设计与分析基础第六章谜题
习题6.1 9.数字填空给定n个不同的整数以及一个包含n个空格的序列,每个空格之间事先给定有不等(>或<)符号,请设计一个算法,将n个整数填入这n个空格中并满足不等式约束.例如,数4,6 ...
算法设计与分析基础第五章谜题
习题5.1 11.Tromino谜题 Tromino是一个由棋盘上的三个1×1方块组成的L型骨牌.我们的问题是,如何用Tromino覆盖一个缺少了一个方块的2n×2n棋盘.除了这个缺失的方块,Trom ...

【算法设计与分析基础】高斯消元

【算法设计与分析基础】高斯消元相关推荐

最新文章

热门文章