C++

Forum

 
cholesky decomposition

naufal (5)

Hi,
I would like to solve a complex hermitian positive definite matrix using c++. This is my code. Luckily, there is no error, but it gives me an accurate answer. Hope to hear from you.
Thank you. 


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
 
 
      #include <iostream>
#include <cmath>
#include<fstream>
#include<complex>
using namespace std;

const int order = 4;


void print(complex<double> M[order][order])
{
	for (int i = 0; i < order; i++)
	{
		for (int j = 0; j < order; j++) cout << M[i][j] << '\t';
		cout << '\n';
	}
}
//Function of cholesky
void cholDec(complex <double> sum[10][10], complex <double> L[4][4], complex <double> Z[4][4], int order)
{

	for (int i = 0; i < order; i++)
	{
		for (int j = 0; j <= i; j++)
		{
			sum[i][j] = 0;

			if (j == i)//summation of diagonals
			{
				for (int k = 0; k < j; k++)

					sum[i][j] += L[j][k] * conj(L[j][k]);
				L[j][j] = sqrt(Z[j][j] - sum[i][j]);
			}

			//Evaluating L(i,j) using L(j,j)
			for (int k = 0; k < j; k++)
				sum[i][j] += (L[i][k] * L[j][k]);
			L[i][j] = (Z[i][j] - sum[i][j]) / L[j][j];

		}

	}
}
	

//Main function
//All the operations is done here
int main(int argc, char** argv)
//int main()
{
	complex <double>  L[4][4], sum[10][10], Lcon[10][10], matrix[10][10], matrixconj[10][10];
	int order = 4;
	complex <double> Z[4][4] = { { 18.44 , -4.2 - 3.8i, 3.5 - 2.3i, 8.0 + 4.2i },
	{ -4.2 - 3.8i, 20.5, 0.3 + 4.0i, 7.8 + 2.2i },{ 3.4 - 10.1i, 8.1 + 3.2i, 17.3, 6.2 - 1.8i },{ 0.3 + 3.8i, 10.4 - 12.1i, 5.4 - 8.4i, 26.3 } };
	
	
	cout << "Original matrix, Z:\n";
	print(Z);
	cout << "\n";


	//call cholesky function
	cholDec(sum, L, Z, order);

	cout << "Lower triangular matrix is\n";
	print(L);
	cout << "\n";

	// Lower Triangular
	cout << "Lower triangular matrix is" << endl;
	for (int i = 0; i < order; i++)
	{
		for (int j = 0; j < order; j++)
		{
			std::cout << L[i][j] << ' ';
		}
		std::cout << std::endl;
		
	}
	cout << "\n";

    // Lower Triangular complex conjugate
    cout << "complex conjugate of Lower triangular  matrix is" << endl;
	for (int i = 0; i < order; i++)
	{
		for (int j = 0; j < order; j++)
		{
			Lcon[i][j] = conj(L[j][i]);
			std::cout << Lcon[i][j] << ' ';
		}
		std::cout << std::endl;
	}

	cout << "\n";
	system("pause");
	return 0;
}
Last edited on
lastchance (6980)
Lines 37 to 40 should be wrapped in an "else".

The last L matrix on line 39 needs conjugating.
naufal (5)
Thank you. It is solved.
naufal (5)
I want to have a new matrix called A. How to write a code in c++ if is created based on the modified version of the lower triangular matrix L?

The code in MATLAB is like this:

A=L;
A(1,4) = L(4,1);
A(3,4) = L(1,1);

May I know how to write in c++?
Thank you.
lastchance (6980)
You will find it considerably more easy to copy whole arrays if you use vectors ... as suggested in your other thread. Otherwise you will need a double loop just to copy an array.

EDIT:
BTW, I've just noticed that your original matrix isn't Hermitian.
Last edited on
Topic archived. No new replies allowed.