高斯消元模板

/*
    列主元高斯消元法:每一次找该列中最大的一行来消元
    3种情况:k:消元后系数矩阵非零行的个数(即有效方程个数),var:求解变量的个数
             r:增广矩阵非0行的个数(r>=k)
    1) k<r  无解  也就是说出现了这样的行(0,0,...,a)a!=0;

    2) k=r时有解
       1. k=var  唯一解
       2. k<var  无穷多解
*/


#include <iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int maxs = 100;
float X[maxs];//解向量
int A[maxs][maxs],B[maxs];//系数矩阵和值
int freeNum;//自由变元个数
int gcd(int a,int b)
{
    if(b==0)
        return a;
    return gcd(b,a%b);
}
//最小公倍数
int lcm(int a,int b)
{
    return a*b/gcd(a,b);
}
/*
    A是系数矩阵,equ是方程的个数,var是变量的个数
    求AX=B
    注意角标都是从1开始
    -1代表无解,0代表无穷多解,1代表唯一解
*/
int gaosi(int equ,int var)
{
    int k,col;//当前处理的行和列
    for(k=1,col=1;k<=equ&&col<=var;k++,col++)
    {
        int max_r=k;
        int maxValue=abs(A[k][col]);
        //往下找到该列中最大行
        for(int i=k+1;i<=equ;i++)
            if(abs(A[i][col])>maxValue)
            {
                maxValue=abs(A[i][col]);
                max_r=i;
            }
        if(max_r!=k)
        {
            //交换两行
            for(int i=1;i<=var;i++)
                swap(A[k][i],A[max_r][i]);
            swap(B[k],B[max_r]);
        }
        //如果值为0,说明该行以下这一列也全为0,那么继续处理该行的下一列进行消元
        /*  系数矩阵是这种情况,那么对第二行的第三个元素进行消元
            3 2 5
            0 0 3
            0 0 2
        */
        if(A[k][col]==0)
        {
            k--;continue;
        }
        //开始进行消元
        for(int i=k+1;i<=equ;i++)
        {
            if(A[i][col]!=0)
            {
                int LCM = lcm(abs(A[k][col]),abs(A[i][col]));
                int tk = LCM/abs(A[k][col]);
                int ti = LCM/abs(A[i][col]);
                if(A[k][col]*A[i][col]<0)
                    ti=-ti;//异号时相加
                for(int j=col;j<=var;j++)
                    A[i][j]=A[i][j]*ti-A[k][j]*tk;
                B[i]=B[i]*ti-B[k]*tk;
            }
        }
    }
    //因为每进行一次消元,那么系数矩阵的非0行就加1,退出循环时多加了一次
    k=k-1;//系数矩阵非0行的个数(k<=var)

    //无解
    for(int i=k+1;i<=equ;i++)
        if(B[i]!=0)
            return -1;
    //唯一解
    if(k==var)
    {
        for(int i=k;i>=1;i--)
        {
            float temp = B[i]*1.0;
            for(int j=var;j>i;j--)
                if(A[i][j]!=0)
                    temp=temp-X[j]*A[i][j];
            X[i]=temp/A[i][i];
        }
        return 1;
    }
    if(k<var)
    {
        freeNum = var-k;
        return 0;
    }
    return 1;
}
int main()
{
    freopen("in.txt","r",stdin);
    int equ,var;
    while(scanf("%d%d",&equ,&var)!=EOF)
    {
        for(int i=1;i<=equ;i++)
            for(int j=1;j<=var;j++)
                scanf("%d",&A[i][j]);
        for(int i=1;i<=equ;i++)
            scanf("%d",&B[i]);
        int flag = gaosi(equ,var);
        if(flag==-1)
            printf("无解\n");
        if(flag==0)
            printf("无穷多解\n自由变元个数: %d\n",freeNum);
        if(flag==1)
        {
            printf("唯一解\n");
            for(int i=1;i<=var;i++)
                printf("x%d = %f\n",i,X[i]);
        }
    }
    return 0;
}

 

posted on 2016-08-19 20:19  wastonl  阅读(178)  评论(0编辑  收藏  举报