examples/omp_smithvalence.C

Valence of sparse matrix over Z or Zp.

/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s

/*
 * examples/omp_smithvalence.C
 *
 * Copyright (C) 2010 J-G Dumas
 *
 * This file is part of LinBox.
 *
 *   LinBox is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU Lesser General Public License as
 *   published by the Free Software Foundation, either version 2 of
 *   the License, or (at your option) any later version.
 *
 *   LinBox is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU Lesser General Public License for more details.
 *
 *   You should have received a copy of the GNU Lesser General Public
 *   License along with LinBox.  If not, see
 *   <http://www.gnu.org/licenses/>.
 */

//#undef _OPENMP
//#define _OPENMP
#include <omp.h>
//#include "linbox-config.h"
#include <iostream>

#include <givaro/givintnumtheo.h>
using namespace Givaro;
#include "linbox/field/gf2.h"
#include "linbox/field/modular.h"
#include "linbox/field/givaro.h"
#include "linbox/field/field-traits.h"
#include "linbox/blackbox/transpose.h"
#include "linbox/blackbox/compose.h"
#include "linbox/blackbox/sparse.h"
#include "linbox/solutions/rank.h"
#include "linbox/solutions/valence.h"
#include "linbox/algorithms/smith-form-sparseelim-local.h"
#include "linbox/util/matrix-stream.h"
#include "linbox/util/timer.h"


template<class Field>
unsigned long& TempLRank(unsigned long& r, char * filename, const Field& F)
{
    std::ifstream input(filename);
    LinBox::MatrixStream< Field > msf( F, input );
    LinBox::SparseMatrix<Field,typename LinBox::Vector<Field>::SparseSeq> FA(msf);
    input.close();
    Timer tim; tim.start();
    LinBox::rankin(r, FA);
    tim.stop();
    F.write(std::cout << "Rank over ") << " is " << r << ' ' << tim << std::endl;
    return r;
}

unsigned long& TempLRank(unsigned long& r, char * filename, const LinBox::GF2& F2)
{
    std::ifstream input(filename);
    LinBox::ZeroOne<LinBox::GF2> A;
    A.read(input);
    input.close();
    LinBox::Timer tim; tim.start();
    LinBox::rankin(r, A, LinBox::Method::SparseElimination() );
    tim.stop();
    F2.write(std::cout << "Rank over ") << " is " << r << ' ' << tim << std::endl;
    return r;
}

unsigned long& LRank(unsigned long& r, char * filename,Givaro::Integer p)
{

Givaro::Integer maxmod16; LinBox::FieldTraits<LinBox::GivaroZpz<Std16> >::maxModulus(maxmod16);
Givaro::Integer maxmod32; LinBox::FieldTraits<LinBox::GivaroZpz<Std32> >::maxModulus(maxmod32);
Givaro::Integer maxmod53; LinBox::FieldTraits<LinBox::Modular<double> >::maxModulus(maxmod53);
Givaro::Integer maxmod64; LinBox::FieldTraits<LinBox::GivaroZpz<Std64> >::maxModulus(maxmod64);
    if (p == 2) {
        LinBox::GF2 F2;
        return TempLRank(r, filename, F2);
    }
    else if (p <= maxmod16) {
        typedef LinBox::GivaroZpz<Std16> Field;
        Field F(p);
        return TempLRank(r, filename, F);
    }
    else if (p <= maxmod32) {
        typedef LinBox::GivaroZpz<Std32> Field;
        Field F(p);
        return TempLRank(r, filename, F);
    }
    else if (p <= maxmod53) {
        typedef LinBox::Modular<double> Field;
        Field F(p);
        return TempLRank(r, filename, F);
    }
    else if (p <= maxmod64) {
        typedef LinBox::GivaroZpz<Std64> Field;
        Field F(p);
        return TempLRank(r, filename, F);
    }
    else {
        typedef LinBox::GivaroZpz<Givaro::Integer> Field;
        Field F(p);
        return TempLRank(r, filename, F);
    }
    return r;
}

std::vector<size_t>& PRank(std::vector<size_t>& ranks, char * filename,Givaro::Integer p, size_t e, size_t intr)
{
Givaro::Integer maxmod;
    LinBox::FieldTraits<LinBox::GivaroZpz<Std64> >::maxModulus(maxmod);
    if (p <= maxmod) {
        typedef LinBox::GivaroZpz<Std64> Ring;
        int64_t lp(p);
    Givaro::Integer q = pow(p,e); int64_t lq(q);
        if (q >Givaro::Integer(lq)) {
            std::cerr << "Power rank might need extra large composite (" << p << '^' << e << ")." << std::endl << "Such a large composite is not yet implemented ..." << std::endl;
            q = p;
            do {
                q *= p; lq = (int64_t)q;
            } while (q ==Givaro::Integer(lq));
            q/=p; lq = (int64_t)q;
            std::cerr << "First trying: " << lq << " (without further warning this will be sufficient)." << std::endl;
        }
        Ring F(lq);
        std::ifstream input(filename);
        LinBox::MatrixStream<Ring> ms( F, input );
        LinBox::SparseMatrix<Ring, LinBox::Vector<Ring>::SparseSeq > A (ms);
        input.close();
        LinBox::PowerGaussDomain< Ring > PGD( F );

        PGD.prime_power_rankin( lq, lp, ranks, A, A.rowdim(), A.coldim(), std::vector<size_t>());
        F.write(std::cout << "Ranks over ") << " are " ;
        for(std::vector<size_t>::const_iterator rit=ranks.begin(); rit != ranks.end(); ++rit)
            std::cout << *rit << ' ';
        std::cout << std::endl;
    }
    else {
        std::cerr << "*** WARNING *** Sorry power rank mod large composite not yet implemented" << std::endl;
        std::cerr << "*** WARNING *** Assuming integer rank, extra factors in the Smith form could be missing" << std::endl;
        ranks.resize(0); ranks.push_back(intr);
    }
    return ranks;
}

using namespace LinBox;

