22 #include <lal/LALConstants.h>
23 #include <lal/LALStdlib.h>
24 #include <lal/AVFactories.h>
25 #include <lal/MatrixUtils.h>
27 #define EIGENINTERNALC_MAXITER 30
32 #define RSS( a, b ) ( fabs( a ) > fabs( b ) ? \
33 fabs( a )*sqrt( 1.0 + ((b)/(a))*((b)/(a)) ) : \
34 ( (b) == 0.0 ? 0.0 : fabs( b )*sqrt( 1.0 + ((a)/(b))*((a)/(b)) ) ) )
111 UINT4 ij, ik, ki, il, li, kl, lk, ii;
145 for ( i = n - 1; i; i-- ) {
154 for ( k = 0, ik = in; k < i; k++, ik++ )
155 scale += fabs(
m[ik] );
163 for ( k = 0, ik = in; k < i; k++, ik++ ) {
168 y = (
x >= 0 ? -sqrt( H ) : sqrt( H ) );
175 for ( k = 0, ik = in, ki = i; k < i; k++, ik++, ki += n ) {
180 for ( l = 0, kl = k*n, il = in; l <= k; l++, kl++, il++ )
182 for ( l = k + 1, lk = k*( n + 1 ) + n, il = in + k + 1;
183 l < i; l++, lk += n, il++ )
193 for ( k = 0, ik = in; k < i; k++, ik++ ) {
195 o[k] =
y = o[k] - K*
x;
197 for ( l = 0, kl = n*k, il = in; l <= k; l++, kl++, il++ )
198 m[kl] -=
x*o[l] +
y*
m[il];
215 for ( i = 0; i < n; i++ ) {
220 for ( k = 0; k < i; k++ ) {
224 for ( l = 0, il = in, lk = k; l < i; l++, il++, lk += n )
226 for ( l = 0, li = i, lk = k; l < i; l++, li += n, lk += n )
236 for ( k = 0, ik = in, ki = i; k < i; k++, ik++, ki += n )
287 for ( i = 1; i < n; i++ )
292 for ( i = 0; i < n; i++ ) {
300 for ( j = i; j < n - 1; j++ ) {
301 d2 = fabs( d[j] ) + fabs( d[j+1] );
310 x = 0.5*( d[i+1] - d[i] )/o[i];
313 x = d[j] - d[i] + o[i]/(
x +
r );
315 x = d[j] - d[i] + o[i]/(
x -
r );
321 for ( k = j - 1; k >= i; k-- ) {
324 o[k+1] = (
r =
RSS(
x,
y ) );
337 r = ( d[k] -
x )*s + 2.0*c*xx;
338 d[k+1] =
x + ( yy = s*
r );
342 for ( l = 0, lk = k; l < n; l++, lk += n ) {
344 m[lk+1] = s*
m[lk] + c*
y;
345 m[lk] = c*
m[lk] - s*
y;
350 if ( r <= LAL_REAL4_MIN && k >= i )
375 UINT4 ij, ik, ki, il, li, kl, lk, ii;
409 for ( i = n - 1; i; i-- ) {
418 for ( k = 0, ik = in; k < i; k++, ik++ )
419 scale += fabs(
m[ik] );
427 for ( k = 0, ik = in; k < i; k++, ik++ ) {
432 y = (
x >= 0 ? -sqrt( H ) : sqrt( H ) );
439 for ( k = 0, ik = in, ki = i; k < i; k++, ik++, ki += n ) {
444 for ( l = 0, kl = k*n, il = in; l <= k; l++, kl++, il++ )
446 for ( l = k + 1, lk = k*( n + 1 ) + n, il = in + k + 1;
447 l < i; l++, lk += n, il++ )
457 for ( k = 0, ik = in; k < i; k++, ik++ ) {
459 o[k] =
y = o[k] - K*
x;
461 for ( l = 0, kl = n*k, il = in; l <= k; l++, kl++, il++ )
462 m[kl] -=
x*o[l] +
y*
m[il];
479 for ( i = 0; i < n; i++ ) {
484 for ( k = 0; k < i; k++ ) {
488 for ( l = 0, il = in, lk = k; l < i; l++, il++, lk += n )
490 for ( l = 0, li = i, lk = k; l < i; l++, li += n, lk += n )
500 for ( k = 0, ik = in, ki = i; k < i; k++, ik++, ki += n )
551 for ( i = 1; i < n; i++ )
556 for ( i = 0; i < n; i++ ) {
564 for ( j = i; j < n - 1; j++ ) {
565 d2 = fabs( d[j] ) + fabs( d[j+1] );
574 x = 0.5*( d[i+1] - d[i] )/o[i];
577 x = d[j] - d[i] + o[i]/(
x +
r );
579 x = d[j] - d[i] + o[i]/(
x -
r );
585 for ( k = j - 1; k >= i; k-- ) {
588 o[k+1] = (
r =
RSS(
x,
y ) );
601 r = ( d[k] -
x )*s + 2.0*c*xx;
602 d[k+1] =
x + ( yy = s*
r );
606 for ( l = 0, lk = k; l < n; l++, lk += n ) {
608 m[lk+1] = s*
m[lk] + c*
y;
609 m[lk] = c*
m[lk] - s*
y;
614 if ( r <= LAL_REAL8_MIN && k >= i )
#define EIGENINTERNALC_MAXITER
static void LALSSymmetricToTriDiagonal(LALStatus *stat, REAL4Vector *diag, REAL4Array *matrix, REAL4Vector *offDiag)
static void LALSTriDiagonalToDiagonal(LALStatus *stat, REAL4Vector *diag, REAL4Array *matrix, REAL4Vector *offDiag)
static void LALDSymmetricToTriDiagonal(LALStatus *stat, REAL8Vector *diag, REAL8Array *matrix, REAL8Vector *offDiag)
static void LALDTriDiagonalToDiagonal(LALStatus *stat, REAL8Vector *diag, REAL8Array *matrix, REAL8Vector *offDiag)
#define ABORT(statusptr, code, mesg)
#define ENDFAIL(statusptr)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define BEGINFAIL(statusptr)
#define MATRIXUTILSH_MSGEDIM
#define MATRIXUTILSH_MSGEITER
#define MATRIXUTILSH_MSGENUL
void LALDSymmetricEigenVectors(LALStatus *stat, REAL8Vector *values, REAL8Array *matrix)
void LALSSymmetricEigenVectors(LALStatus *stat, REAL4Vector *values, REAL4Array *matrix)
#define LAL_REAL8_EPS
Difference between 1 and the next resolvable REAL8 2^-52.
#define LAL_REAL4_MIN
Smallest normalized REAL4 number 2^-126.
#define LAL_REAL8_MIN
Smallest normalized REAL8 number 2^-1022.
#define LAL_REAL4_EPS
Difference between 1 and the next resolvable REAL4 2^-23.
double REAL8
Double precision real floating-point number (8 bytes).
uint32_t UINT4
Four-byte unsigned integer.
int32_t INT4
Four-byte signed integer.
float REAL4
Single precision real floating-point number (4 bytes).
#define MATRIXUTILSH_ENUL
Unexpected null pointer in arguments.
#define MATRIXUTILSH_EDIM
Bad matrix dimensions.
#define MATRIXUTILSH_EITER
Did not converge after maximum iterations.
void LALDCreateVector(LALStatus *, REAL8Vector **, UINT4)
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
LAL status structure, see The LALStatus structure for more details.
struct tagLALStatus * statusPtr
Pointer to the next node in the list; NULL if this function is not reporting a subroutine error.
Multidimentional array of REAL4, see DATATYPE-Array types for more details.
UINT4Vector * dimLength
Vector of array dimensions.
REAL4 * data
Pointer to the data array.
Vector of type REAL4, see DATATYPE-Vector types for more details.
REAL4 * data
Pointer to the data array.
UINT4 length
Number of elements in array.
Multidimentional array of REAL8, see DATATYPE-Array types for more details.
UINT4Vector * dimLength
Vector of array dimensions.
REAL8 * data
Pointer to the data array.
Vector of type REAL8, see DATATYPE-Vector types for more details.
REAL8 * data
Pointer to the data array.
UINT4 length
Number of elements in array.
UINT4 length
Number of elements in array.
UINT4 * data
Pointer to the data array.