Dear Sir/Madam,
Could you help me with the selftest.c (see the attached file).
I've received different results on sparc and intel for 32-b
Running the GSL self test suite I am seeing the following error when
running
on x86 32-bit:
make[3]: Entering directory `$(@D)/linalg'
FAIL: LQ_update m(5,3) [88]
FAIL: LQ Rank-1 Update [97]
FAIL: test
==================
1 of 1 test failed
==================
The same test works just fine on 32-bit SPARC:
make[3]: Entering directory `$(@D)/linalg'
Completed [233845/233845]
PASS: test
=============
1 test passed
=============
The test failed with the following results:
Expected: 0 1
Got : 6.82121026329696178e-13 1.00000000000056843
i.e. an absolute error in the case of 0 expected is ~ 6.8e-13 (instead
of ~2.5e-13); /2 * 512.0 * 2.2204460492503131e-16
a relative error in the case of 1 expected is ~ 5.7e-13 (instead of
~2.5e-13); /2 * 512.0 * 2.2204460492503131e-16
The implementation in the cases of sparc and intel is different, because
of the different assembler code for sparc and intel.
Rounding also can be different.
Therefore, the absolute or relative error of the Strassen algorithm can
be different on different platforms.
The question is, which max value of error is acceptable.
I do not understand, why the barrier value is chosen
2 * 512.0 * 2.2204460492503131e-16
2.2204460492503131e-16 is a machine eps for double.
Why 2 * 512.0 * eps?
Maybe, this is too strong condition?
The calculating error should depend on the size of source matrices.
Could you explain this to me?
Thsnk you and best wishes,
Elena Ivanova.
typedef unsigned int size_t;
struct gsl_block_struct
{
size_t size;
double *data;
};
typedef struct gsl_block_struct gsl_block;
typedef struct
{
size_t size1;
size_t size2;
size_t tda;
double * data;
gsl_block * block;
int owner;
} gsl_matrix;
gsl_matrix *gsl_matrix_alloc (const size_t n1, const size_t n2);
void gsl_matrix_free (gsl_matrix * m);
inline
double
gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j)
{
return m->data[i * m->tda + j] ;
}
inline
void
gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x)
{
m->data[i * m->tda + j] = x ;
}
int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src);
extern double cos ( double );
extern double sin ( double );
extern double fabs ( double );
typedef struct
{
size_t size;
size_t stride;
double *data;
gsl_block *block;
int owner;
}
gsl_vector;
gsl_vector *gsl_vector_alloc (const size_t n);
void gsl_vector_free (gsl_vector * v);
inline
void
gsl_vector_set (gsl_vector * v, const size_t i, double x)
{
v->data[i * v->stride] = x;
}
void gsl_test (int status, const char *test_description, ...);
int gsl_blas_dger (double alpha,
const gsl_vector * X,
const gsl_vector * Y,
gsl_matrix * A);
int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);
int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
gsl_matrix * Q, gsl_matrix * L);
enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
typedef enum CBLAS_TRANSPOSE CBLAS_TRANSPOSE_t;
int gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA,
double alpha,
const gsl_matrix * A,
const gsl_vector * X,
double beta,
gsl_vector * Y);
int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
const gsl_vector * v, gsl_vector * w);
int gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA,
CBLAS_TRANSPOSE_t TransB,
double alpha,
const gsl_matrix * A,
const gsl_matrix * B,
double beta,
gsl_matrix * C);
int
check (double x, double actual, double eps)
{
if (x == actual)
{
return 0;
}
else if (actual == 0)
{
return fabs(x) > eps;
}
else
{
return (fabs(x - actual)/fabs(actual)) > eps;
}
}
gsl_matrix *
create_general_matrix(unsigned long size1, unsigned long size2)
{
unsigned long i, j;
gsl_matrix * m = gsl_matrix_alloc(size1, size2);
for(i=0; i<size1; i++) {
for(j=0; j<size2; j++) {
gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
}
}
return m;
}
gsl_matrix * m53;
int
test_LQ_update_dim(const gsl_matrix * m, double eps)
{
int s = 0;
unsigned long i,j, M = m->size1, N = m->size2;
gsl_matrix * lq1 = gsl_matrix_alloc(N,M);
gsl_matrix * lq2 = gsl_matrix_alloc(N,M);
gsl_matrix * q1 = gsl_matrix_alloc(M,M);
gsl_matrix * l1 = gsl_matrix_alloc(N,M);
gsl_matrix * q2 = gsl_matrix_alloc(M,M);
gsl_matrix * l2 = gsl_matrix_alloc(N,M);
gsl_vector * d2 = gsl_vector_alloc( ( ( M ) < ( N ) ? ( M ) : ( N ) ));
gsl_vector * u = gsl_vector_alloc(M);
gsl_vector * v = gsl_vector_alloc(N);
gsl_vector * w = gsl_vector_alloc(M);
gsl_matrix_transpose_memcpy(lq1,m);
gsl_matrix_transpose_memcpy(lq2,m);
for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
gsl_blas_dger(1.0, v, u, lq1);
s += gsl_linalg_LQ_decomp(lq2, d2);
s += gsl_linalg_LQ_unpack(lq2, d2, q2, l2);
gsl_blas_dgemv(CblasNoTrans, 1.0, q2, u, 0.0, w);
s += gsl_linalg_LQ_update(q2, l2, v, w);
gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,l2,q2,0.0,lq2);
for(i=0; i<N; i++) {
for(j=0; j<M; j++) {
double s1 = gsl_matrix_get(lq1, i, j);
double s2 = gsl_matrix_get(lq2, i, j);
int foo = check(s1, s2, eps);
s += foo;
}
}
gsl_vector_free(d2);
gsl_vector_free(u);
gsl_vector_free(v);
gsl_vector_free(w);
gsl_matrix_free(lq1);
gsl_matrix_free(lq2);
gsl_matrix_free(q1);
gsl_matrix_free(l1);
gsl_matrix_free(q2);
gsl_matrix_free(l2);
return s;
}
int test_LQ_update(void)
{
int f;
int s = 0;
f = test_LQ_update_dim(m53, 2 * 512.0 * 2.2204460492503131e-16 );
gsl_test(f, " LQ_update m(5,3)");
s += f;
return s;
}
int main(void)
{
m53 = create_general_matrix(5,3);
gsl_test(test_LQ_update(), "LQ Rank-1 Update");
}
