b2ac/ext/src/Eigenmethods_double.c
void inverseIteration_double(double *A, double* eigenvector,
		double eigenvalue, int32_t numberOfIterations)
{
	// Storage room for the inverse to form.
	int i;
	double determinant;
	double inverse[9];
	double tmp_storage[3];
	double norm_value = sqrt((double) 3);

	// Reset the eigenvector to a uniform unit vector.
	eigenvector[0] = 1 / norm_value;
	eigenvector[1] = eigenvector[0];
	eigenvector[2] = eigenvector[0];

	// Subtract the eigenvalue from the matrix.
	A[0] -= eigenvalue;
	A[4] -= eigenvalue;
	A[8] -= eigenvalue;

    // By not dividing by the determinant, but letting the normalization
    // in the end take care of that "scaling parameter", we are
    // more numerically stable against singular matrices.
	determinant = matrixInverse3by3_double(A, inverse);

	for (i = 0; i < numberOfIterations; i++)
	{
		// Perform matrix-vector multiplication.
		tmp_storage[0] = inverse[0] * eigenvector[0] +
						 inverse[1] * eigenvector[1] +
						 inverse[2] * eigenvector[2];
		tmp_storage[1] = inverse[3] * eigenvector[0] +
						 inverse[4] * eigenvector[1] +
						 inverse[5] * eigenvector[2];
		tmp_storage[2] = inverse[6] * eigenvector[0] +
						 inverse[7] * eigenvector[1] +
						 inverse[8] * eigenvector[2];

		// Normalise the eigenvector candidate.
		norm_value = sqrt(tmp_storage[0] * tmp_storage[0] +
						  tmp_storage[1] * tmp_storage[1] +
						  tmp_storage[2] * tmp_storage[2]);
		eigenvector[0] = tmp_storage[0] / norm_value;
		eigenvector[1] = tmp_storage[1] / norm_value;
		eigenvector[2] = tmp_storage[2] / norm_value;
	}

	// Add the eigenvalue back to the matrix.
	A[0] += eigenvalue;
	A[4] += eigenvalue;
	A[8] += eigenvalue;
}
/**
 * Inverse iteration for finding eigenvector.
 *
 * @brief Given a matrix and a known eigenvalue, this function
 * multiplies the matrix (A - mu*I)^{-1} with a iteratively improved
 * eigenvector approximation.
 *
 * Note that for a accurate eigenvalue, one one iteration is required for
 * a reasonable approximation.
 *
 * References:
 * [1] http://en.wikipedia.org/wiki/Inverse_iteration
 * [2] Gene H. Golub and Charles F. Van Loan. Matrix computations (3rd
 *     ed.). Johns Hopkins University Press, Baltimore, MD, USA, 1996.
 *
 */
Inverse iteration for finding eigenvector.
@brief Given a matrix and a known eigenvalue, this function
multiplies the matrix (A - mu*I)^{-1} with a iteratively improved
eigenvector approximation.

Note that for a accurate eigenvalue, one one iteration is required for
a reasonable approximation.


Inverse iteration for finding eigenvector . @brief Given a matrix and a known eigenvalue this function multiplies the matrix ( A - mu * I ) ^ { - 1 } with a iteratively improved eigenvector approximation . Note that for a accurate eigenvalue one one iteration is required for a reasonable approximation .
void inverseIteration_double(double *A, double* eigenvector,
		double eigenvalue, int32_t numberOfIterations)
{
	int i;
	double determinant;
	double inverse[9];
	double tmp_storage[3];
	double norm_value = sqrt((double) 3);
	eigenvector[0] = 1 / norm_value;
	eigenvector[1] = eigenvector[0];
	eigenvector[2] = eigenvector[0];
	A[0] -= eigenvalue;
	A[4] -= eigenvalue;
	A[8] -= eigenvalue;
	determinant = matrixInverse3by3_double(A, inverse);
	for (i = 0; i < numberOfIterations; i++)
	{
		tmp_storage[0] = inverse[0] * eigenvector[0] +
						 inverse[1] * eigenvector[1] +
						 inverse[2] * eigenvector[2];
		tmp_storage[1] = inverse[3] * eigenvector[0] +
						 inverse[4] * eigenvector[1] +
						 inverse[5] * eigenvector[2];
		tmp_storage[2] = inverse[6] * eigenvector[0] +
						 inverse[7] * eigenvector[1] +
						 inverse[8] * eigenvector[2];
		norm_value = sqrt(tmp_storage[0] * tmp_storage[0] +
						  tmp_storage[1] * tmp_storage[1] +
						  tmp_storage[2] * tmp_storage[2]);
		eigenvector[0] = tmp_storage[0] / norm_value;
		eigenvector[1] = tmp_storage[1] / norm_value;
		eigenvector[2] = tmp_storage[2] / norm_value;
	}
	A[0] += eigenvalue;
	A[4] += eigenvalue;
	A[8] += eigenvalue;
}
void inverseIteration_double ( double * A , double * eigenvector , double eigenvalue , int32_t numberOfIterations ) { int i ; double determinant ; double inverse [ 9 ] ; double tmp_storage [ 3 ] ; double norm_value = sqrt ( ( double ) 3 ) ; eigenvector [ 0 ] = 1 / norm_value ; eigenvector [ 1 ] = eigenvector [ 0 ] ; eigenvector [ 2 ] = eigenvector [ 0 ] ; A [ 0 ] -= eigenvalue ; A [ 4 ] -= eigenvalue ; A [ 8 ] -= eigenvalue ; determinant = matrixInverse3by3_double ( A , inverse ) ; for ( i = 0 ; i < numberOfIterations ; i ++ ) { tmp_storage [ 0 ] = inverse [ 0 ] * eigenvector [ 0 ] + inverse [ 1 ] * eigenvector [ 1 ] + inverse [ 2 ] * eigenvector [ 2 ] ; tmp_storage [ 1 ] = inverse [ 3 ] * eigenvector [ 0 ] + inverse [ 4 ] * eigenvector [ 1 ] + inverse [ 5 ] * eigenvector [ 2 ] ; tmp_storage [ 2 ] = inverse [ 6 ] * eigenvector [ 0 ] + inverse [ 7 ] * eigenvector [ 1 ] + inverse [ 8 ] * eigenvector [ 2 ] ; norm_value = sqrt ( tmp_storage [ 0 ] * tmp_storage [ 0 ] + tmp_storage [ 1 ] * tmp_storage [ 1 ] + tmp_storage [ 2 ] * tmp_storage [ 2 ] ) ; eigenvector [ 0 ] = tmp_storage [ 0 ] / norm_value ; eigenvector [ 1 ] = tmp_storage [ 1 ] / norm_value ; eigenvector [ 2 ] = tmp_storage [ 2 ] / norm_value ; } A [ 0 ] += eigenvalue ; A [ 4 ] += eigenvalue ; A [ 8 ] += eigenvalue ; }
Inverse iteration for finding eigenvector.
Inverse iteration for finding eigenvector .
// Storage room for the inverse to form. // Reset the eigenvector to a uniform unit vector. // Subtract the eigenvalue from the matrix. // By not dividing by the determinant, but letting the normalization // in the end take care of that "scaling parameter", we are // more numerically stable against singular matrices. // Perform matrix-vector multiplication. // Normalise the eigenvector candidate. // Add the eigenvalue back to the matrix.