C++ Bc. 33 cpp

#include <iostream>
#include <matvec/symmat.h>

void cholesky(GNU_gama::Mat<>& A);
void cholesky(GNU_gama::SymMat<>& A);

int main()
{
  { /* reseni s plnou matici */
    GNU_gama::Mat<> A(5,5);
    
    A =
      225, 210, 180, 135,  75, 
      210, 365, 311, 230, 122, 
      180, 311, 365, 266, 134, 
      135, 230, 266, 230, 110, 
       75, 122, 134, 110,  55; 

    std::cout.width(3);
    std::cout << A;
    cholesky(A);
    std::cout.width(3);
    std::cout << A;
  }

  { /* reseni se symetrickou matici */
    GNU_gama::SymMat<> A(5);
    
    A = 
      225,
      210, 365,
      180, 311, 365,
      135, 230, 266, 230,
       75, 122, 134, 110,  55; 
    
    std::cout.width(3);
    std::cout << A;
    cholesky(A);
    std::cout.width(3);
    std::cout << A;
  }
}


/* reseni s plnou matici */
void cholesky(GNU_gama::Mat<>& A)
{
  const int N = A.rows();

  for (int i=1; i<=N; i++)
    {
      // odmocnina z diagonalniho prvku
      const double d = std::sqrt(A(i,i));
      A(i,i) = d;

      // vynulovani diagonalniho radku
      for (int j=i+1; j<=N; j++) A(i,j) = 0; 

      // vypocet diagonalniho sloupce 
      for (int j=i+1; j<=N; j++) A(j,i) /= d;

      // vypocet submatice B
      for (int r=i+1; r<=N; r++)
        for (int s=i+1; s<=N; s++)
          A(r,s) -= A(r,i)*A(s,i);
    }
}


/* reseni se symetrickou matici */
void cholesky(GNU_gama::SymMat<>& A)
{
  const int N = A.rows();

  for (int i=1; i<=N; i++)
    {
      // odmocnina z diagonalniho prvku
      const double d = std::sqrt(A(i,i));
      A(i,i) = d;

      // vypocet diagonalniho sloupce 
      for (int j=i+1; j<=N; j++) A(j,i) /= d;

      // vypocet submatice B
      for (int r=i+1; r<=N; r++)
        for (int s=i+1; s<=r; s++)
          A(r,s) -= A(r,i)*A(s,i);
    }
}

  5   5

225 210 180 135  75 
210 365 311 230 122 
180 311 365 266 134 
135 230 266 230 110 
 75 122 134 110  55 
  5   5

 15   0   0   0   0 
 14  13   0   0   0 
 12  11  10   0   0 
  9   8   7   6   0 
  5   4   3   2   1 
  5

225 
210 365 
180 311 365 
135 230 266 230 
 75 122 134 110  55 
  5

 15 
 14  13 
 12  11  10 
  9   8   7   6 
  5   4   3   2   1 

[ Zpět ]