Skip to content

Instantly share code, notes, and snippets.

@lynsei
Forked from Yegorsh/cholesky_decomposition.cpp
Created December 21, 2021 03:27
Show Gist options
  • Save lynsei/c36a7b9873d613523c7c6823be68fce1 to your computer and use it in GitHub Desktop.
Save lynsei/c36a7b9873d613523c7c6823be68fce1 to your computer and use it in GitHub Desktop.
Cholesky decomposition/Разложение Холецкого на C++
#include <iostream>
#include <cmath>
using namespace std;
void FillMatrix(double **matrix, const int size){
for (int i=0; i<size; i++){
for (int j=0; j<size; j++){
cout << "a[" << i+1 << "][" << j+1 << "] = ";
cin >> matrix[i][j];
}
}
}
void PrintMatrix(double **matrix, const int size){
for (int i=0; i<size; i++){
for (int j=0; j<size; j++){
cout << matrix[i][j] << "\t";
}
cout << endl;
}
}
void clearMemory(double** matrix, int size){
for (int i = 0; i < size; i++) {
delete[] matrix[i];
}
delete [] matrix;
}
void Transpose(double **matrix, double **transposed, const int size){
for (int i=0; i<size; i++){
for (int j=0; j<size; j++){
transposed[i][j] = matrix[j][i];
}
}
}
int main(){
int size = 3;
double **A = new double *[size];
for (int i = 0; i < size; i++) {
A[i] = new double [size];
}
double **Transposed = new double *[size];
for (int i = 0; i < size; i++) {
Transposed[i] = new double [size];
}
double **L = new double *[size];
for (int i = 0; i < size; i++) {
L[i] = new double [size];
}
//FillMatrix(A, size);
A[0][0] =4; A[0][1] =12; A[0][2] =-16;
A[1][0] =12; A[1][1] =37; A[1][2] =-43;
A[2][0] =-16; A[2][1] =-43; A[2][2] =98;
Transpose(A, Transposed, size);
cout << "A =" << endl;
PrintMatrix(A, size);
for (int i=0; i<size; i++){
for (int j=0; j<size; j++){
if (Transposed[i][j] != A[i][j]){
cout << "Матрица не является симметричной." << endl;
exit(0);
}
}
}
for (int i = 0; i < size; i++) {
for (int j = 0; j <= i; j++) {
double sum = 0;
for (int k = 0; k < j; k++) {
sum += L[i][k] * L[j][k];
}
if (i == j) {
L[i][j] = sqrt(A[i][i] - sum);
if ((L[i][i] <= 0) || isnan(L[i][i])) {
cout << "Матрица не положительно определенная." << endl;
exit(0);
}
}
else {
L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
}
}
}
cout << endl << "L =" << endl;
PrintMatrix(L, size);
clearMemory(A,size);
clearMemory(Transposed, size);
clearMemory(L, size);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment