Advertisement
Jade39

multigrid

Apr 17th, 2014
58
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 8.50 KB | None | 0 0
  1. #include <iostream>
  2. #include <fstream>
  3. #include <cmath>
  4. #include <conio.h>
  5.  
  6. using namespace std;
  7.  
  8. typedef double** vector;
  9.  
  10. const double Star[3][3] = { 0.25, 0.5, 0.25,
  11.                             0.5,  1,   0.5,
  12.                             0.25, 0.5, 0.25 };
  13. int NSize = 1023; // размер блочной матрицы
  14. int NCell = 1023; // размер блока
  15. int NSmoothIters = 1; // количество сглаживающих итераций
  16. int MaxDeep = 2; // глубина рекурсии
  17. double epsS = 1e-8;
  18. double epsM = 1e-8;
  19. vector f, tmp;
  20. vector *v, *r;
  21. fstream fout("out.txt",ios::out);
  22.  
  23. double _abs(double a)
  24. {
  25.     return (a > 0 ? a : -a);
  26. }
  27.  
  28. void Set(vector a, vector b, int nS, int nC)
  29. {
  30.     int i,j;
  31.     for (i=0;i<nS;i++)
  32.         for (j=0;j<nC;j++)
  33.             a[i][j]=b[i][j];
  34. }
  35.  
  36. void Allocate()
  37. {
  38.     int i,j,k;
  39.     int nS = NSize, nC = NCell;
  40.     v = new vector[MaxDeep + 1];
  41.     r = new vector[MaxDeep + 1];
  42.     f = new double* [NSize];
  43.     tmp = new double* [NSize];
  44.     for (i = 0; i <= MaxDeep; i++)
  45.     {  
  46.         v[i] = new double* [nS];
  47.         r[i] = new double* [nS];
  48.         for (j = 0; j < nS; j++)
  49.         {  
  50.             v[i][j] = new double [nC];
  51.             for (k = 0; k < nC; k++)
  52.                 v[i][j][k] = 0;
  53.             r[i][j] = new double [nC];
  54.             for (k = 0; k < nC; k++)
  55.                 r[i][j][k] = 0;
  56.         }
  57.        
  58.         nS = nS / 2;
  59.         nC = nC / 2;
  60.     }
  61.     for (i = 0; i < NSize; i++)
  62.     {
  63.         f[i] = new double [NCell];
  64.         tmp[i] = new double [NCell];
  65.     }
  66. }
  67.  
  68. void Restriction(vector a, vector b, int nS, int nC) // блочный вектор a[nS][nC] переводим
  69. // в блочный вектор b[nS/2][nC/2];
  70. {
  71.     int i,j,k,l,ii,jj, nS2 = nS/2, nC2 = nC/2;
  72.     for (i = 0; i < nS2; i++)
  73.         for (j = 0; j < nC2; j++)
  74.         {
  75.             ii = 2 * i + 1;
  76.             jj = 2 * j + 1;
  77.             b[i][j] = 0;
  78.             for (k = -1; k < 2; k++)
  79.                 for (l = -1; l < 2; l++)
  80.                     if ((ii + k < nS) && (ii + k >= 0) && (jj + l < nC) && (jj + l >= 0))
  81.                         b[i][j] += a[ii + k][jj + l] * Star [k + 1][l + 1];
  82.         }
  83. }
  84.  
  85. void Prolongation(vector a, vector b, int nS, int nC) //блочный вектор a[nS][nC] переводим
  86. // в блочный вектор b[2*nS + 1][2*nC + 1];
  87. {
  88.     int i,j,k,l,ii,jj, nS2 = 2*nS + 1, nC2 = 2*nC + 1;
  89.     for (i = 0; i < nS2; i++)
  90.         for (j = 0; j < nC2; j++)
  91.             b[i][j] = 0;
  92.     for (i = 0; i < nS; i++)
  93.         for (j = 0; j < nC; j++)
  94.         {
  95.             ii = 2 * i + 1;
  96.             jj = 2 * j + 1;
  97.             b[ii][jj] = 0;
  98.             for (k = -1; k < 2; k++)
  99.                 for (l = -1; l < 2; l++)
  100.                     if ((ii + k < nS2) && (ii + k >= 0) && (jj + l < nC2) && (jj + l >= 0))
  101.                         b[ii + k][jj + l] += a[i][j] * Star [k + 1][l + 1];
  102.         }
  103. }
  104.  
  105. void ZeidelIteration1(vector u, vector f, int nS, int nC)
  106. {
  107.     int i, j;
  108.     for (i = 0; i < nS; i++)
  109.         for (j = 0; j < nC; j++)
  110.         {
  111.             u[i][j] = f[i][j];
  112.             if (i != 0) u[i][j] += u[i-1][j];
  113.             if (i != nS - 1) u[i][j] += u[i+1][j];
  114.             if (j != 0) u[i][j] += u[i][j-1];
  115.             if (j != nC - 1) u[i][j] += u[i][j+1];
  116.             u[i][j] = u[i][j] / 4;
  117.         }
  118. }
  119.  
  120. void ZeidelIteration2(vector u, vector f, int nS, int nC)
  121. {
  122.     int i, j;
  123.     for (i = nS-1; i >= 0; i--)
  124.         for (j = nC - 1; j >= 0; j--)
  125.         {
  126.             u[i][j] = f[i][j];
  127.             if (i != 0) u[i][j] += u[i-1][j];
  128.             if (i != nS - 1) u[i][j] += u[i+1][j];
  129.             if (j != 0) u[i][j] += u[i][j-1];
  130.             if (j != nC - 1) u[i][j] += u[i][j+1];
  131.             u[i][j] = u[i][j] / 4;
  132.         }
  133. }
  134.  
  135. void ZeidelMethod(vector u, vector f, double eps, int nS, int nC)
  136. {
  137.     int i,j;
  138.     double max, tmp;
  139.     while (true)
  140.     {
  141.         max = 0;
  142.         for (i = 0; i < nS; i++)
  143.             for (j = 0; j < nC; j++)
  144.             {
  145.                 tmp = f[i][j];
  146.                 if (i != 0) tmp += u[i-1][j];
  147.                 if (i != nS - 1) tmp += u[i+1][j];
  148.                 if (j != 0) tmp += u[i][j-1];
  149.                 if (j != nC - 1) tmp += u[i][j+1];
  150.                 tmp = tmp / 4;
  151.                 if (_abs(tmp - u[i][j]) > max) max = _abs(tmp - u[i][j]);
  152.                 u[i][j] = tmp;
  153.             }
  154.             if (max < eps)
  155.                 return;
  156.     }
  157. }
  158.  
  159. void SmoothingIteration(vector u, vector f, int nS, int nC) // сглаживающая итерация
  160. {
  161. //  ZeidelIteration(u,f,nS,nC);
  162. }
  163.  
  164. void RoughSolution(vector u, vector f, double eps, int nS, int nC) // решение на грубой сетке
  165. {
  166.     ZeidelMethod(u,f,eps,nS,nC);
  167. }
  168.  
  169. void Remainder(vector b, vector u, vector res, int nS, int nC) // res = b - Au
  170. {
  171.     int i, j;
  172.     for (i = 0; i < nS; i++)
  173.         for (j = 0; j < nC; j++)
  174.         {
  175.             res[i][j] = b[i][j] - 4 * u[i][j];
  176.             if (i != 0) res[i][j] += u[i-1][j];
  177.             if (i != nS - 1) res[i][j] += u[i+1][j];
  178.             if (j != 0) res[i][j] += u[i][j-1];
  179.             if (j != nC - 1) res[i][j] += u[i][j+1];
  180.         }
  181. }
  182.  
  183. void GetAuDecrementToRes(vector u, vector res, int nS, int nC) // res = b - Au
  184. {
  185.     int i, j;
  186.     for (i = 0; i < nS; i++)
  187.         for (j = 0; j < nC; j++)
  188.         {
  189.             res[i][j] -=  4 * u[i][j];
  190.             if (i != 0) res[i][j] += u[i-1][j];
  191.             if (i != nS - 1) res[i][j] += u[i+1][j];
  192.             if (j != 0) res[i][j] += u[i][j-1];
  193.             if (j != nC - 1) res[i][j] += u[i][j+1];
  194.         }
  195. }
  196.  
  197. void Sum(vector a, vector b, int nS, int nC) // результат в a
  198. {
  199.     for (int i = 0; i < nS; i++)
  200.         for (int j = 0; j < nC; j++)
  201.             a[i][j] += b[i][j];
  202. }
  203.  
  204. void Diff(vector a, vector b, int nS, int nC) // результат в a
  205. {
  206.     for (int i = 0; i < nS; i++)
  207.         for (int j = 0; j < nC; j++)
  208.             a[i][j] -= b[i][j];
  209. }
  210.  
  211. void PrintVector(vector a, int nS, int nC)
  212. {
  213.     int i,j;
  214.     for (i = 0; i < nS; i++)
  215.     {  
  216.         for (j = 0; j < nC; j++)
  217.         {  
  218.             fout.width(10);
  219.             fout.precision(4);
  220.             fout << a[i][j];
  221.         }
  222.         fout << endl;
  223.     }
  224. }
  225.  
  226. void Zero(vector a, int nS, int nC)
  227. {
  228.     int i,j;
  229.     for (i = 0; i<nS;i++)
  230.         for (j = 0; j< nC; j++)
  231.             a[i][j] = 0;
  232. }
  233.  
  234. void MultiGridIteration(int deep, vector u, vector f, int nS, int nC)
  235. {
  236.     if (deep == MaxDeep)
  237.     {
  238.         RoughSolution(u,f,epsS,nS,nC);
  239.         return;
  240.     }  
  241.     int i;
  242.     for (i = 0; i < NSmoothIters; i++)
  243.         ZeidelIteration1(u,f,nS,nC);
  244.     //fout << "PreZejdel at " << deep << ":" << endl;
  245.     //PrintVector(u,nS,nC);
  246.     Remainder(f,u,tmp,nS,nC);
  247.     Restriction(tmp,r[deep + 1],nS,nC);
  248.     //Restriction(f,r[deep + 1],nS,nC);
  249.     //fout << "Restriction at " << deep << ":" << endl;
  250.     //PrintVector(r[deep + 1],nS / 2,nC / 2);
  251.     Zero(v[deep + 1],nS /2 ,nC / 2);
  252.     MultiGridIteration(deep + 1, v[deep + 1], r[deep + 1], nS / 2, nC / 2);
  253.     //fout << "Solution at " << deep <<":" << endl;
  254.     //PrintVector(v[deep + 1],nS / 2,nC / 2);
  255.     Prolongation(v[deep + 1], tmp, nS / 2, nC / 2);
  256.     //fout << "Prolongation at " << deep <<":" << endl;
  257.     //PrintVector(tmp,nS,nC);
  258.     Sum(u,tmp,nS,nC);
  259.     //Set(u,tmp,nS,nC);
  260.     for (i = 0; i < NSmoothIters; i++)
  261.         ZeidelIteration2(u,f,nS,nC);
  262.     //fout << "PostZejdel at " << deep << ":" << endl;
  263.     //PrintVector(u,nS,nC);
  264.     /*
  265.     Remainder(f,u,tmp,nS,nC);
  266.     Restriction(tmp,r[deep + 1],nS,nC);
  267.     MultiGridIteration(deep + 1, v[deep + 1], r[deep + 1], nS / 2, nC / 2);
  268.     Prolongation(v[deep + 1], tmp,nS / 2, nC / 2);
  269.     Sum(u,tmp,nS,nC);*/
  270. }
  271.  
  272. double Norm(vector u, int nS, int nC)
  273. {
  274.     int i, j;
  275.     double max = 0;
  276.     for (i = 0; i < nS; i++)
  277.         for (j = 0; j < nC; j++)
  278.             if (_abs(u[i][j]) > max) max = _abs(u[i][j]);
  279.     return max;
  280. }
  281.  
  282. void FormF()
  283. {
  284.     int i,j;
  285.     for (i = 0; i < NSize; i++)
  286.         for (j = 0; j < NCell; j++)
  287.         /*{
  288.             if ((i == 0) || (i == NSize - 1))
  289.             {
  290.                 if ((j == 0) || (j == NCell - 1)) f[i][j] = 2;
  291.                 else f[i][j] = 1;
  292.             }
  293.             else
  294.             {
  295.                 if ((j == 0) || (j == NCell - 1)) f[i][j] = 1;
  296.                 else f[i][j] = 0;
  297.             }
  298.         }*/
  299.         f[i][j] = 16;
  300. }
  301.  
  302. void Multigrid()
  303. {
  304.     Allocate();
  305.     vector solo  = new double* [NSize];
  306.     for (int i = 0; i < NSize; i++)
  307.         solo[i] = new double[NCell];
  308.     for (int i = 0; i< NSize;i++)
  309.         for (int j=0; j<NCell; j++)
  310.             solo[i][j] = 0;
  311.     FormF();
  312.     int iter = 0;
  313.     double NormOld, NormNew;
  314.     NormNew = Norm(f,NSize,NCell);
  315.     while (true)
  316.     {
  317.         Zero(v[0],NSize,NCell);
  318.         MultiGridIteration(0,v[0],f,NSize,NCell);
  319.         Sum(solo, v[0], NSize, NCell);
  320.         GetAuDecrementToRes(v[0],f,NSize,NCell);
  321.         iter++;
  322.         NormOld = NormNew;
  323.         NormNew = Norm(f,NSize,NCell);
  324.         cout << "Iteration " << iter << ": N1 / N2 = " << NormNew/NormOld << endl;
  325.         cout << "Remainer: " << NormNew << endl;
  326.         //fout << "Current solution:" << endl;
  327.         //PrintVector(solo,NSize,NCell);
  328.         //fout << "Remainder:" << endl;
  329.         /*PrintVector(f,NSize,NCell);
  330.         fout << "---------------------------------------" << endl;
  331.         fout << endl;
  332.         fout << endl;
  333.         fout << endl;
  334.         fout << endl;
  335.         fout << endl;
  336.         fout << endl;*/
  337.         //if (NormNew < epsM)
  338.         //  break;
  339.         if (iter > 30)
  340.             break;
  341.     }
  342.     fout << "Remainder : " << Norm(tmp,NSize,NCell);
  343.     //fout << endl << "Iterations of MultiGrid Method: " << iter << endl;
  344.     PrintVector(solo,NSize,NCell);
  345. }
  346.  
  347. int main()
  348. {
  349.     MaxDeep = (int) (log(NSize + 1.0) / log(2.0)) - 1;
  350.     cout << "Deep = " << MaxDeep << endl;
  351.     Multigrid();
  352.     cout << "Done!" << endl;
  353.     _getch();
  354. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement