Package: librsb Version: 1.2.0-rc5-3 Severity: important librsb-rc5 (shipped in Debian 9) suffers a few severe bugs leading to numerically wrong results. These were solved in -rc7. I attach a program detecting the main bug and a patch to solve it and the others, by essentially bringing those files to be like in -rc7.
/* This program demonstrates a bug fixed in librsb-1.2-rc7. */
#include <rsb.h>
#include <stdio.h>
int main(const int argc, char * const argv[])
{
rsb_err_t errval = RSB_ERR_NO_ERROR;
int i;
const int N = 2;
double complex x[N], y[N];
const double complex alpha = 1.0;
double complex VA[]={1,1,2};
int IA[]={0,0,1},JA[]={0,1,1},nnz=3;
struct rsb_mtx_t * mtxAp;
// matrix:
// (1,0) (1,0)
// (1,0) (1,0)
// represented as upper hermitian:
// (1,0) (1,0)
// (0,0) (1,0)
// multiplicand vector:
// (0,1)
// (0,0)
// result shall be:
// (0,1)
// (0,1)
errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS);
if(errval != RSB_ERR_NO_ERROR)
goto err;
mtxAp = rsb_mtx_alloc_from_coo_begin(2,
RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX, N, N,
RSB_FLAG_DEFAULT_STORAGE_FLAGS|RSB_FLAG_UPPER_HERMITIAN , &errval);
errval = rsb_mtx_set_vals(mtxAp, VA, IA, JA, nnz, RSB_FLAG_NOFLAGS);
if(errval != RSB_ERR_NO_ERROR)
goto err;
errval = rsb_mtx_alloc_from_coo_end(&mtxAp);
if(errval != RSB_ERR_NO_ERROR)
goto err;
for (i = 0; i < N; i++) x[i] = y[i] = 0.0;
x[0] = I;
errval = rsb_spmv(RSB_TRANSPOSITION_N, &alpha, mtxAp, x, 1, NULL, y, 1);
if(errval != RSB_ERR_NO_ERROR)
goto err;
for (i = 0; i < N; i++)
printf("(%g,%g)\t", creal(y[i]), cimag(y[i]));
printf("\n");
if(cimag(y[1])!=1)
printf("Result seems incorrect -- are you using <
librsb-1.2-rc7?\n");
else
printf("Result seems correct -- are you using >=
librsb-1.2-rc7?\n");
rsb_mtx_free(mtxAp);
errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS);
if(errval != RSB_ERR_NO_ERROR)
goto err;
return 0;
err:
rsb_perror(NULL,errval);
return -1;
}
Index: rsb_util.m4
===================================================================
--- rsb_util.m4 (revision 3505)
+++ rsb_util.m4 (revision 3516)
@@ -2708,6 +2708,7 @@
--ti;
vT[ti] += vB[bi] + vC[ci];
++bi;
+ ++ci;
++ti;
++onz;
}
@@ -2715,16 +2716,16 @@
t0 = ti;
if ( bi<nnzB && ci<nnzC &&
RSB_COO_GT(IB[bi],JB[bi],IC[ci],JC[ci]) )
{
- IT[ti] = IB[bi];
- JT[ti] = JB[bi];
- vT[ti] = vB[bi] + vC[ci];
- ++bi,++ci,++ti;
+ IT[ti] = IC[ci];
+ JT[ti] = JC[ci];
+ vT[ti] = vC[ci];
+ ++ci,++ti;
}
while( bi<nnzB && ci<nnzC &&
RSB_COO_GT(IB[bi],JB[bi],IC[ci],JC[ci]) && ti > 0 &&
RSB_COO_EQ(IC[ci],JC[ci],IT[ti-1],JT[ti-1]) )
{
--ti;
- vT[ti] += vB[bi] + vC[ci];
+ vT[ti] += vC[ci];
++ci;
++ti;
++onz;
@@ -2733,6 +2734,8 @@
if( ci < nnzC && bi < nnzB )
goto again`_'RSB_M4_CHOPSPACES(mtype);
+again`_once_'RSB_M4_CHOPSPACES(mtype):
+
if ( bi<nnzB && ci==nnzC )
{
IT[ti] = IB[bi];
@@ -2761,14 +2764,14 @@
while( ci<nnzC && bi==nnzB && ti > 0 &&
RSB_COO_EQ(IC[ci],JC[ci],IT[ti-1],JT[ti-1]) )
{
--ti;
- IT[ti] = IC[ci];
- JT[ti] = JC[ci];
vT[ti]+= vC[ci];
++ci;
++ti;
++onz;
}
+ if( ci < nnzC || bi < nnzB )
+ goto again`_once_'RSB_M4_CHOPSPACES(mtype);
}
else
#endif
@@ -2825,7 +2828,7 @@
mtypeb*tdst = dst;
ifelse(RSB_M4_AND(RSB_M4_IS_COMPLEX_TYPE(mtypeb)),1,`dnl
if(RSB_DOES_CONJUGATE(transA))
- for(nzi=0;nzi<nnz;++nzi)
RSB_M4_ASSIGN(mtypeb,mtypea,`tdst[nzi]',`RSB_M4_CONJ(`(mtypeb)(alpha*tsrc[nzi])',mtypeb,RSB_M4_TRANS_C,RSB_M4_SYMBOL_UNSYMMETRIC)')
+ for(nzi=0;nzi<nnz;++nzi)
RSB_M4_ASSIGN(mtypeb,mtypea,`tdst[nzi]',`(mtypeb)(alpha*RSB_M4_CONJ(`tsrc[nzi]',mtypeb,RSB_M4_TRANS_C,RSB_M4_SYMBOL_UNSYMMETRIC))')
else
')dnl
for(nzi=0;nzi<nnz;++nzi)
RSB_M4_ASSIGN(mtypeb,mtypea,`tdst[nzi]',`(mtypeb)(alpha*tsrc[nzi])')
Index: rsb_krnl_bcss_macros.m4
===================================================================
--- rsb_krnl_bcss_macros.m4 (revision 3505)
+++ rsb_krnl_bcss_macros.m4 (revision 3516)
@@ -922,9 +922,9 @@
k=fk;
if(k==lk)continue;
j=bindx[k];
- cacc +=
RSB_M4_CONJ(VA[k]*rhs[tcolsu*j*xstride],mtype,ntransposition,k_symmetry);
+ cacc +=
RSB_M4_CONJ(VA[k],mtype,ntransposition,k_symmetry)*rhs[tcolsu*j*xstride];
if(roff!=coff || (j!=i))
-
tout[(tcolsu)*(j)*ystride]+=RSB_M4_CONJ(VA[k]*bt,mtype,ttransposition,k_symmetry);
+
tout[(tcolsu)*(j)*ystride]+=RSB_M4_CONJ(VA[k],mtype,ttransposition,k_symmetry)*bt;
++k;
dnl RSB_M4_SIMPLE_LOOP_UNROLL_2S..
RSB_M4_SIMPLE_LOOP_UNROLL_2S_J(`k',`LI',`fk+1',`lk-1',`dnl
@@ -933,7 +933,7 @@
`const mtype
b_'``''LI`'=rhs[tcolsu*(`j_'``''LI`')*xstride];
`const mtype a_'``''LI`'=VA[k+LI];
`mtype
c_'``''LI`'=RSB_M4_CONJ_SYM(mtype,ttransposition,k_symmetry)( `a_'``''LI)*bt;
-dnl `mtype c_'``''LI`'=RSB_M4_CONJ(( `a_'``''LI *bt
),mtype,transposition,k_symmetry);
+dnl `mtype c_'``''LI`'=RSB_M4_CONJ((
`a_'``''LI),mtype,transposition,k_symmetry) *bt ;
dnl
',`dnl
cacc +=
RSB_M4_CONJ_SYM(mtype,ntransposition,k_symmetry)(`a_'``''LI)*b_``''LI;
@@ -942,9 +942,9 @@
if(k<lk)
{
j=bindx[k];
- cacc +=
RSB_M4_CONJ(VA[k]*rhs[trowsu*j*xstride],mtype,ntransposition,k_symmetry);
+ cacc +=
RSB_M4_CONJ(VA[k],mtype,ntransposition,k_symmetry)*rhs[trowsu*j*xstride];
if(roff!=coff || (j!=i))
-
tout[(tcolsu)*(j)*ystride]+=RSB_M4_CONJ(VA[k]*bt,mtype,ttransposition,k_symmetry);
+
tout[(tcolsu)*(j)*ystride]+=RSB_M4_CONJ(VA[k],mtype,ttransposition,k_symmetry)*bt;
++k;
}
popdef(`ntransposition')dnl
@@ -1047,7 +1047,7 @@
dnl Fixed for Hermitian k_symmetry.
dnl
ifelse(is_diag_d_spsv_kernel,1,`dnl
-
out[trowsu*bri]-=RSB_M4_CONJ(*a*ax_0,mtype,transposition,k_symmetry);
+
out[trowsu*bri]-=RSB_M4_CONJ(*a,mtype,transposition,k_symmetry)*ax_0;
',`dnl
{RSB_M4_KERNEL_FUNCTION_BODY(`row',`rows',b_rows,`column',`columns',b_columns,mtype,,mop,unrolling,RSB_M4_SYMBOL_UNSYMMETRIC)}
')dnl
Index: rsb_libspblas_handle.c
===================================================================
--- rsb_libspblas_handle.c (revision 3505)
+++ rsb_libspblas_handle.c (revision 3516)
@@ -1,6 +1,6 @@
/*
-Copyright (C) 2008-2015 Michele Martone
+Copyright (C) 2008-2017 Michele Martone
This file is part of librsb.
@@ -1527,9 +1527,15 @@
\ingroup gr_internals
*/
rsb_err_t errval = RSB_ERR_NO_ERROR;
+#if 1
+ struct rsb_blas_sparse_matrix_t * bsm = NULL;
+ if( ( bsm = rsb__BLAS_matrix_retrieve(A)) != NULL && bsm->base ==
blas_one_base )
+ i-=1, j-=1;
+ errval = rsb__do_get_coo_element(bsm->mtxAp,v,i,j);
+#else
struct rsb_mtx_t * mtxAp = rsb__BLAS_inner_matrix_retrieve(A);
-
errval = rsb__do_get_coo_element(mtxAp,v,i,j);
+#endif
return RSB_ERROR_TO_BLAS_ERROR(errval);
}
signature.asc
Description: PGP signature
