root / tests / test-krylov-elim.C

/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
//--------------------------------------------------------------------------
//          Test for the krylov-elimination
//--------------------------------------------------------------------------
// usage: test-krylov-elim p A, to compute the rank profile of the (n+m)xn matrix B
// formed by the n identity vectors and the mxn matrix A over Z/pZ
//-------------------------------------------------------------------------

//-------------------------------------------------------------------------
#define DEBUG 0

#include <iostream>
#include "Matio.h"
#include "timer.h"
using namespace std;
#include "fflas-ffpack/modular-balanced.h"
#include "fflas-ffpack/ffpack.h"


typedef Modular<double> Field;

template<class T>
std::ostream& printvect(std::ostream& o, T* vect, size_t dim){
        for(size_t i=0; i<dim; ++i)
                o << vect[i] << " " ;
        return o << std::endl;
}

int main(int argc, char** argv){

        int m,n;

        
        if (argc!=3){
                cerr<<"usage : test-lqup <p> <A>"<<endl
                    <<"         to compute the rank profile of the (n+m)xn matrix B formed by the n identity vectors and the mxn matrix A over Z/pZ"
                    <<endl;
                exit(-1);
        }
        Field F(atoi(argv[1]));
        Field::Element one;
        F.init(one, 1UL);
        Field::Element * A = read_field<Field> (F,argv[2],&m,&n);

        Field::Element * B = new Field::Element[(m+n)*n];
        for (size_t i=0; i<(n+m)*n;++i) *(B+i)=0;

        size_t deg = (n-1)/m+1;
        size_t curr_row = 0;
        size_t it_idx = 0;
        size_t bk_idx = 0;
        for (size_t i=0; i<m; ++i){
                for (size_t j=0; j<deg; ++j){
                        if (curr_row < n+m -1){
                                F.assign( *(B + curr_row*n + n-1 - it_idx), one);
                                curr_row++;
                                it_idx++;
                        }
                }
                for (size_t j=0; j<n; ++j)
                        *(B + curr_row*n + j) = *(A + bk_idx*n + j);
                bk_idx++;
                curr_row++;
        }
        write_field (F, cout<<"A = "<<endl, A, m,n,n);
        write_field (F, cout<<"B = "<<endl, B, m+n,n,n);

        size_t *rp = new size_t[n];

        FFPACK::SpecRankProfile(F, m, n, A, n, deg,rp);

        size_t * P = new size_t[n];
        size_t * Q = new size_t[n+m];
        FFPACK::LUdivine(F, FFLAS::FflasNonUnit, m+n, n, B, n, P, Q, FFPACK::FfpackLQUP);

        printvect (cout<<"RankProfile (A) = "<<endl, rp, n)<<endl;

        printvect (cout<<"RankProfile (B) = "<<endl, Q, n)<<endl;


        
        
        return 0;              
}