[eclib] 01/01: Imported Upstream version 20150408
Julien Puydt
julien.puydt at laposte.net
Wed Apr 8 18:48:48 UTC 2015
This is an automated email from the git hooks/post-receive script.
jpuydt-guest pushed a commit to annotated tag upstream/20150408
in repository eclib.
commit 057c210ea336ca42fbd25f9140507eba00ba5def
Author: Julien Puydt <julien.puydt at laposte.net>
Date: Wed Apr 8 20:31:56 2015 +0200
Imported Upstream version 20150408
---
configure.ac | 6 +-
libsrc/eclib/homspace.h | 18 +++-
libsrc/eclib/newforms.h | 8 ++
libsrc/eclib/smat.h | 2 +
libsrc/eclib/splitbase.h | 4 +
libsrc/eclib/xsplit.h | 13 +--
libsrc/eclib/xsplit_data.h | 7 +-
libsrc/homspace.cc | 218 ++++++++++++++++++++++++++++++++++++++++++++-
libsrc/mat.cc | 56 ++++++++----
libsrc/smat.cc | 152 +++++++++++++++++++++++++++++++
libsrc/xsplit.cc | 180 +++++++++++++++++++++++++++----------
libsrc/xsplit_data.cc | 26 ++++--
12 files changed, 604 insertions(+), 86 deletions(-)
diff --git a/configure.ac b/configure.ac
index 9888784..c06041a 100644
--- a/configure.ac
+++ b/configure.ac
@@ -3,7 +3,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.65])
-AC_INIT([eclib], [20150323], [john.cremona at gmail.com])
+AC_INIT([eclib], [20150408], [john.cremona at gmail.com])
AM_INIT_AUTOMAKE([-Wall])
AC_MSG_NOTICE([Configuring eclib...])
AC_CONFIG_SRCDIR([libsrc])
@@ -23,8 +23,8 @@ LT_INIT
# revision to zero and increment age.
# 3. If interfaces were removed (breaks backward compatibility): increment
# current, and set both revision and age to zero.
-LT_CURRENT=1
-LT_REVISION=1
+LT_CURRENT=2
+LT_REVISION=0
LT_AGE=0
AC_SUBST(LT_CURRENT)
AC_SUBST(LT_REVISION)
diff --git a/libsrc/eclib/homspace.h b/libsrc/eclib/homspace.h
index c32aaec..eca6836 100644
--- a/libsrc/eclib/homspace.h
+++ b/libsrc/eclib/homspace.h
@@ -56,9 +56,13 @@ public:
vector<long> eigrange(long i);
long op_prime(int i); // the i'th operator prime for Tp or Wq
mat opmat(int i, int dual, int verb=0);
+ vec opmat_col(int i, int j, int verb=0);
+ mat opmat_cols(int i, const vec& jlist, int verb=0);
mat opmat_restricted(int i,const subspace& s, int dual, int verb=0);
// versions returning an smat:
smat s_opmat(int i,int dual,int verb=0);
+ svec s_opmat_col(int i, int j, int verb=0);
+ smat s_opmat_cols(int i, const vec& jlist, int verb=0);
smat s_opmat_restricted(int i,const ssubspace& s, int dual,int verb=0);
// Extend a dual vector of length rk to one of length nsymb:
@@ -118,26 +122,38 @@ public:
// {return applyop(mlist,m.beta())-applyop(mlist,m.alpha());}
mat calcop(string opname, long p, const matop& mlist, int dual, int display=0) const;
+ vec calcop_col(string opname, long p, int j, const matop& mlist, int display=0) const;
+ mat calcop_cols(string opname, long p, const vec& jlist, const matop& mlist, int display=0) const;
mat calcop_restricted(string opname, long p, const matop& mlist, const subspace& s, int dual, int display=0) const;
smat s_calcop(string opname, long p, const matop& mlist, int dual, int display=0) const;
+ svec s_calcop_col(string opname, long p, int j, const matop& mlist, int display=0) const;
+ smat s_calcop_cols(string opname, long p, const vec& jlist, const matop& mlist, int display=0) const;
smat s_calcop_restricted(string opname, long p, const matop& mlist, const ssubspace& s, int dual, int display=0) const;
public:
mat heckeop(long p, int dual, int display=0) const;
+ vec heckeop_col(long p, int j, int display=0) const;
+ mat heckeop_cols(long p, const vec& jlist, int display=0) const;
mat heckeop_restricted(long p, const subspace& s, int dual, int display=0) const;
smat s_heckeop(long p, int dual, int display=0) const;
+ svec s_heckeop_col(long p, int j, int display=0) const;
+ smat s_heckeop_cols(long p, const vec& jlist, int display=0) const;
smat s_heckeop_restricted(long p, const ssubspace& s, int dual, int display=0) const;
mat newheckeop(long p, int dual, int display=0) const;
mat wop(long q, int dual, int display=0) const;
smat s_wop(long q, int dual, int display=0) const;
mat fricke(int dual, int display=0) const;
- mat conj(int dual,int display=0) const;
+ mat conj(int dual, int display=0) const;
+ vec conj_col(int j, int display=0) const;
+ mat conj_cols(const vec& jlist, int display=0) const;
mat conj_restricted(const subspace& s, int dual,int display=0) const;
smat s_conj(int dual, int display=0) const;
+ svec s_conj_col(int j, int display=0) const;
+ smat s_conj_cols(const vec& jlist, int display=0) const;
smat s_conj_restricted(const ssubspace& s, int dual, int display=0) const;
vec maninvector(long p) const;
vec projmaninvector(long p) const;
diff --git a/libsrc/eclib/newforms.h b/libsrc/eclib/newforms.h
index 2f14002..73a8fc9 100644
--- a/libsrc/eclib/newforms.h
+++ b/libsrc/eclib/newforms.h
@@ -113,10 +113,18 @@ private:
// already defined.
mat opmat(int i, int d, int v=0)
{return h1->opmat(i,d,v);}
+ vec opmat_col(int i, int j, int v=0)
+ {return h1->opmat_col(i,j,v);}
+ mat opmat_cols(int i, const vec& jlist, int v=0)
+ {return h1->opmat_cols(i,jlist,v);}
mat opmat_restricted(int i, const subspace& s, int d, int v=0)
{return h1->opmat_restricted(i,s,d,v);}
smat s_opmat(int i, int d, int v=0)
{return h1->s_opmat(i,d,v);}
+ svec s_opmat_col(int i, int j, int v=0)
+ {return h1->s_opmat_col(i,j,v);}
+ smat s_opmat_cols(int i, const vec& jlist, int v=0)
+ {return h1->s_opmat_cols(i,jlist,v);}
smat s_opmat_restricted(int i, const ssubspace& s, int d, int v=0)
{return h1->s_opmat_restricted(i,s,d,v);}
long matdim(void) {return h1->dimension;}
diff --git a/libsrc/eclib/smat.h b/libsrc/eclib/smat.h
index 6891e81..f21c46b 100644
--- a/libsrc/eclib/smat.h
+++ b/libsrc/eclib/smat.h
@@ -78,9 +78,11 @@ public:
friend svec operator* ( const smat& A, const svec& v );
friend svec operator* ( const svec& v, const smat& A );
friend svec mult_mod_p( const smat& A, const svec& v, const scalar& p );
+ friend vec mult_mod_p( const smat& A, const vec& v, const scalar& p );
friend svec mult_mod_p( const svec& v, const smat& A, const scalar& p );
friend smat operator* ( const smat& A, const smat& B );
friend smat mult_mod_p ( const smat& A, const smat& B, const scalar& p );
+ friend smat mult_mod_p_flint ( const smat& A, const smat& B, const scalar& p );
friend smat transpose(const smat&);
friend int operator==(const smat&, const smat&);
// Equality mod p:
diff --git a/libsrc/eclib/splitbase.h b/libsrc/eclib/splitbase.h
index 54d3e02..77262b3 100644
--- a/libsrc/eclib/splitbase.h
+++ b/libsrc/eclib/splitbase.h
@@ -50,8 +50,12 @@
class splitter_base {
public:
virtual mat opmat(int,int,int=0) = 0;
+ virtual vec opmat_col(int i, int j, int verb=0) = 0;
+ virtual mat opmat_cols(int i, const vec& jlist, int verb=0) = 0;
virtual mat opmat_restricted(int,const subspace& s, int,int=0) = 0;
virtual smat s_opmat(int,int,int=0) = 0;
+ virtual svec s_opmat_col(int i, int j, int verb=0) = 0;
+ virtual smat s_opmat_cols(int i, const vec& jlist, int verb=0) = 0;
virtual smat s_opmat_restricted(int,const ssubspace& s, int, int=0) = 0;
virtual long matdim(void) = 0;
virtual long matden(void) = 0;
diff --git a/libsrc/eclib/xsplit.h b/libsrc/eclib/xsplit.h
index 67af30a..efb1699 100644
--- a/libsrc/eclib/xsplit.h
+++ b/libsrc/eclib/xsplit.h
@@ -78,17 +78,17 @@ class form_finder {
vector< vector<long> > gaplist; // Vector to hold all (sub)eiglists
vector<vec> gbplus, gbminus; // Vector to hold all bplus/bminus
- int* havemat;
- vector<string> opfilenames; // Temp filenames
-
ff_data* root; // Always points to root data node
-
+
void make_opmat(long i, ff_data &data); // Puts it in the_opmat
void make_submat(ff_data &data);
+ smat make_nested_submat(long ip, ff_data &data);
void go_down(ff_data &data, long eig, int last=0);
void go_up( ff_data &data );
void make_basis(ff_data &data);
- vec getbasis1(const ssubspace* s); // Assuming dim(s)=1, get basis vector
+ vec make_basis1(ff_data &data); // current space has dim. 1
+ vec make_basis2(ff_data &data, const vec& v); // current space has dim. >1 and
+ // relative basis vector v
#ifdef ECLIB_MULTITHREAD
threadpool pool; // Job queue
@@ -96,4 +96,7 @@ class form_finder {
#endif
};
+vec getbasis1(const ssubspace* s); // Assuming dim(s)=1, get basis vector
+vec lift(const vec& v); // Lift basis vector
+
#endif
diff --git a/libsrc/eclib/xsplit_data.h b/libsrc/eclib/xsplit_data.h
index d26b58b..36e3e5e 100644
--- a/libsrc/eclib/xsplit_data.h
+++ b/libsrc/eclib/xsplit_data.h
@@ -66,7 +66,8 @@ class ff_data {
// Getters (to maintain consistency)
nodestatus status(); // Return status of current node
- ssubspace* submats(); // Return parent subspace
+ ssubspace* abs_space(); // Return parent absolute subspace
+ ssubspace* rel_space(); // Return parent relative subspace
long depth(); // Return current depth
long subdim(); // Return subdimension
long eig(); // Return associated eigenvalue
@@ -111,13 +112,15 @@ class ff_data {
long eigenvalue_; // Corresponding eigenvalue
vector< long > eiglist_; // Sequence of eigenvalues leading to current
vec bplus_, bminus_;
- ssubspace* nest_; // Current subspace (dynamically created)
+ ssubspace* abs_space_; // Current absolute subspace (dynamically created)
+ ssubspace* rel_space_; // Current relative subspace (dynamically created)
smat conjmat_; // Used only when plus==0 and bigmats==1
smat the_opmat_;
smat submat_;
ff_data* parent_; // Pointer to parent data node
vector<ff_data*> children_; // Pointers to corresponding data nodes
+ ff_data* child_; // Pointer to favoured child
int numChildren_; // Store number of children
vector<childstatus> completedChildren_; // Flags for child completion
int submatUsage_; // Counter for submat
diff --git a/libsrc/homspace.cc b/libsrc/homspace.cc
index f2c5569..5fe6a15 100644
--- a/libsrc/homspace.cc
+++ b/libsrc/homspace.cc
@@ -739,7 +739,57 @@ mat homspace::calcop(string opname, long p, const matop& mlist,
}
return m;
}
-
+
+vec homspace::calcop_col(string opname, long p, int j, const matop& mlist,
+ int display) const
+{
+ // j counts from 1
+ vec colj = applyop(mlist,freemods[j-1]).as_vec();
+ if (display)
+ {
+ cout << "Image of "<<j<<"-th generator under " << opname << "(" << p << ") = "
+ << colj << endl;
+ }
+ return colj;
+}
+
+mat homspace::calcop_cols(string opname, long p, const vec& jlist, const matop& mlist,
+ int display) const
+{
+ int i, d = dim(jlist);
+ mat m(d,rk);
+ for (i=1; i<=d; i++)
+ {
+ m.setrow(i, applyop(mlist,freemods[jlist[i]-1]).as_vec());
+ }
+ return m;
+}
+
+svec homspace::s_calcop_col(string opname, long p, int j, const matop& mlist,
+ int display) const
+{
+ // j counts from 1
+ svec colj = applyop(mlist,freemods[j-1]);
+ if (display)
+ {
+ cout << "Image of "<<j<<"-th generator under " << opname << "(" << p << ") = "
+ << colj.as_vec() << endl;
+ }
+ return colj;
+}
+
+smat homspace::s_calcop_cols(string opname, long p, const vec& jlist, const matop& mlist,
+ int display) const
+{
+ int i, d = dim(jlist);
+ smat m(d,rk);
+ for (i=1; i<=d; i++)
+ {
+ m.setrow(i, applyop(mlist,freemods[jlist[i]-1]));
+ }
+ return m;
+}
+
smat homspace::s_calcop(string opname, long p, const matop& mlist,
int dual, int display) const
{
@@ -782,6 +832,34 @@ smat homspace::s_heckeop(long p, int dual, int display) const
return s_calcop(name,p,matlist,dual,display);
}
+vec homspace::heckeop_col(long p, int j, int display) const
+{
+ matop matlist(p,modulus);
+ string name = ((modulus%p) ? T_opname : W_opname);
+ return calcop_col(name,p,j,matlist,display);
+}
+
+mat homspace::heckeop_cols(long p, const vec& jlist, int display) const
+{
+ matop matlist(p,modulus);
+ string name = ((modulus%p) ? T_opname : W_opname);
+ return calcop_cols(name,p,jlist,matlist,display);
+}
+
+svec homspace::s_heckeop_col(long p, int j, int display) const
+{
+ matop matlist(p,modulus);
+ string name = ((modulus%p) ? T_opname : W_opname);
+ return s_calcop_col(name,p,j,matlist,display);
+}
+
+smat homspace::s_heckeop_cols(long p, const vec& jlist, int display) const
+{
+ matop matlist(p,modulus);
+ string name = ((modulus%p) ? T_opname : W_opname);
+ return s_calcop_cols(name,p,jlist,matlist,display);
+}
+
mat homspace::heckeop_restricted(long p,const subspace& s,
int dual, int display) const
{
@@ -848,6 +926,50 @@ mat homspace::conj(int dual, int display) const
return m;
}
+vec homspace::conj_col(int j, int display) const
+{
+ // j counts from 1
+ symb s = symbol(freegens[j-1]);
+ vec colj = coords_cd(-s.cee(),s.dee()).as_vec();
+ if (display) cout << "Column "<<j<<" of matrix of conjugation = " << colj << endl;
+ return colj;
+}
+
+mat homspace::conj_cols(const vec& jlist, int display) const
+{
+ int d = dim(jlist);
+ mat m(d,rk);
+ int i;
+ for (i=1; i<=d; i++)
+ {
+ symb s = symbol(freegens[jlist[i]-1]);
+ m.setrow(i, coords_cd(-s.cee(),s.dee()).as_vec());
+ }
+ return m;
+}
+
+svec homspace::s_conj_col(int j, int display) const
+{
+ // j counts from 1
+ symb s = symbol(freegens[j-1]);
+ svec colj = coords_cd(-s.cee(),s.dee());
+ if (display) cout << "Column "<<j<<" of matrix of conjugation = " << colj.as_vec() << endl;
+ return colj;
+}
+
+smat homspace::s_conj_cols(const vec& jlist, int display) const
+{
+ int d = dim(jlist);
+ smat m(d,rk);
+ int i;
+ for (i=1; i<=d; i++)
+ {
+ symb s = symbol(freegens[jlist[i]-1]);
+ m.setrow(i, coords_cd(-s.cee(),s.dee()));
+ }
+ return m;
+}
+
smat homspace::s_conj(int dual, int display) const
{
smat m(rk,rk);
@@ -942,6 +1064,100 @@ mat homspace::opmat(int i, int dual, int v)
else return heckeop(p,dual,0); // Automatically chooses W or T
}
+vec homspace::opmat_col(int i, int j, int v)
+{
+ if(i==-1) return conj_col(j,v);
+ if((i<0)||(i>=nap))
+ {
+ cout<<"Error in homspace::opmat_col(): called with i = " << i << endl;
+ ::abort();
+ return vec(dimension); // shouldn't happen
+ }
+ long p = op_prime(i);
+ if(v)
+ {
+ cout<<"Computing col "<<j<<" of " << ((::divides(p,modulus)) ? W_opname : T_opname)
+ <<"("<<p<<")..."<<flush;
+ }
+ vec ans = heckeop_col(p,j,0); // Automatically chooses W or T
+ if(v)
+ {
+ cout<<"done."<<endl;
+ }
+ return ans;
+}
+
+mat homspace::opmat_cols(int i, const vec& jlist, int v)
+{
+ if(i==-1) return conj_cols(jlist,v);
+ int d = dim(jlist);
+ if((i<0)||(i>=nap))
+ {
+ cout<<"Error in homspace::opmat_cols(): called with i = " << i << endl;
+ ::abort();
+ return mat(d,rk); // shouldn't happen
+ }
+ long p = op_prime(i);
+ if(v)
+ {
+ cout<<"Computing "<<d<<" cols of " << ((::divides(p,modulus)) ? W_opname : T_opname)
+ <<"("<<p<<")..."<<flush;
+ }
+ mat ans = heckeop_cols(p,jlist,0); // Automatically chooses W or T
+ if(v)
+ {
+ cout<<"done."<<endl;
+ }
+ return ans;
+}
+
+svec homspace::s_opmat_col(int i, int j, int v)
+{
+ if(i==-1) return s_conj_col(j,v);
+ if((i<0)||(i>=nap))
+ {
+ cout<<"Error in homspace::opmat_col(): called with i = " << i << endl;
+ ::abort();
+ return svec(dimension); // shouldn't happen
+ }
+ long p = op_prime(i);
+ if(v)
+ {
+ cout<<"Computing col "<<j<<" of " << ((::divides(p,modulus)) ? W_opname : T_opname)
+ <<"("<<p<<")..."<<flush;
+ }
+ svec ans = s_heckeop_col(p,j,0); // Automatically chooses W or T
+ if(v)
+ {
+ cout<<"done."<<endl;
+ }
+ return ans;
+}
+
+smat homspace::s_opmat_cols(int i, const vec& jlist, int v)
+{
+ if(i==-1) return s_conj_cols(jlist,v);
+ int d = dim(jlist);
+ if((i<0)||(i>=nap))
+ {
+ cout<<"Error in homspace::opmat_col(): called with i = " << i << endl;
+ ::abort();
+ return smat(d,rk); // shouldn't happen
+ }
+ long p = op_prime(i);
+ if(v)
+ {
+ cout<<"Computing "<<d<<" cols of " << ((::divides(p,modulus)) ? W_opname : T_opname)
+ <<"("<<p<<")..."<<flush;
+ }
+ smat ans = s_heckeop_cols(p,jlist,0); // Automatically chooses W or T
+ if(v)
+ {
+ cout<<"done."<<endl;
+ }
+ return ans;
+}
+
smat homspace::s_opmat(int i, int dual, int v)
{
if(i==-1) return s_conj(dual,v);
diff --git a/libsrc/mat.cc b/libsrc/mat.cc
index 9e1392b..2c560b2 100644
--- a/libsrc/mat.cc
+++ b/libsrc/mat.cc
@@ -1484,12 +1484,16 @@ mat echmodp(const mat& entries, vec& pcols, vec& npcols,
#undef mod_mat_init
#undef mod_mat_clear
#undef mod_mat_entry
+#undef mod_mat_nrows
+#undef mod_mat_ncols
#undef mod_mat_rref
#define uscalar hlimb_t // unsigned int
#define mod_mat hmod_mat_t // uses unsigned ints
#define mod_mat_init hmod_mat_init
#define mod_mat_clear hmod_mat_clear
#define mod_mat_entry hmod_mat_entry
+#define mod_mat_nrows hmod_mat_nrows
+#define mod_mat_ncols hmod_mat_ncols
#define mod_mat_rref hmod_mat_rref
#else
#include "flint/nmod_mat.h"
@@ -1498,12 +1502,16 @@ mat echmodp(const mat& entries, vec& pcols, vec& npcols,
#undef mod_mat_init
#undef mod_mat_clear
#undef mod_mat_entry
+#undef mod_mat_nrows
+#undef mod_mat_ncols
#undef mod_mat_rref
#define uscalar mp_limb_t // unsigned long
#define mod_mat nmod_mat_t // uses unsigned longs
#define mod_mat_init nmod_mat_init
#define mod_mat_clear nmod_mat_clear
#define mod_mat_entry nmod_mat_entry
+#define mod_mat_nrows nmod_mat_nrows
+#define mod_mat_ncols nmod_mat_ncols
#define mod_mat_rref nmod_mat_rref
#endif
#include "flint/profiler.h"
@@ -1516,7 +1524,7 @@ mat echmodp(const mat& entries, vec& pcols, vec& npcols,
// present, which should have been determined by the configure script.
// The unsigned scalar types are #define'd as uscalar.
-mat ref_via_flint(const mat& M, scalar pr)
+void mod_mat_from_mat(mod_mat& A, const mat& M, scalar pr)
{
long nr=M.nrows(), nc=M.ncols();
long i, j;
@@ -1525,11 +1533,33 @@ mat ref_via_flint(const mat& M, scalar pr)
uscalar mod = (uscalar)pr;
// create flint matrix copy of M:
- mod_mat A;
mod_mat_init(A, nr, nc, mod);
for(i=0; i<nr; i++)
for(j=0; j<nc; j++)
- mod_mat_entry(A,i,j) = M(i+1,j+1);
+ mod_mat_entry(A,i,j) = (uscalar)posmod(M(i+1,j+1),pr);
+}
+
+mat mat_from_mod_mat(const mod_mat& A, scalar a) // scalar just to fix return type
+{
+ long nr=mod_mat_nrows(A), nc=mod_mat_ncols(A);
+
+ // create matrix copy of A:
+ mat M(nr, nc);
+ long i, j;
+ for(i=0; i<nr; i++)
+ for(j=0; j<nc; j++)
+ M(i+1,j+1) = mod_mat_entry(A,i,j);
+ return M;
+}
+
+mat ref_via_flint(const mat& M, scalar pr)
+{
+ long nr=M.nrows(), nc=M.ncols();
+ long i, j;
+
+ // create flint matrix copy of M:
+ mod_mat A;
+ mod_mat_from_mat(A,M,pr);
// reduce A to rref:
#ifdef TRACE_FLINT_RREF
@@ -1544,11 +1574,9 @@ mat ref_via_flint(const mat& M, scalar pr)
#endif
// copy back to a new matrix for return:
- mat ans(nr,nc);
- for(i=0; i<nr; i++)
- for(j=0; j<nc; j++)
- ans(i+1,j+1) = mod_mat_entry(A,i,j);
+ mat ans = mat_from_mod_mat(A, pr);
+ // clear the flint matrix and return:
mod_mat_clear(A);
return ans;
}
@@ -1572,15 +1600,9 @@ mat ref_via_flint(const mat& M, vec& pcols, vec& npcols,
// cout << "Size of uscalar = "<<8*sizeof(uscalar)<<" bits"<<endl;
#endif
- // copy of the modulus
- uscalar mod = (uscalar)pr;
-
// create flint matrix copy of M:
mod_mat A;
- mod_mat_init(A, nr, nc, mod);
- for(i=0; i<nr; i++)
- for(j=0; j<nc; j++)
- mod_mat_entry(A,i,j) = (uscalar)posmod(M(i+1,j+1),pr);
+ mod_mat_from_mat(A,M,pr);
#ifdef TRACE_FLINT_RREF
timeit_t t;
@@ -1618,11 +1640,9 @@ mat ref_via_flint(const mat& M, vec& pcols, vec& npcols,
}
// copy back to a new matrix for return:
- mat ans(rk,nc);
- for(i=0; i<rk; i++)
- for(j=0; j<nc; j++)
- ans(i+1,j+1) = mod_mat_entry(A,i,j);
+ mat ans = mat_from_mod_mat(A,pr).slice(rk,nc);
+ // clear the flint matrix and return:
mod_mat_clear(A);
return ans;
}
diff --git a/libsrc/smat.cc b/libsrc/smat.cc
index 3d55194..3da9279 100644
--- a/libsrc/smat.cc
+++ b/libsrc/smat.cc
@@ -566,6 +566,36 @@ svec mult_mod_p( const svec& v, const smat& A, const scalar& p )
return prod;
}
+svec mult_mod_p( const smat& A, const svec& v, const scalar& p )
+{
+ if( v.d != A.ncols() )
+ {
+ cout << "incompatible sizes in A*v\n";
+ cout << "Dimensions "<<dim(A)<<" and "<<v.d<<endl;
+ abort();
+ }
+ svec w(A.nrows());
+ int i;
+ for(i=1; i<=A.nrows(); i++)
+ w.set(i,dotmodp(A.row(i),v,p));
+ return w;
+}
+
+vec mult_mod_p( const smat& A, const vec& v, const scalar& p )
+{
+ if( dim(v) != A.ncols() )
+ {
+ cout << "incompatible sizes in A*v\n";
+ cout << "Dimensions "<<dim(A)<<" and "<<dim(v)<<endl;
+ abort();
+ }
+ vec w(A.nrows());
+ int i;
+ for(i=1; i<=A.nrows(); i++)
+ w.set(i,dotmodp(A.row(i),v,p));
+ return w;
+}
+
#if(1)
smat operator* ( const smat& A, const smat& B )
{
@@ -977,3 +1007,125 @@ int liftmats_chinese(const smat& m1, scalar pr1, const smat& m2, scalar pr2,
}
return 1;
}
+
+// Possible FLINT_LEVEL values are as follows.
+//
+// 0: no FLINT support (or a version <2.3)
+// 1: support for 64-bit nmod_mat (standard from version 2.3)
+// 2: support for 32-bit hmod_mat (non-standard, from version 2.3)
+//
+// The configure script should have detected whether the functions
+// nmod_mat_rref and/or hmod_mat_rref are available to be used here.
+//
+#if FLINT_LEVEL!=0
+//#define TRACE_FLINT_RREF
+#include "flint/fmpz.h"
+#if (SCALAR_OPTION==1)&&(FLINT_LEVEL==2)
+#include "flint/hmod_mat.h"
+#undef uscalar
+#undef mod_mat
+#undef mod_mat_init
+#undef mod_mat_clear
+#undef mod_mat_entry
+#undef mod_mat_nrows
+#undef mod_mat_ncols
+#undef mod_mat_rref
+#undef mod_mat_mul
+#define uscalar hlimb_t // unsigned int
+#define mod_mat hmod_mat_t // uses unsigned ints
+#define mod_mat_init hmod_mat_init
+#define mod_mat_clear hmod_mat_clear
+#define mod_mat_entry hmod_mat_entry
+#define mod_mat_nrows hmod_mat_nrows
+#define mod_mat_ncols hmod_mat_ncols
+#define mod_mat_rref hmod_mat_rref
+#define mod_mat_mul hmod_mat_mul
+#else
+#include "flint/nmod_mat.h"
+#undef uscalar
+#undef mod_mat
+#undef mod_mat_init
+#undef mod_mat_clear
+#undef mod_mat_entry
+#undef mod_mat_nrows
+#undef mod_mat_ncols
+#undef mod_mat_rref
+#undef mod_mat_mul
+#define uscalar mp_limb_t // unsigned long
+#define mod_mat nmod_mat_t // uses unsigned longs
+#define mod_mat_init nmod_mat_init
+#define mod_mat_clear nmod_mat_clear
+#define mod_mat_entry nmod_mat_entry
+#define mod_mat_nrows nmod_mat_nrows
+#define mod_mat_ncols nmod_mat_ncols
+#define mod_mat_rref nmod_mat_rref
+#define mod_mat_mul nmod_mat_mul
+#endif
+#include "flint/profiler.h"
+#include "eclib/timer.h"
+
+// FLINT has two types for modular matrices: standard in FLINT-2.3 has
+// nmod_mat_t with entries of type mp_limb_t (unsigned long);
+// non-standard (in an optional branch) is hmod_mat_t, with entries
+// hlimb_t (unsigned int). We use the former when scalar=long and the
+// latter when scalar=int, provided that the optional functions are
+// present, which should have been determined by the configure script.
+// The unsigned scalar types are #define'd as uscalar.
+
+void mod_mat_from_smat(mod_mat& A, const smat& M, scalar pr)
+{
+ long nr=M.nrows(), nc=M.ncols();
+ long i, j;
+
+ // copy of the modulus for FLINT
+ uscalar mod = (uscalar)pr;
+
+ // create flint matrix copy of M:
+ mod_mat_init(A, nr, nc, mod);
+ for(i=0; i<nr; i++)
+ for(j=0; j<nc; j++)
+ mod_mat_entry(A,i,j) = (uscalar)posmod(M.elem(i+1,j+1),pr);
+}
+
+smat smat_from_mod_mat(const mod_mat& A, const scalar& p) //scalar just to fix return type
+{
+ long nr=mod_mat_nrows(A), nc=mod_mat_ncols(A);
+
+ // create matrix copy of A:
+ smat M(nr, nc);
+ long i, j;
+ for(i=0; i<nr; i++)
+ {
+ svec rowi(nc);
+ for(j=0; j<nc; j++)
+ rowi.set(j+1, mod_mat_entry(A,i,j));
+ M.setrow(i+1,rowi);
+ }
+ return M;
+}
+
+smat mult_mod_p_flint ( const smat& A, const smat& B, const scalar& pr )
+{
+ if( A.ncols() != B.nrows() )
+ {
+ cerr << "incompatible smats in operator *\n";
+ abort();
+ }
+ mod_mat AA, BB, CC;
+ mod_mat_from_smat(AA,A,pr);
+ mod_mat_from_smat(BB,B,pr);
+ mod_mat_init(CC, A.nrows(), B.ncols(), pr);
+ // timer T;
+ // T.start();
+ // mod_mat_mul(CC,AA,BB);
+ // T.stop();
+ // cout<<"mult_mod_p_flint time (size "<<dim(A)<<"x"<<dim(B)<<"): ";
+ // T.show();
+ smat C = smat_from_mod_mat(CC, pr);
+ mod_mat_clear(AA);
+ mod_mat_clear(BB);
+ mod_mat_clear(CC);
+ return C;
+}
+
+#endif
diff --git a/libsrc/xsplit.cc b/libsrc/xsplit.cc
index 63a8af0..c850f24 100644
--- a/libsrc/xsplit.cc
+++ b/libsrc/xsplit.cc
@@ -29,7 +29,6 @@
#define ECLIB_RECURSION_DIM_LIMIT 5821
#include <eclib/logger.h>
#include <eclib/xsplit.h>
-
#include <eclib/smatrix_elim.h>
subspace sparse_combine(const subspace& s1, const subspace& s2);
mat sparse_restrict(const mat& m, const subspace& s);
@@ -70,6 +69,7 @@ form_finder::~form_finder(void) {
delete root;
}
+// This is only used when bigmats==1 and we compute opmats on the entire ambient space:
void form_finder::make_opmat(long i, ff_data &data) {
data.the_opmat_ = h -> s_opmat(i,dual,verbose);
}
@@ -81,25 +81,73 @@ void form_finder::make_submat( ff_data &data ) {
if( bigmats ) {
// fetch the_opmat from file, or compute
make_opmat(depth,data);
-
+
if( depth == 0 ) data.submat_ = data.the_opmat_;
else {
ECLOG(1) << "restricting the_opmat to subspace...";
- data.submat_ = restrict_mat(data.the_opmat_,*(data.nest_));
+ data.submat_ = restrict_mat(data.the_opmat_,*(data.abs_space_));
ECLOG(1) << "done." << endl;
}
-
+
data.the_opmat_ = smat(0,0); // releases its space
}
else {
- if( data.submat_.nrows() == 0 ) {
- if( depth == 0 ) data.submat_ = h -> s_opmat(depth,1,verbose);
- else data.submat_ = h -> s_opmat_restricted(depth,*(data.nest_),1,verbose);
- }
+ if( data.submat_.nrows() == 0 ) // else we have it already
+ {
+ if( depth == 0 )
+ data.submat_ = h -> s_opmat(depth,1,verbose);
+ else
+ {
+ //data.submat_ = h -> s_opmat_restricted(depth,*(data.abs_space_),1,verbose);
+ data.submat_ = make_nested_submat(depth,data);
+ }
+ }
}
}
/**
+ * make_nested_submat()
+ *
+ * Computes and returns the submat -- new nested version
+ */
+
+smat form_finder::make_nested_submat(long ip, ff_data &data)
+{
+ long depth = data.depth_; // current depth
+ long subdim = data.subdim_; // current dimension
+ ff_data *d = &data; // Pointer to nodes
+ int i, j, level;
+
+ ECLOG(1) << "Computing operator of size " << subdim
+ << " at depth " << depth << "..." << flush;
+
+ // first we go up the chain, composing pivotal indices
+
+ vec jlist = iota((scalar)subdim);
+ smat b = d->rel_space_->bas();
+ level = depth;
+ while (level--)
+ {
+ ECLOG(2) << "["<<level<<"]" << flush;
+ jlist = d->rel_space_->pivs()[jlist];
+ d->parent_->child_ = d;
+ d = d->parent_;
+ if(level) b = mult_mod_p(d->rel_space_->bas(), b, MODULUS);
+ //if(level) b = mult_mod_p_flint(d->rel_space_->bas(), b, MODULUS);
+ }
+
+ // now compute the matrix of images of the j'th generator for j in jlist
+ ECLOG(2) << " basis done... " << flush;
+ smat m = h -> s_opmat_cols(ip, jlist, 0);
+ ECLOG(2) << "opmat done... " << flush;
+ m = mult_mod_p(m,b,MODULUS);
+ //m = mult_mod_p_flint(m,b,MODULUS);
+ ECLOG(1) <<"done."<<endl;
+ return m;
+}
+
+
+/**
* go_down()
*
* Initiates creation of new subspace; data stored in new
@@ -121,8 +169,6 @@ void form_finder::go_down(ff_data &data, long eig, int last) {
ECLOG(1) << "Increasing depth to " << depth+1 << ", "
<< "trying eig = " << eig << "..."
<< "after scaling, eig = " << eig2 << "..." << endl;
- // if(depth) eig2*= denom(*nest[depth]); // else latter is 1 anyway
- // ^ data.nest_
ssubspace s(0);
vector<int> submat_dim = dim(data.submat_);
@@ -131,10 +177,10 @@ void form_finder::go_down(ff_data &data, long eig, int last) {
ECLOG(1) << "Using sparse elimination (size = [ "
<< submat_dim_ss.str() << "], density ="
- << density(data.submat_) << ")..." << endl;
+ << density(data.submat_) << ")..." << flush;
ECLOG(3) << "submat = " << data.submat_ << flush;
- s = eigenspace(data.submat_,eig2);
+ s = eigenspace(data.submat_,eig2); // the relative eigenspace
// Increment data usage counter for parent node
data.increaseSubmatUsage();
@@ -148,16 +194,19 @@ void form_finder::go_down(ff_data &data, long eig, int last) {
// data.submat_ = smat(0,0);
//}
- ECLOG(1) << "done (dim = " << dim(s) << "), combining subspaces..." << flush;
-
- if( depth == 0 ) child -> nest_ = new ssubspace(s);
- else child -> nest_ = new ssubspace(combine( *(data.nest_),s ));
-
- ECLOG(1) << "done." << endl;
+ ECLOG(1) << "done (dim = " << dim(s) << ")"<<endl;
+ // ECLOG(1) << ", combining subspaces..." << flush;
+
+ child -> rel_space_ = new ssubspace(s);
+ // if( depth == 0 )
+ // child -> abs_space_ = new ssubspace(s);
+ // else
+ // child -> abs_space_ = new ssubspace(combine( *(data.abs_space_),s ));
+ // ECLOG(1) << "done." << endl;
depth++; // Local depth increment (does not effect data nodes)
- child -> subdim_ = dim( *(child -> nest_) );
+ child -> subdim_ = dim( *(child -> rel_space_) );
ECLOG(1) << "Eigenvalue " << eig
<< " has multiplicity " << child -> subdim_ << endl;
@@ -212,20 +261,24 @@ void form_finder::make_basis( ff_data &data ) {
if(plusflag) {
// must treat separately since we did not
- // define nest[0] in order to save space
+ // define abs_space[0] in order to save space
if(depth==0) {
- data.bplus_ = vec(dimen);
+ data.bplus_ = vec(dimen);
data.bplus_[1] = 1;
}
else {
- data.bplus_ = getbasis1(data.nest_);
+ //data.bplus_ = getbasis1(data.abs_space_);
+ data.bplus_ = make_basis1(data);
}
-
+
return;
}
- ssubspace* s = data.nest_; // only used when depth>0
- ssubspace *spm0, *spm;
+ ssubspace* s;
+ if( bigmats ) {
+ s = data.abs_space_; // only used when depth>0
+ }
+ ssubspace *spm_rel, *spm_abs;
SCALAR eig = denom1;
// if(depth) eig*=denom(*s);
smat subconjmat; // only used when depth>0
@@ -234,51 +287,80 @@ void form_finder::make_basis( ff_data &data ) {
// will only be a 2x2 in this case (genus 1 only!)
}
else {
- subconjmat = h->s_opmat_restricted(-1,*s,1,verbose);
+ //subconjmat = h->s_opmat_restricted(-1,*s,1,verbose);
+ subconjmat = make_nested_submat(-1,data);
}
// C++11 loop over two variables (similar to python)
// for( int b : { -1,+1 } ) { /* use b as -1 or +1 */ }
for(long signeig=+1; signeig>-2; signeig-=2) {
- SCALAR seig;
+ SCALAR seig;
seig = eig;
-
+
if(signeig<0) seig =- eig;
-
+
if(depth) {
- spm0 = new ssubspace(eigenspace(subconjmat,seig));
- spm = new ssubspace(combine(*s,*spm0));
- delete spm0;
+ spm_rel = new ssubspace(eigenspace(subconjmat,seig));
+ //spm_abs = new ssubspace(combine(*s,*spm_rel));
}
else {
- spm = new ssubspace(eigenspace(subconjmat,seig));
+ spm_rel = new ssubspace(eigenspace(subconjmat,seig));
+ //spm_abs = spm_rel;
}
-
- if(dim(*spm)!=1) {
+
+ if(dim(*spm_rel)!=1) {
cout << "error in form_finder::makebasis; ";
cout << "\nfinal (";
-
- if(signeig>0) cout << "+";
+
+ if(signeig>0) cout << "+";
else cout << "-";
-
- cout << ") subspace has dimension " << dim(*spm) << endl;
+
+ cout << ") subspace has dimension " << dim(*spm_rel) << endl;
cout << "aborting this branch!" << endl;
- delete spm;
+ //delete spm_abs;
+ delete spm_rel;
return;
}
-
- if(signeig>0) data.bplus_ = getbasis1(spm);
- else data.bminus_ = getbasis1(spm);
- delete spm;
+ //vec v = getbasis1(spm_abs);
+ vec w = make_basis2(data, spm_rel->bas().as_mat().col(1));
+
+ if(signeig>0) data.bplus_ = w;
+ else data.bminus_ = w;
+
+ //delete spm_abs;
+ delete spm_rel;
}
}
-vec form_finder::getbasis1(const ssubspace* s)
+vec form_finder::make_basis2(ff_data &data, const vec& v)
+{
+ ff_data *d = &data;
+ int level = data.depth_;
+ vec w = v;
+ while (level--)
+ {
+ w = mult_mod_p(d->rel_space_->bas(), w, MODULUS);
+ d = d->parent_;
+ }
+ return lift(w);
+}
+
+vec form_finder::make_basis1(ff_data &data)
+{
+ vec v(1); v.set(1,1);
+ return make_basis2(data, v);
+}
+
+vec getbasis1(const ssubspace* s)
+{
+ return lift(basis(*s).as_mat().col(1));
+}
+
+vec lift(const vec& v)
{
- VEC b = basis(*s).as_mat().col(1);
#ifdef MODULAR
- VEC bb;
+ VEC b=v, bb;
if(lift(b,MODULUS,bb))
b = bb;
else
@@ -298,7 +380,7 @@ void form_finder::recover(vector< vector<long> > eigs) {
for(unsigned int iform=0; iform<eigs.size(); iform++) {
if(verbose) {
cout << "Form number " << iform+1 << " with eigs ";
-
+
int n = eigs[iform].size();
if(n>10) n = 10;
@@ -525,7 +607,7 @@ void form_finder::store(vec bp, vec bm, vector<long> eigs) {
gnfcount++;
// Inform about newform count
- ECLOG(0) << "Current newform subtotal count at " << gnfcount << endl;
+ ECLOG(1) << "Current newform subtotal count at " << gnfcount << endl;
}
#if (METHOD==2)
diff --git a/libsrc/xsplit_data.cc b/libsrc/xsplit_data.cc
index c0cffe2..d534c42 100644
--- a/libsrc/xsplit_data.cc
+++ b/libsrc/xsplit_data.cc
@@ -39,7 +39,8 @@ ff_data::ff_data( form_finder* ff )
subdim_( 0 ),
eigenvalue_( 0 ),
eiglist_(),
- nest_( NULL ),
+ abs_space_( NULL ),
+ rel_space_( NULL ),
conjmat_(),
the_opmat_(),
parent_( NULL ),
@@ -55,8 +56,10 @@ ff_data::ff_data( form_finder* ff )
*/
ff_data::~ff_data() {
// Delete dynamically created objects
- delete nest_;
- nest_ = NULL;
+ delete abs_space_;
+ delete rel_space_;
+ abs_space_ = NULL;
+ rel_space_ = NULL;
}
#ifdef ECLIB_MULTITHREAD
@@ -98,12 +101,21 @@ nodestatus ff_data::status() {
}
/**
- * submats()
+ * abs_space()
*
- * Returns relevant subspace of current depth.
+ * Returns absolute subspace of current depth.
*/
-ssubspace* ff_data::submats() {
- return nest_;
+ssubspace* ff_data::abs_space() {
+ return abs_space_;
+}
+
+/**
+ * rel_space()
+ *
+ * Returns relative subspace of current depth.
+ */
+ssubspace* ff_data::rel_space() {
+ return rel_space_;
}
/**
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/eclib.git
More information about the debian-science-commits
mailing list