← Back to branch summary

~ubuntu-branches/ubuntu/saucy/eclib/saucy

« back to all changes in this revision

Viewing changes to libsrc/mat.cc

  • Committer: Package Import Robot
  • Author(s): Julien Puydt
  • Date: 2013-02-13 17:19:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130213171942-bxhm258djn8ev4wq
Tags: 2013-01-01-1
Upgraded to latest upstream.

Show diffs side-by-side

added added

removed removed

Lines of Context:
libsrc/mat.cc
1547
1547
 return m.slice(rk,nc);
1548
1548
}
1549
1549
 
1550
 
// The following function computes the upper-triangular echelon form of m modulo the prime pr.
 
1550
// Possible FLINT_LEVEL values are as follows.
 
1551
//
 
1552
// 0: no FLINT support (or a version <2.3)
 
1553
// 1: support for 64-bit nmod_mat (standard from version 2.3)
 
1554
// 2: support for 32-bit hmod_mat (non-standard, from version 2.3)
 
1555
//
 
1556
// The configure script should have detected whether the functions
 
1557
// nmod_mat_rref and/or hmod_mat_rref are available to be used here.
 
1558
//
 
1559
#if FLINT_LEVEL!=0
 
1560
//#define TRACE_FLINT_RREF
 
1561
#include "flint/fmpz.h"
 
1562
#if (SCALAR_OPTION==1)&&(FLINT_LEVEL==2)
 
1563
#include "flint/hmod_mat.h"
 
1564
#undef uscalar
 
1565
#undef mod_mat
 
1566
#undef mod_mat_init
 
1567
#undef mod_mat_clear
 
1568
#undef mod_mat_entry
 
1569
#undef mod_mat_rref
 
1570
#define uscalar hlimb_t // unsigned int
 
1571
#define mod_mat hmod_mat_t // uses unsigned ints
 
1572
#define mod_mat_init hmod_mat_init
 
1573
#define mod_mat_clear hmod_mat_clear
 
1574
#define mod_mat_entry hmod_mat_entry
 
1575
#define mod_mat_rref hmod_mat_rref
 
1576
#else
 
1577
#include "flint/nmod_mat.h"
 
1578
#undef uscalar
 
1579
#undef mod_mat
 
1580
#undef mod_mat_init
 
1581
#undef mod_mat_clear
 
1582
#undef mod_mat_entry
 
1583
#undef mod_mat_rref
 
1584
#define uscalar mp_limb_t // unsigned long
 
1585
#define mod_mat nmod_mat_t // uses unsigned longs
 
1586
#define mod_mat_init nmod_mat_init
 
1587
#define mod_mat_clear nmod_mat_clear
 
1588
#define mod_mat_entry nmod_mat_entry
 
1589
#define mod_mat_rref nmod_mat_rref
 
1590
#endif
 
1591
#include "flint/profiler.h"
 
1592
 
 
1593
// FLINT has two types for modular matrices: standard in FLINT-2.3 has
 
1594
// nmod_mat_t with entries of type mp_limb_t (unsigned long);
 
1595
// non-standard (in an optional branch) is hmod_mat_t, with entries
 
1596
// hlimb_t (unsigned int).  We use the former when scalar=long and the
 
1597
// latter when scalar=int, provided that the optional functions are
 
1598
// present, which should have been determined by the configure script.
 
1599
// The unsigned scalar types are #define'd as uscalar.
 
1600
 
 
1601
mat ref_via_flint(const mat& M, scalar pr)
 
1602
{
 
1603
  long nr=nrows(M), nc=ncols(M);
 
1604
  long i, j, n=max(nr,nc);
 
1605
 
 
1606
  // copy of the modulus for FLINT
 
1607
  uscalar mod = (uscalar)pr;
 
1608
 
 
1609
  // create flint matrix copy of M:
 
1610
  mod_mat A;
 
1611
  mod_mat_init(A, nr, nc, mod);
 
1612
  for(i=0; i<nr; i++)
 
1613
    for(j=0; j<nc; j++)
 
1614
      mod_mat_entry(A,i,j) = M(i+1,j+1);
 
1615
 
 
1616
  // reduce A to rref:
 
1617
#ifdef TRACE_FLINT_RREF
 
1618
  timeit_t t;
 
1619
  timeit_start(t);
 
1620
  cerr<<"(nr,nc)=("<<nr<<","<<nc<<"): "<<flush;
 
1621
#endif
 
1622
  long rank = mod_mat_rref(A);
 
1623
#ifdef TRACE_FLINT_RREF
 
1624
  timeit_stop(t);
 
1625
  cerr<<" cpu = "<<(t->cpu)<<" ms, wall = "<<(t->wall)<<" ms"<<endl;
 
1626
#endif
 
1627
 
 
1628
  // copy back to a new matrix for return:
 
1629
  mat ans(nr,nc);
 
1630
  for(i=0; i<nr; i++)
 
1631
    for(j=0; j<nc; j++)
 
1632
      ans(i+1,j+1) = mod_mat_entry(A,i,j);
 
1633
 
 
1634
  mod_mat_clear(A);
 
1635
  return ans;
 
1636
}
 
1637
 
 
1638
// The following function computes the reduced echelon form
 
1639
// of M modulo the prime pr, calling FLINT's nmod_mat_rref function.
 
1640
 
 
1641
mat ref_via_flint(const mat& M, vec& pcols, vec& npcols,
 
1642
                                  long& rk, long& ny, scalar pr)
 
1643
{
 
1644
  long nr=nrows(M), nc=ncols(M);
 
1645
  long i, j, k, n=max(nr,nc);
 
1646
 
 
1647
#ifdef TRACE_FLINT_RREF
 
1648
#if (SCALAR_OPTION==1)
 
1649
  cout << "In ref_via_flint(M) with M having "<<nr<<" rows and "<<nc<<" columns, using hmod_mat and modulus "<<pr<<"."<<endl;
 
1650
#else
 
1651
  cout << "In ref_via_flint(M) with M having "<<nr<<" rows and "<<nc<<" columns, using nmod_mat and modulus "<<pr<<"."<<endl;
 
1652
#endif
 
1653
  //  cout << "Size of  scalar = "<<8*sizeof(scalar)<<" bits"<<endl;
 
1654
  //  cout << "Size of uscalar = "<<8*sizeof(uscalar)<<" bits"<<endl;
 
1655
#endif
 
1656
 
 
1657
  // copy of the modulus
 
1658
  uscalar mod = (uscalar)pr;
 
1659
 
 
1660
  // create flint matrix copy of M:
 
1661
  mod_mat A;
 
1662
  mod_mat_init(A, nr, nc, mod);
 
1663
  for(i=0; i<nr; i++)
 
1664
    for(j=0; j<nc; j++)
 
1665
      mod_mat_entry(A,i,j) = (uscalar)posmod(M(i+1,j+1),pr);
 
1666
 
 
1667
#ifdef TRACE_FLINT_RREF
 
1668
  timeit_t t;
 
1669
  timeit_start(t);
 
1670
  cerr<<"(nr,nc)=("<<nr<<","<<nc<<"): "<<flush;
 
1671
#endif
 
1672
 
 
1673
  // reduce A to rref:
 
1674
  rk = mod_mat_rref(A);
 
1675
#ifdef TRACE_FLINT_RREF
 
1676
  timeit_stop(t);
 
1677
  cerr<<"rank = "<<rk<<". cpu = "<<(t->cpu)<<" ms, wall = "<<(t->wall)<<" ms"<<endl;
 
1678
#endif
 
1679
 
 
1680
  // construct vectors of pivotal and non-pivotal columns
 
1681
  ny = nc-rk;
 
1682
  pcols.init(rk);
 
1683
  npcols.init(ny);
 
1684
  for (i = j = k = 0; i < rk; i++)
 
1685
    {
 
1686
      while (mod_mat_entry(A, i, j) == 0UL)
 
1687
        {
 
1688
          npcols[k+1] = j+1;
 
1689
          k++;
 
1690
          j++;
 
1691
        }
 
1692
      pcols[i+1] = j+1;
 
1693
      j++;
 
1694
    }
 
1695
  while (k < ny)
 
1696
    {
 
1697
      npcols[k+1] = j+1;
 
1698
      k++;
 
1699
      j++;
 
1700
    }
 
1701
 
 
1702
  // copy back to a new matrix for return:
 
1703
  mat ans(rk,nc);
 
1704
  for(i=0; i<rk; i++)
 
1705
    for(j=0; j<nc; j++)
 
1706
      ans(i+1,j+1) = mod_mat_entry(A,i,j);
 
1707
 
 
1708
  mod_mat_clear(A);
 
1709
  return ans;
 
1710
}
 
1711
#endif // FLINT_LEVEL
 
1712
 
 
1713
// The following function computes the upper-triangular echelon form
 
1714
// of m modulo the prime pr.
1551
1715
 
1552
1716
mat echmodp_uptri(const mat& entries, vec& pcols, vec& npcols,
1553
1717
                                  long& rk, long& ny, scalar pr)
1554
1718
{
1555
 
// cout << "In echmodp_uptri with entries = " << entries;
 
1719
// cout << "In echmodp_uptri with matrix = " << entries;
1556
1720
 long nc,nr,r,c,r2,r3,rmin;
1557
1721
 scalar min, mr2c;
1558
1722
 nr=entries.nro, nc=entries.nco;
1568
1732
     mij=m.entries+(r-1)*nc+c-1;
1569
1733
     min = *mij;   rmin = r;
1570
1734
     for (r2=r+1, mij+=nc; (r2<=nr)&&(min==0); r2++, mij+=nc)
1571
 
       { 
 
1735
       {
1572
1736
         mr2c = *mij;
1573
1737
         if (0!=mr2c) { min=mr2c; rmin=r2 ;}
1574
1738
       }
1575
1739
     if (min==0) npcols[++ny] = c;
1576
 
     else       
 
1740
     else
1577
1741
       {
1578
1742
         pcols[++rk] = c;
1579
1743
         if (rmin>r) m.swaprows(r,rmin);