int main (int argc, char **argv)
{
std::cout << "num process: " << omp_get_num_procs() << std::endl;
std::cout << "num threads: " << omp_get_num_threads() << std::endl;
std::cout << "max threads: " << omp_get_max_threads() << std::endl;
omp_set_dynamic(10);
std::cout << "dyn threads: " << omp_get_dynamic() << std::endl;
    //     commentator.setMaxDetailLevel (-1);
    //     commentator.setMaxDepth (-1);
    //     commentator.setReportStream (std::cerr);


    if (argc < 2 || argc > 4) {
        std::cerr << "Usage: valence <matrix-file-in-supported-format> [-ata|-aat|valence] [coprime]" << std::endl;
        return -1;
    }

    std::ifstream input (argv[1]);
    if (!input) { std::cerr << "Error opening matrix file " << argv[1] << std::endl; return -1; }

    PID_integer ZZ;
    MatrixStream< PID_integer > ms( ZZ, input );
    typedef SparseMatrix<PID_integer>  Blackbox;
    Blackbox A (ms);
    input.close();

    std::cout << "A is " << A.rowdim() << " by " << A.coldim() << std::endl;

    PID_integer::Element val_A;

    LinBox::Timer chrono; chrono.start();
    if (argc >= 3) {
        Transpose<Blackbox> T(&A);
        if (strcmp(argv[2],"-ata") == 0) {
            Compose< Transpose<Blackbox>, Blackbox > C (&T, &A);
            std::cout << "A^T A is " << C.rowdim() << " by " << C.coldim() << std::endl;
            valence(val_A, C);
        }
        else if (strcmp(argv[2],"-aat") == 0) {
            Compose< Blackbox, Transpose<Blackbox> > C (&A, &T);
            std::cout << "A A^T is " << C.rowdim() << " by " << C.coldim() << std::endl;
            valence(val_A, C);
        }
        else {
            std::cout << "Suppose primes are contained in " << argv[2] << std::endl;
            val_A = LinBox::Integer(argv[2]);
        }
    }
    else {
        if (A.rowdim() != A.coldim()) {
            std::cerr << "Valence works only on square matrices, try either to change the dimension in the matrix file, or to compute the valence of A A^T or A^T A, via the -aat or -ata options."  << std::endl;
            exit(0);
        }
        else
            valence (val_A, A);
    }

    std::cout << "Valence is " << val_A << std::endl;

    std::vector<Givaro::Integer> Moduli;
    std::vector<size_t> exponents;
    IntFactorDom<> FTD;

    typedef std::pair<Givaro::Integer,unsigned long> PairIntRk;
    std::vector< PairIntRk > smith;


Givaro::Integer coprimeV=2;
    if (argc >= 4) {
        coprimeV =Givaro::Integer(argv[3]);
    }
    while ( gcd(val_A,coprimeV) > 1 ) {
        FTD.nextprimein(coprimeV);
    }

    if (argc >= 4) {
        std::cout << "Suppose " << argv[3] << " is coprime with Smith form" << std::endl;
    }

    std::cout << "Integer rank: " << std::endl;

    unsigned long coprimeR; LRank(coprimeR, argv[1], coprimeV);
    smith.push_back(PairIntRk(coprimeV, coprimeR));
    //         std::cerr << "Rank mod " << coprimeV << " is " << coprimeR << std::endl;

    std::cout << "Some factors (50000 factoring loop bound): ";
    FTD.set(Moduli, exponents, val_A, 50000);
    std::vector<size_t>::const_iterator eit=exponents.begin();
    for(std::vector<Givaro::Integer>::const_iterator mit=Moduli.begin();
        mit != Moduli.end(); ++mit,++eit)
        std::cout << *mit << '^' << *eit << ' ';
    std::cout << std::endl;

    std::vector<Givaro::Integer> SmithDiagonal(coprimeR,Givaro::Integer(1));

        omp_set_num_threads(omp_get_max_threads());    
    std::cout << "num procs: " << omp_get_num_procs() << std::endl;
    std::cout << "max threads: " << omp_get_max_threads() << std::endl;
    std::cout << "num threads: " << omp_get_num_threads() << std::endl;
#pragma omp parallel for shared(SmithDiagonal, Moduli, coprimeR)
        
    for(size_t j=0; j<Moduli.size(); ++j) {
        omp_set_num_threads(omp_get_max_threads());    
    std::cout << "for num procs: " << omp_get_num_procs() << std::endl;
    std::cout << "for max threads: " << omp_get_max_threads() << std::endl;
    std::cout << "for num threads: " << omp_get_num_threads() << std::endl;
        unsigned long r; LRank(r, argv[1], Moduli[j]);
        std::cerr << "Rank mod " << Moduli[j] << " is " << r << " on thread: " << omp_get_thread_num() << std::endl;
        smith.push_back(PairIntRk( Moduli[j], r));
        for(size_t i=r; i < coprimeR; ++i)
            SmithDiagonal[i] *= Moduli[j];
    }


    /*
       for(std::vector<Givaro::Integer>::const_iterator mit=Moduli.begin();
       mit != Moduli.end(); ++mit) {
       unsigned long r; LRank(r, argv[1], *mit);
       std::cerr << "Rank mod " << *mit << " is " << r << std::endl;
       smith.push_back(PairIntRk(*mit, r));
       for(size_t i=r; i < coprimeR; ++i)
       SmithDiagonal[i] *= *mit;
       }
       */

    eit=exponents.begin();
    std::vector<PairIntRk>::const_iterator sit=smith.begin();
    for( ++sit; sit != smith.end(); ++sit, ++eit) {
        if (sit->second != coprimeR) {
            std::vector<size_t> ranks;
            ranks.push_back(sit->second);
            if (*eit > 1) {
                PRank(ranks, argv[1], sit->first, *eit, coprimeR);
            }
            else {
                PRank(ranks, argv[1], sit->first, 2, coprimeR);
            }
            if (ranks.size() == 1) ranks.push_back(coprimeR);
            for(size_t expo = (*eit)<<1; ranks.back() < coprimeR; expo<<=1) {
                PRank(ranks, argv[1], sit->first, expo, coprimeR);
                if (ranks.size() < expo) {
                    std::cerr << "*** WARNING *** Larger prime power not yet implemented" << std::endl;
                    break;
                }
            }
            std::vector<size_t>::const_iterator rit=ranks.begin();
            unsigned long modrank = *rit;
            for(++rit; rit!= ranks.end(); ++rit) {
                if ((*rit)>= coprimeR) break;
                for(size_t i=(*rit); i < coprimeR; ++i)
                    SmithDiagonal[i] *= sit->first;
                modrank = *rit;
            }
        }
    }

Givaro::Integer si=1;
    size_t num=0;
    for( std::vector<Givaro::Integer>::const_iterator dit=SmithDiagonal.begin();
         dit != SmithDiagonal.end(); ++dit) {
        if (*dit == si) ++num;
        else {
            std::cerr << '[' << si << ',' << num << "] ";
            num=1;
            si = *dit;
        }
    }
    std::cerr << '[' << si << ',' << num << "] " << std::endl;
    chrono.stop();
    std::cerr << chrono << std::endl;


    return 0;
}