【HDU】3359 - Kind of a Blur(高斯消元 & 矩阵)

题目链接:这里写链接内容


这里写图片描述


高斯消元模板题。


代码如下:

#include<queue>
#include<cmath>
#include<stack>
#include<cstdio>
#include<vector>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
#define INF 0x3f3f3f3f
#define CLR(a,b) memset(a,b,sizeof(a))
#define PI acos(-1.0)
const double eps = 1e-8;
struct Matrix
{
    int w,h;
    double m[105][105];

    void init(int th,int tw)
    {
        h = th;
        w = tw;
        CLR(m,0);
    }

    void Guass()
    {
        for (int i = 0 ; i < h-1 ; i++)
        {
            //交换行(减小误差) 
            int t = i;
            for (int j = i+1 ; j < h ; j++)
            {
                if (fabs(m[j][i]) > fabs(m[t][i]))
                    t = j;
            }
            if (t != i)
            {
                for (int j = 0 ; j < w ; j++)
                    swap(m[t][j],m[i][j]);
            }

            //消元
            for (int j = i+1 ; j < h ; j++)
            {
                if (fabs(m[j][i]) < eps)
                    continue;
                double mul = m[j][i] / m[i][i];
                for (int k = i ; k < w ; k++)
                    m[j][k] -= m[i][k] * mul;
            }
        }

        //反向求解
        for (int i = h-1 ; i >= 0 ; i--)
        {
            if (fabs(m[i][i]) < eps)
                continue;
            for (int j = i+1 ; j < h ; j++)
                m[i][w-1] -= m[i][j]*m[j][w-1];
            m[i][w-1] /= m[i][i];
        }
    }
};
int main()
{
    int w,h,d;
    double a[15][15];
    bool isFirst = true;
    while (~scanf ("%d %d %d",&w,&h,&d) && h)
    {
        for (int i = 0 ; i < h ; i++)
        {
            for (int j = 0 ; j < w ; j++)
                scanf ("%lf",&a[i][j]);
        }
        if (isFirst)
            isFirst = false;
        else
            printf ("\n");
        Matrix ans;
        ans.init(h*w,h*w+1);

        //构造增广矩阵 
        for (int i = 0 ; i < h ; i++)
        {
            for (int j = 0 ; j < w ; j++)
            {
                int cnt = 0;        //计数 
                for (int k = i-d ; k <= i+d ; k++)      //枚举行 
                {
                    if (k < 0 || k >= h)
                        continue;
                    int last = d - abs(i-k);        //剩余步数 
                    for (int l = j-last ; l <= j+last ; l++)        //枚举列 
                    {
                        if (l < 0 || l >= w)
                            continue;
                        cnt++;
                        ans.m[i*w+j][k*w+l] = 1;        //对应的是i*w+j个方程 
                    }
                }
                ans.m[i*w+j][w*h] = cnt*a[i][j];
            }
        }
        ans.Guass();

//      puts("=======");
//      for (int i = 0 ; i < w*h ; i++)
//      {
//          for (int j = 0 ; j < w*h+1 ; j++)
//              printf ("%-8.2lf ",ans.m[i][j]);
//          printf ("\n");
//      }
//      puts("=======");

        //输出结果
        for (int i = 0 ; i < h ; i++)
        {
            for (int j = 0 ; j < w ; j++)
                printf ("%8.2lf",ans.m[i*w+j][ans.w-1]);
            printf ("\n");
        }
    }
    return 0;
}

猜你喜欢

转载自blog.csdn.net/wyg1997/article/details/78254392