Compiling solve using Wiedemann's method

[Jonathan]

Hi all;
I'm having difficulty compiling a program to use Wiedemann's method.
I'm using version 1.1.5 rc0. I'm specifying a WiedemannTraits
structure as

Method::Wiedemann MW;
MW.singular(WiedemannTraits::SINGULAR);
MW.preconditioner(WiedemannTraits::SPARSE);
MW.certificate(false);
MW.maxTries(10);

I've tried using the solve routine from the solutions interface
similar to the examples:
LinBox::solve(soln, Mat, in, MW);
but when I compile I get a list of over 200 errors.

I tried instantiating the WiedemannSolver explicitly:
WiedemannSolver<Field> ws(F, MW);
ws.solve(Mat, soln, in, cert);
but it gives the same list of errors.

Using a different routine in WiedemannSolver, it compiles fine:
ws.findNullspaceElement(soln, Mat);
However, findNullspaceElement does not seem to apply the
preconditioning method specified in the WiedemannTraits, so it is
unable to compute a solution most of the time with the large sparse
systems I'm working with.

I've included below a short test program that gives the compile
errors.

Any help would be much appreciated.

Thanks,
Jonathan




#include <iostream>
#include <vector>
#include <NTL/ZZ_p.h>
#include <NTL/vec_ZZ.h>
#include <linbox/field/modular.h>
#include <linbox/field/ntl-ZZ_p.h>
#include <linbox/blackbox/sparse.h>
#include <linbox/solutions/solve.h>
#include <linbox/solutions/methods.h>
#include <linbox/field/PIR-ntl-ZZ_p.h>
#include <linbox/field/unparametric.h>

using namespace std;
using namespace NTL;
using namespace LinBox;

//#define USE_NULLSPACE_ROUTINE
#define USE_EXPLICIT_SOLVER

void get_vector(vec_long v, long length) {
long isZero;
for (long i = 0; i < length; i++) {
isZero = RandomBnd(2);
if (isZero) {
v[i] = 0;
} else {
v[i] = RandomBnd(100);
}
}
cout << v << endl;
}

int main() {
typedef UnparametricField<ZZ_p> Field;

ZZ p = to_ZZ(27);

ZZ_p::init(p);
std::string pstring;
stringstream strstream;
strstream << p << ends;
pstring = strstream.str();
Field F(pstring.data());

long numrows = 10, numcols = 12;
vector<Field::Element> soln(numcols);
SparseMatrix<Field> Mat(F,numrows,numcols);
vec_ZZ kervec;
kervec.SetLength(numcols);

ZZ_p p2;
vec_long tempvec;
tempvec.SetLength(numrows);
for (long i = 0; i < numcols; i++) {
get_vector(tempvec, numrows);
for(long j = 0; j < numrows; j++) {
if(tempvec[j] != 0) {
conv(p2, tempvec[j]);
Mat.setEntry(j, i, p2);
}
}
}

Method::Wiedemann MW;
MW.singular(WiedemannTraits::SINGULAR);
MW.preconditioner(WiedemannTraits::SPARSE);
MW.certificate(false);
MW.maxTries(10);

try {
// Compute a vector in the kernel of Mat
#ifdef USE_NULLSPACE_ROUTINE
enum WiedemannSolver<Field>::ReturnStatus rs;
WiedemannSolver<Field> ws(F, MW);
rs = ws.findNullspaceElement(soln, Mat);
if (rs != WiedemannSolver<Field>::OK) {
cerr << "WiedemannSolver returned " << rs << endl;
return -1;
}
for(long i = 0; i < numcols; i++)
kervec[i] = (rep(soln[i]));
#else
vector<Field::Element> in(numrows);
vector<Field::Element> v(numcols);
for (long i = 0; i < numcols; i++)
v[i] = random_ZZ_p();
Mat.apply(in, v);

#ifdef USE_EXPLICIT_SOLVER
vector<Field::Element> cert(numcols);
enum WiedemannSolver<Field>::ReturnStatus rs;
WiedemannSolver<Field> ws(F, MW);
rs = ws.solve(Mat, soln, in, cert);
if (rs != WiedemannSolver<Field>::OK) {
cerr << "WiedemannSolver returned " << rs << endl;
return -1;
}
#else
LinBox::solve(soln, Mat, in, MW);
#endif // USE_EXPLICIT_SOLVER

for(long i = 0; i < numcols; i++)
kervec[i] = (rep(soln[i] - v[i]));

#endif // USE_NULLSPACE_ROUTINE
} catch(LinBox::SolveFailed) {
cout << "caught SolveFailed" << endl;
return -1;
}

cout << "Kernel vector = " << kervec << endl;
return 0;
}

[Clement Pernet]

Hi Jonathan,

The compilation fails, because the field your are using is not supported
by givaro-extension, in order to run Wiedemann algorithm (extensions are
used by the preconditioner).

The way you create a finite field is incorrect:
Unparametric<X> is meant to be a characteristic 0 ring.
And I have the impression that you are trying to use it for the finite
field GF(27).

Can you tell us more about what is you computation domain?

For GF(27), you should use
GivaroExtension F (3,3)

Alternatively, you can also use other base field implementation, eg
Modular<int>, for GF(3)
and create the extension as follows

Modular<int> F(3);
GivaroExtension<Modular<int> > G(F, 3);

Let me know if I am completely wrong about your intensions.

Cheers,
Clément

Jonathan a écrit :

[Jonathan F. Hammell]

Actually, I want to solve over (Z/nZ), where n is a large composite (ZZ
type) with known prime factorization. My plan was to solve the system
modulo each prime factor of n and then use CRT to combine the results.

I haven't found too much documentation for Givaro, so I'm just trying to
extrapolate the constructions from the examples. This probably isn't
the best forum to ask, but are the examples the best source?

Thanks for your help.

Jonathan

--
Jonathan F. Hammell
Graduate Student
Dept. of Computer Science
University of Calgary

+1(403)210-9418
jham...@cpsc.ucalgary.ca

[Jonathan]

Sorry, I was actually looking at the wrong module. The linbox doc
directory under field seems to contain the required information,
wrapping the givaro methods.

Thanks for the GivaroExtension suggestion. I'll look into it some
more.

Jonathan

On Mar 6, 1:19 pm, "Jonathan F. Hammell" <jhamm...@cpsc.ucalgary.ca>
wrote:

> jhamm...@cpsc.ucalgary.ca

[Jonathan]

Okay, I'm now able to compile Wiedemann's method, but I'm once again
getting a segmentation fault each time I try to solve a system.

In this case I'm trying to solve modulo a small prime, so from
previous suggestions I'm using
typedef Modular<int> Field;

I pasted my short example program below.

Jonathan


#include <iostream>
#include <vector>
#include <NTL/ZZ_p.h>
#include <NTL/vec_ZZ.h>
#include <linbox/field/modular.h>
#include <linbox/field/ntl-ZZ_p.h>
#include <linbox/blackbox/sparse.h>
#include <linbox/solutions/solve.h>
#include <linbox/solutions/methods.h>

using namespace std;
using namespace NTL;
using namespace LinBox;

void get_vector(int p, vec_long v, long length) {

long isZero;
for (long i = 0; i < length; i++) {

isZero = RandomBnd(3);

if (isZero) {
v[i] = 0;
} else {

v[i] = RandomBnd(p);

}
}
cout << v << endl;
}

int main() {

typedef Modular<int> Field;

int p = 17;
Field F(p);

long numrows = 10, numcols = 12;
vector<Field::Element> soln(numcols);
SparseMatrix<Field> Mat(F,numrows,numcols);

vec_long kervec;
kervec.SetLength(numcols);

vec_long tempvec;
tempvec.SetLength(numrows);
for (long i = 0; i < numcols; i++) {

get_vector(p, tempvec, numrows);

for(long j = 0; j < numrows; j++) {

Mat.setEntry(j, i, tempvec[j]);

}
}

Method::Wiedemann MW;
MW.singular(WiedemannTraits::SINGULAR);
MW.preconditioner(WiedemannTraits::SPARSE);
MW.certificate(false);
MW.maxTries(10);

try {
// Compute a vector in the kernel of Mat

vector<Field::Element> in(numrows);
vector<Field::Element> v(numcols);
for (long i = 0; i < numcols; i++)

v[i] = RandomBnd(p);
Mat.apply(in, v);

LinBox::solve(soln, Mat, in, MW);

for(long i = 0; i < numcols; i++)

kervec[i] = soln[i] - v[i];

} catch(LinBox::SolveFailed) {
cout << "caught SolveFailed" << endl;
return -1;
}

cout << "Kernel vector = " << kervec << endl;
return 0;
}

> > >> ....

[Clement Pernet]

The sparse matrix that your are building is the 0 matrix, since your
function get_vector does not modify v outside of its scope.
LinBox wiedemann algorithm segfaults with a matrix of rank 0 (to be fixed).

Simply replacing

void get_vector(int p, vec_long v, long length) {

by
void get_vector(int p, vec_long& v, long length) {

makes it work the way you want.

Now I also sugest that you do not set the zero entries, in order to
keep the matrix sparse.

After these modifications, the execution raises a Solve Failed
exception, which I still did not investigated.

However, I found that we do have a function sampling a vector from the
kernel in WiedemannSolver (algorithms/wiedemann.h).

So doing the following
WiedemannSolver<Field> WS(F, MW);
WS.findNullspaceElement (soln, Mat);

for(long i = 0; i < numcols; i++)

kervec[i] = soln[i];
seems to work for me. It finds very often the 0 vector when p is small,
but this gets better for larger p.
A work around could be to loop while you get the 0 vector, although
there might be a more clever technique, which is probably not yet
implemented.

Sorry for not having driven you directly to this function, since I am
not so familiar with this part of the library.

Attached is my version of your code, summarizing the mods I have just
described.

Clément

Jonathan a écrit :

[Jonathan]

What a silly mistake! Sorry about that. Thanks for pointing it out.

It took me awhile, but I was finally able to get a reproducible example
of how it is failing in my actual system. When I use the
findNullspaceElement routine, I occasionally will get a incorrect
answer, i.e. it returns (without error) a vector that is *not* in the
kernel. That is, after calling

ws.findNullspaceElement(soln, Mat);
Mat.apply(test, soln);

the check that test is the zero vector fails.

I've attached the code with a matrix that always return an incorrect
answer on my system with version 1.1.5rc0.

I appreciate the help.