ulvis.paste.net

Paste Search Dynamic
Recent pastes
rowswap
  1. #include<iostream>
  2. #include<cmath>
  3.  
  4. using namespace std;
  5.  
  6. void RowSwap(double *A, double *b, int i, int j, int n)
  7. {
  8.         double temp;
  9.         for (int l = 0; l < n; l++)
  10.         {
  11.                 temp = A[j + n * l];
  12.                 A[j + n * l] = A[i + n * l];
  13.                 A[i + n * l] = temp;
  14.         }
  15.  
  16.         temp = b[j];
  17.         b[j] = b[i];
  18.         b[i] = temp;
  19.  
  20.         return;
  21. }
  22.  
  23. void Matvecmul(double *A, double *v, int n)
  24. {
  25.         double *t = new double[n];
  26.         for (int i = 0; i < n; i++)
  27.         {
  28.                 t[i] = 0;
  29.                 for (int j = 0; j < n; j++)
  30.                 {
  31.                         t[i] += A[i + n * j] * v[j];
  32.                 }
  33.         }
  34.  
  35.         for (int i = 0; i < n; i++)
  36.         {
  37.                 v[i] = t[i];
  38.         }
  39.         delete[] t;
  40.         return;
  41. }
  42.  
  43. void Matmatmul(double *A, double *B, int n)
  44. {
  45.         double *C = new double[n*n];
  46.         for (int i = 0; i < n; i++)
  47.         {
  48.                 for (int j = 0; j < n; j++)
  49.                 {
  50.                         C[i + n * j] = 0;
  51.                         for (int k = 0; k < n; k++)
  52.                         {
  53.                                 C[i + n * j] += A[i + n * k] * B[k + n * j];
  54.                         }
  55.                 }
  56.         }
  57.  
  58.         for (int i = 0; i < n; i++)
  59.         {
  60.                 for (int j = 0; j < n; j++)
  61.                 {
  62.                         B[i + n * j] = C[i + n * j];
  63.                 }
  64.         }
  65.        
  66.         delete[] C;
  67.         return;
  68. }
  69.  
  70. void Householder(double *A, double *v, int n)
  71. {
  72.         double sum = 0;
  73.         for (int i = 0; i < n; i++)
  74.         {
  75.                 sum += v[i] * v[i];
  76.         }
  77.  
  78.         for (int i = 0; i < n; i++)
  79.         {
  80.                 for (int j = 0; j < n; j++)
  81.                 {
  82.                         A[i + n * j] = - 2 * v[i] * v[j] / sum;
  83.                         if (i == j)
  84.                                 A[i + n * j] += 1;
  85.                 }
  86.         }
  87. }
  88.  
  89. int main()
  90. {
  91.         int n = 100;
  92.         double *A = new double[n*n];
  93.         double *b = new double[n];
  94.        
  95.         for (int i = 0; i < n; i++)
  96.         {
  97.                 for (int j = 0; j < n; j++)
  98.                 {
  99.                         A[i + n * j] = 0;
  100.                         if (abs(i - j) == 1)
  101.                                 A[i + n * j] = -1.0 * (n + 1) * (n + 1) ;
  102.                         if (i == j)
  103.                                 A[i + n * j] = 2.0 * (n + 1) * (n + 1);
  104.                 }
  105.                 b[i] = 0;
  106.                 if (i == 0)
  107.                         b[i] = 1;
  108.         }
  109.        
  110.         for (int i = 0; i < n - 1; i++)
  111.         {
  112.                 int maxindex = i;
  113.                 for (int k = i; k < n; k++)
  114.                 {
  115.                         if (abs(A[maxindex + n * i]) < abs(A[k + n * i]))
  116.                                 maxindex = k;
  117.                 }
  118.                 RowSwap(A, b, i, maxindex, n);
  119.                
  120.                 double *v = new double[n - i];
  121.                 double *H = new double[(n - i)*(n - i)];
  122.                 double sum = 0;
  123.                 for (int j = 0; j < n - i; j++)
  124.                 {
  125.                         sum += A[j + i + n * i] * A[j + i + n * i];
  126.                         v[j] = A[j + i + n * i];
  127.                 }
  128.                 v[0] -= sqrt(sum);
  129.        
  130.                 Householder(H, v, n - i);
  131.                 delete[] v;
  132.                
  133.                 double *w1 = new double[n - i];
  134.                 double *w2 = new double[n - i];
  135.                 double *K = new double[(n - i)*(n - i)];
  136.  
  137.                 for (int j = 0; j < n - i; j++)
  138.                 {
  139.                         w1[j] = A[j + i + n * i];
  140.                         w2[j] = b[j + i];
  141.                 }
  142.  
  143.                 for (int j = 0; j < n - i; j++)
  144.                 {
  145.                         for (int k = 0; k < n - i; k++)
  146.                         {
  147.                                 K[j + (n - i) * k] = A[j + i + n * (k + i)];
  148.                         }
  149.                 }
  150.                
  151.                 Matvecmul(H, w1, n - i);
  152.                 Matvecmul(H, w2, n - i);
  153.                 Matmatmul(H, K, n - i);
  154.                
  155.                 for (int j = 0; j < n - i; j++)
  156.                 {
  157.                         A[j + i + n * i] = w1[j];
  158.                         b[j + i] = w2[j];
  159.                 }
  160.  
  161.                 for (int j = 0; j < n - i; j++)
  162.                 {
  163.                         for (int k = 0; k < n - i; k++)
  164.                         {
  165.                                 A[j + i + n * (k + i)] = K[j + (n - i) * k];
  166.                         }
  167.                 }
  168.                        
  169.                 delete[] w1;
  170.                 delete[] w2;
  171.                 delete[] K;
  172.                 delete[] H;
  173.         }
  174.  
  175.         double *x = new double[n];
  176.  
  177.         for (int i = n - 1; i >= 0; i--)
  178.         {
  179.                 double sum = 0;
  180.                 for (int j = i + 1; j < n ; j++)
  181.                 {
  182.                         sum += A[i + n * j] * x[j];
  183.                 }
  184.                 x[i] = (b[i] - sum) / A[i + n * i];
  185.         }
  186.  
  187.         for (int i = 0; i < n; i++)
  188.         {
  189.                 cout << "x_" << i + 1 << " = " << x[i] << endl;
  190.         }
  191.  
  192.         return 0;
  193. }
  194.  
Parsed in 0.024 seconds