Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <fstream>
- #include <cmath>
- #include <conio.h>
- using namespace std;
- typedef double** vector;
- const double Star[3][3] = { 0.25, 0.5, 0.25,
- 0.5, 1, 0.5,
- 0.25, 0.5, 0.25 };
- int NSize = 1023; // размер блочной матрицы
- int NCell = 1023; // размер блока
- int NSmoothIters = 1; // количество сглаживающих итераций
- int MaxDeep = 2; // глубина рекурсии
- double epsS = 1e-8;
- double epsM = 1e-8;
- vector f, tmp;
- vector *v, *r;
- fstream fout("out.txt",ios::out);
- double _abs(double a)
- {
- return (a > 0 ? a : -a);
- }
- void Set(vector a, vector b, int nS, int nC)
- {
- int i,j;
- for (i=0;i<nS;i++)
- for (j=0;j<nC;j++)
- a[i][j]=b[i][j];
- }
- void Allocate()
- {
- int i,j,k;
- int nS = NSize, nC = NCell;
- v = new vector[MaxDeep + 1];
- r = new vector[MaxDeep + 1];
- f = new double* [NSize];
- tmp = new double* [NSize];
- for (i = 0; i <= MaxDeep; i++)
- {
- v[i] = new double* [nS];
- r[i] = new double* [nS];
- for (j = 0; j < nS; j++)
- {
- v[i][j] = new double [nC];
- for (k = 0; k < nC; k++)
- v[i][j][k] = 0;
- r[i][j] = new double [nC];
- for (k = 0; k < nC; k++)
- r[i][j][k] = 0;
- }
- nS = nS / 2;
- nC = nC / 2;
- }
- for (i = 0; i < NSize; i++)
- {
- f[i] = new double [NCell];
- tmp[i] = new double [NCell];
- }
- }
- void Restriction(vector a, vector b, int nS, int nC) // блочный вектор a[nS][nC] переводим
- // в блочный вектор b[nS/2][nC/2];
- {
- int i,j,k,l,ii,jj, nS2 = nS/2, nC2 = nC/2;
- for (i = 0; i < nS2; i++)
- for (j = 0; j < nC2; j++)
- {
- ii = 2 * i + 1;
- jj = 2 * j + 1;
- b[i][j] = 0;
- for (k = -1; k < 2; k++)
- for (l = -1; l < 2; l++)
- if ((ii + k < nS) && (ii + k >= 0) && (jj + l < nC) && (jj + l >= 0))
- b[i][j] += a[ii + k][jj + l] * Star [k + 1][l + 1];
- }
- }
- void Prolongation(vector a, vector b, int nS, int nC) //блочный вектор a[nS][nC] переводим
- // в блочный вектор b[2*nS + 1][2*nC + 1];
- {
- int i,j,k,l,ii,jj, nS2 = 2*nS + 1, nC2 = 2*nC + 1;
- for (i = 0; i < nS2; i++)
- for (j = 0; j < nC2; j++)
- b[i][j] = 0;
- for (i = 0; i < nS; i++)
- for (j = 0; j < nC; j++)
- {
- ii = 2 * i + 1;
- jj = 2 * j + 1;
- b[ii][jj] = 0;
- for (k = -1; k < 2; k++)
- for (l = -1; l < 2; l++)
- if ((ii + k < nS2) && (ii + k >= 0) && (jj + l < nC2) && (jj + l >= 0))
- b[ii + k][jj + l] += a[i][j] * Star [k + 1][l + 1];
- }
- }
- void ZeidelIteration1(vector u, vector f, int nS, int nC)
- {
- int i, j;
- for (i = 0; i < nS; i++)
- for (j = 0; j < nC; j++)
- {
- u[i][j] = f[i][j];
- if (i != 0) u[i][j] += u[i-1][j];
- if (i != nS - 1) u[i][j] += u[i+1][j];
- if (j != 0) u[i][j] += u[i][j-1];
- if (j != nC - 1) u[i][j] += u[i][j+1];
- u[i][j] = u[i][j] / 4;
- }
- }
- void ZeidelIteration2(vector u, vector f, int nS, int nC)
- {
- int i, j;
- for (i = nS-1; i >= 0; i--)
- for (j = nC - 1; j >= 0; j--)
- {
- u[i][j] = f[i][j];
- if (i != 0) u[i][j] += u[i-1][j];
- if (i != nS - 1) u[i][j] += u[i+1][j];
- if (j != 0) u[i][j] += u[i][j-1];
- if (j != nC - 1) u[i][j] += u[i][j+1];
- u[i][j] = u[i][j] / 4;
- }
- }
- void ZeidelMethod(vector u, vector f, double eps, int nS, int nC)
- {
- int i,j;
- double max, tmp;
- while (true)
- {
- max = 0;
- for (i = 0; i < nS; i++)
- for (j = 0; j < nC; j++)
- {
- tmp = f[i][j];
- if (i != 0) tmp += u[i-1][j];
- if (i != nS - 1) tmp += u[i+1][j];
- if (j != 0) tmp += u[i][j-1];
- if (j != nC - 1) tmp += u[i][j+1];
- tmp = tmp / 4;
- if (_abs(tmp - u[i][j]) > max) max = _abs(tmp - u[i][j]);
- u[i][j] = tmp;
- }
- if (max < eps)
- return;
- }
- }
- void SmoothingIteration(vector u, vector f, int nS, int nC) // сглаживающая итерация
- {
- // ZeidelIteration(u,f,nS,nC);
- }
- void RoughSolution(vector u, vector f, double eps, int nS, int nC) // решение на грубой сетке
- {
- ZeidelMethod(u,f,eps,nS,nC);
- }
- void Remainder(vector b, vector u, vector res, int nS, int nC) // res = b - Au
- {
- int i, j;
- for (i = 0; i < nS; i++)
- for (j = 0; j < nC; j++)
- {
- res[i][j] = b[i][j] - 4 * u[i][j];
- if (i != 0) res[i][j] += u[i-1][j];
- if (i != nS - 1) res[i][j] += u[i+1][j];
- if (j != 0) res[i][j] += u[i][j-1];
- if (j != nC - 1) res[i][j] += u[i][j+1];
- }
- }
- void GetAuDecrementToRes(vector u, vector res, int nS, int nC) // res = b - Au
- {
- int i, j;
- for (i = 0; i < nS; i++)
- for (j = 0; j < nC; j++)
- {
- res[i][j] -= 4 * u[i][j];
- if (i != 0) res[i][j] += u[i-1][j];
- if (i != nS - 1) res[i][j] += u[i+1][j];
- if (j != 0) res[i][j] += u[i][j-1];
- if (j != nC - 1) res[i][j] += u[i][j+1];
- }
- }
- void Sum(vector a, vector b, int nS, int nC) // результат в a
- {
- for (int i = 0; i < nS; i++)
- for (int j = 0; j < nC; j++)
- a[i][j] += b[i][j];
- }
- void Diff(vector a, vector b, int nS, int nC) // результат в a
- {
- for (int i = 0; i < nS; i++)
- for (int j = 0; j < nC; j++)
- a[i][j] -= b[i][j];
- }
- void PrintVector(vector a, int nS, int nC)
- {
- int i,j;
- for (i = 0; i < nS; i++)
- {
- for (j = 0; j < nC; j++)
- {
- fout.width(10);
- fout.precision(4);
- fout << a[i][j];
- }
- fout << endl;
- }
- }
- void Zero(vector a, int nS, int nC)
- {
- int i,j;
- for (i = 0; i<nS;i++)
- for (j = 0; j< nC; j++)
- a[i][j] = 0;
- }
- void MultiGridIteration(int deep, vector u, vector f, int nS, int nC)
- {
- if (deep == MaxDeep)
- {
- RoughSolution(u,f,epsS,nS,nC);
- return;
- }
- int i;
- for (i = 0; i < NSmoothIters; i++)
- ZeidelIteration1(u,f,nS,nC);
- //fout << "PreZejdel at " << deep << ":" << endl;
- //PrintVector(u,nS,nC);
- Remainder(f,u,tmp,nS,nC);
- Restriction(tmp,r[deep + 1],nS,nC);
- //Restriction(f,r[deep + 1],nS,nC);
- //fout << "Restriction at " << deep << ":" << endl;
- //PrintVector(r[deep + 1],nS / 2,nC / 2);
- Zero(v[deep + 1],nS /2 ,nC / 2);
- MultiGridIteration(deep + 1, v[deep + 1], r[deep + 1], nS / 2, nC / 2);
- //fout << "Solution at " << deep <<":" << endl;
- //PrintVector(v[deep + 1],nS / 2,nC / 2);
- Prolongation(v[deep + 1], tmp, nS / 2, nC / 2);
- //fout << "Prolongation at " << deep <<":" << endl;
- //PrintVector(tmp,nS,nC);
- Sum(u,tmp,nS,nC);
- //Set(u,tmp,nS,nC);
- for (i = 0; i < NSmoothIters; i++)
- ZeidelIteration2(u,f,nS,nC);
- //fout << "PostZejdel at " << deep << ":" << endl;
- //PrintVector(u,nS,nC);
- /*
- Remainder(f,u,tmp,nS,nC);
- Restriction(tmp,r[deep + 1],nS,nC);
- MultiGridIteration(deep + 1, v[deep + 1], r[deep + 1], nS / 2, nC / 2);
- Prolongation(v[deep + 1], tmp,nS / 2, nC / 2);
- Sum(u,tmp,nS,nC);*/
- }
- double Norm(vector u, int nS, int nC)
- {
- int i, j;
- double max = 0;
- for (i = 0; i < nS; i++)
- for (j = 0; j < nC; j++)
- if (_abs(u[i][j]) > max) max = _abs(u[i][j]);
- return max;
- }
- void FormF()
- {
- int i,j;
- for (i = 0; i < NSize; i++)
- for (j = 0; j < NCell; j++)
- /*{
- if ((i == 0) || (i == NSize - 1))
- {
- if ((j == 0) || (j == NCell - 1)) f[i][j] = 2;
- else f[i][j] = 1;
- }
- else
- {
- if ((j == 0) || (j == NCell - 1)) f[i][j] = 1;
- else f[i][j] = 0;
- }
- }*/
- f[i][j] = 16;
- }
- void Multigrid()
- {
- Allocate();
- vector solo = new double* [NSize];
- for (int i = 0; i < NSize; i++)
- solo[i] = new double[NCell];
- for (int i = 0; i< NSize;i++)
- for (int j=0; j<NCell; j++)
- solo[i][j] = 0;
- FormF();
- int iter = 0;
- double NormOld, NormNew;
- NormNew = Norm(f,NSize,NCell);
- while (true)
- {
- Zero(v[0],NSize,NCell);
- MultiGridIteration(0,v[0],f,NSize,NCell);
- Sum(solo, v[0], NSize, NCell);
- GetAuDecrementToRes(v[0],f,NSize,NCell);
- iter++;
- NormOld = NormNew;
- NormNew = Norm(f,NSize,NCell);
- cout << "Iteration " << iter << ": N1 / N2 = " << NormNew/NormOld << endl;
- cout << "Remainer: " << NormNew << endl;
- //fout << "Current solution:" << endl;
- //PrintVector(solo,NSize,NCell);
- //fout << "Remainder:" << endl;
- /*PrintVector(f,NSize,NCell);
- fout << "---------------------------------------" << endl;
- fout << endl;
- fout << endl;
- fout << endl;
- fout << endl;
- fout << endl;
- fout << endl;*/
- //if (NormNew < epsM)
- // break;
- if (iter > 30)
- break;
- }
- fout << "Remainder : " << Norm(tmp,NSize,NCell);
- //fout << endl << "Iterations of MultiGrid Method: " << iter << endl;
- PrintVector(solo,NSize,NCell);
- }
- int main()
- {
- MaxDeep = (int) (log(NSize + 1.0) / log(2.0)) - 1;
- cout << "Deep = " << MaxDeep << endl;
- Multigrid();
- cout << "Done!" << endl;
- _getch();
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement