zoukankan      html  css  js  c++  java
  • 高斯消元模板

    /*
        列主元高斯消元法:每一次找该列中最大的一行来消元
        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("无解
    ");
            if(flag==0)
                printf("无穷多解
    自由变元个数: %d
    ",freeNum);
            if(flag==1)
            {
                printf("唯一解
    ");
                for(int i=1;i<=var;i++)
                    printf("x%d = %f
    ",i,X[i]);
            }
        }
        return 0;
    }
    
  • 相关阅读:
    AcWing 852. spfa判断负环 边权可能为负数。
    面试题 02.01. 移除重复节点
    1114. 按序打印
    剑指 Offer 38. 字符串的排列
    557. 反转字符串中的单词 III
    645. 错误的集合
    面试题 05.03. 翻转数位
    1356. 根据数字二进制下 1 的数目排序
    748. 最短完整词
    剑指 Offer 32
  • 原文地址:https://www.cnblogs.com/wt20/p/5788899.html
Copyright © 2011-2022 走看看