cholesky decomposition

[naufal]

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.

 
      #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;
}
    

[lastchance]

Lines 37 to 40 should be wrapped in an "else".

The last L matrix on line 39 needs conjugating.

[naufal]

Thank you. It is solved.

[naufal]

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]

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.