b2ac/ext/src/Eigenmethods_double.c
void QRAlgorithm_double(double *A, double *eigenvalues)
{
	double mu;
	int k;
	int m = 2;
	int mMapping[3] = {0, 4, 8};

	double B[9];
	// Copy matrix, since we destroy it in the process of
	// finding eigenvalues.
	for (k = 0; k < 9; k++)
	{
		B[k] = A[k];
	}

	// First, convert the A to Hessenberg form.
	matrixToHessenbergForm_Givens_double(B);

	while (m > 0)
	{
		// Apply the shift to improve algorithm convergence.
		// The index (m * 2) refer to the diagonal entries of the
		// B. Set the shift value to the "lower right corner" value.
		mu = B[mMapping[m]];
		for (k = 0; k <= m; k++)
		{
			B[mMapping[k]] -= mu;
		}

		// Run the Givens rotation producing RQ from B.
		performQRStep_Givens_double(B, m);

		// Add the shift to the diagonal elements again.
		for (k = 0; k <= m; k++)
		{
			B[mMapping[k]] += mu;
		}

		// If B[m, m-1] is sufficiently small, an eigenvalue can be
		// deemed to have been found.

		if (fabs(B[mMapping[m] - 1]) < QR_ALGORITHM_TOLERANCE_DOUBLE)
		{
			eigenvalues[m] = B[mMapping[m]];
			m--;
		}
	}

	// Two eigenvalues have been found in the loop above. The final one is
	// left in the A A top left corner.
	eigenvalues[0] = B[0];

	// Also sort the eigenvalues with simple Bubblesort.
	// Use mu as temporary variable for sorting.
	if (eigenvalues[0] > eigenvalues[1])
	{
		mu = eigenvalues[1];
		eigenvalues[1] = eigenvalues[0];
		eigenvalues[0] = mu;
	}
	if (eigenvalues[1] > eigenvalues[2])
	{
		mu = eigenvalues[2];
		eigenvalues[2] = eigenvalues[1];
		eigenvalues[1] = mu;
	}
	if (eigenvalues[0] > eigenvalues[1])
	{
		mu = eigenvalues[1];
		eigenvalues[1] = eigenvalues[0];
		eigenvalues[0] = mu;
	}
}
/**
 * The QR Algorithm, for finding eigenvalues.
 *
 * @brief This version of the QR algorithms is implemented
 * purely using Givens rotations. It fills the input
 * vector eigenvalues with the eigenvalues it detects.
 */
The QR Algorithm, for finding eigenvalues.
@brief This version of the QR algorithms is implemented
purely using Givens rotations. It fills the input
vector eigenvalues with the eigenvalues it detects.
The QR Algorithm for finding eigenvalues . @brief This version of the QR algorithms is implemented purely using Givens rotations . It fills the input vector eigenvalues with the eigenvalues it detects .
void QRAlgorithm_double(double *A, double *eigenvalues)
{
	double mu;
	int k;
	int m = 2;
	int mMapping[3] = {0, 4, 8};
	double B[9];
	for (k = 0; k < 9; k++)
	{
		B[k] = A[k];
	}
	matrixToHessenbergForm_Givens_double(B);
	while (m > 0)
	{
		mu = B[mMapping[m]];
		for (k = 0; k <= m; k++)
		{
			B[mMapping[k]] -= mu;
		}
		performQRStep_Givens_double(B, m);
		for (k = 0; k <= m; k++)
		{
			B[mMapping[k]] += mu;
		}
		if (fabs(B[mMapping[m] - 1]) < QR_ALGORITHM_TOLERANCE_DOUBLE)
		{
			eigenvalues[m] = B[mMapping[m]];
			m--;
		}
	}
	eigenvalues[0] = B[0];
	if (eigenvalues[0] > eigenvalues[1])
	{
		mu = eigenvalues[1];
		eigenvalues[1] = eigenvalues[0];
		eigenvalues[0] = mu;
	}
	if (eigenvalues[1] > eigenvalues[2])
	{
		mu = eigenvalues[2];
		eigenvalues[2] = eigenvalues[1];
		eigenvalues[1] = mu;
	}
	if (eigenvalues[0] > eigenvalues[1])
	{
		mu = eigenvalues[1];
		eigenvalues[1] = eigenvalues[0];
		eigenvalues[0] = mu;
	}
}
void QRAlgorithm_double ( double * A , double * eigenvalues ) { double mu ; int k ; int m = 2 ; int mMapping [ 3 ] = { 0 , 4 , 8 } ; double B [ 9 ] ; for ( k = 0 ; k < 9 ; k ++ ) { B [ k ] = A [ k ] ; } matrixToHessenbergForm_Givens_double ( B ) ; while ( m > 0 ) { mu = B [ mMapping [ m ] ] ; for ( k = 0 ; k <= m ; k ++ ) { B [ mMapping [ k ] ] -= mu ; } performQRStep_Givens_double ( B , m ) ; for ( k = 0 ; k <= m ; k ++ ) { B [ mMapping [ k ] ] += mu ; } if ( fabs ( B [ mMapping [ m ] - 1 ] ) < QR_ALGORITHM_TOLERANCE_DOUBLE ) { eigenvalues [ m ] = B [ mMapping [ m ] ] ; m -- ; } } eigenvalues [ 0 ] = B [ 0 ] ; if ( eigenvalues [ 0 ] > eigenvalues [ 1 ] ) { mu = eigenvalues [ 1 ] ; eigenvalues [ 1 ] = eigenvalues [ 0 ] ; eigenvalues [ 0 ] = mu ; } if ( eigenvalues [ 1 ] > eigenvalues [ 2 ] ) { mu = eigenvalues [ 2 ] ; eigenvalues [ 2 ] = eigenvalues [ 1 ] ; eigenvalues [ 1 ] = mu ; } if ( eigenvalues [ 0 ] > eigenvalues [ 1 ] ) { mu = eigenvalues [ 1 ] ; eigenvalues [ 1 ] = eigenvalues [ 0 ] ; eigenvalues [ 0 ] = mu ; } }
The QR Algorithm, for finding eigenvalues.
The QR Algorithm for finding eigenvalues .
// Copy matrix, since we destroy it in the process of // finding eigenvalues. // First, convert the A to Hessenberg form. // Apply the shift to improve algorithm convergence. // The index (m * 2) refer to the diagonal entries of the // B. Set the shift value to the "lower right corner" value. // Run the Givens rotation producing RQ from B. // Add the shift to the diagonal elements again. // If B[m, m-1] is sufficiently small, an eigenvalue can be // deemed to have been found. // Two eigenvalues have been found in the loop above. The final one is // left in the A A top left corner. // Also sort the eigenvalues with simple Bubblesort. // Use mu as temporary variable for sorting.