Метод вращений решения линейных систем

        s[k][l] = a[l][k] / (sqrt( a[k][k]*a[k][k] + a[l][k]*a[l][k] ));

 

           // Умножение матрицы a[][] на T[k][l]

           akk=a[k][k];  alk=a[l][k];  akl=a[k][l];  all=a[l][l];

           a[k][k] =  akk*c[k][l] + alk*s[k][l];

           a[k][l] =  akl*c[k][l] + all*s[k][l];

           a[l][k] = -akk*s[k][l] + alk*c[k][l];

           a[l][l] = -akl*s[k][l] + all*c[k][l];

 

           // Вектор свободных членов умножается на T[k][l]

           bk = b[k]; bl = b[l];

           b[k] =  bk*c[k][l] + bl*s[k][l];

           b[l] = -bk*s[k][l] + bl*c[k][l];

       }

      }

      // Теперь матрица a[][] -- верхняя  диагональная.

 

      // Обратный ход

      float h;

      x[n-1]=b[n-1]/a[n-1][n-1];

      for (l=(n-1);l>=1;l--) {

          h=b[l-1];

          for (k=(l+1);k<=n;k++) h=h-x[k-1]*a[l-1][k-1];

          x[l-1]=h/a[l-1][l-1];

      }

 

    }

 
 

    // ------------------------------------------------------

    // Main

    // ------------------------------------------------------

    void main() {

     clrscr();

     matrix A(_N,_M);

     for (int i=0; i<_N; i++) for (int j=0; j<_M; j++) A.seta(i,j,i+j);

     for (i=0; i<_N; i++) A.setb(i,1);

      printf("Решение СЛАУ методом вращения (QR-разложение)\n");

     printf("=Исходная система========================================\n");

      A.print_slau();

      printf("\n=Решение системы=========================================\n");

      A.qr();

     A.print_x();

      getch();

    }