diff --git a/.gitignore b/.gitignore index a7b4f2be..89df69dd 100644 --- a/.gitignore +++ b/.gitignore @@ -109,6 +109,7 @@ tests/test-charpoly tests/test-maxdelayeddim tests/test-quasisep tests/test-sss +tests/test-rrr test-det-check test-fdot test-fspmm-recint diff --git a/benchmarks/Makefile.am b/benchmarks/Makefile.am index 43324f64..727ed9ea 100755 --- a/benchmarks/Makefile.am +++ b/benchmarks/Makefile.am @@ -33,7 +33,7 @@ endif PERFPUBLISHERFILE=benchmarks-report.xml -FFLA_BENCH = benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm benchmark-fgesv benchmark-ftrsv benchmark-ftrtri benchmark-inverse benchmark-fsytrf benchmark-fsyrk benchmark-lqup benchmark-fsyr2k benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-quasisep benchmark-sss benchmark-storage-transpose benchmark-qscomp +FFLA_BENCH = benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm benchmark-fgesv benchmark-ftrsv benchmark-ftrtri benchmark-inverse benchmark-fsytrf benchmark-fsyrk benchmark-lqup benchmark-fsyr2k benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-quasisep benchmark-sss benchmark-storage-transpose benchmark-qscomp benchmark-rrrgen benchmark-rrroperations benchmark-SSSvsRRR BLAS_BENCH = benchmark-sgemm$(EXEEXT) benchmark-dgemm benchmark-dtrsm LAPA_BENCH = benchmark-dtrtri benchmark-dgetri benchmark-dgetrf benchmark-dsytrf @@ -77,6 +77,9 @@ benchmark_fsyrk_SOURCES = benchmark-fsyrk.C benchmark_quasisep_SOURCES = benchmark-quasisep.C benchmark_qscomp_SOURCES = benchmark-qscomp.C benchmark_sss_SOURCES = benchmark-sss.C +benchmark_SSSvsRRR_SOURCES = benchmark-SSSvsRRR.C +benchmark_rrrgen_SOURCES = benchmark-rrrgen.C +benchmark_rrroperations_SOURCES = benchmark-rrroperations.C benchmark_charpoly_SOURCES = benchmark-charpoly.C benchmark_charpoly_mp_SOURCES = benchmark-charpoly-mp.C benchmark_lqup_SOURCES = benchmark-lqup.C diff --git a/benchmarks/benchmark-SSSvsRRR.C b/benchmarks/benchmark-SSSvsRRR.C new file mode 100644 index 00000000..64a07d7a --- /dev/null +++ b/benchmarks/benchmark-SSSvsRRR.C @@ -0,0 +1,271 @@ +/* Copyright (c) FFLAS-FFPACK + * Written by Hippolyte Signargout + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK 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.1 of the License, or (at your option) any later version. + * + * This library 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 this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + */ + +// Template from benchmark-quasisep.C + +#define __FFLASFFPACK_OPENBLAS_NT_ALREADY_SET 1 + +#include "fflas-ffpack/fflas-ffpack-config.h" +#include +#include +#include +#include "fflas-ffpack/fflas-ffpack.h" +#include "fflas-ffpack/utils/timer.h" +#include "fflas-ffpack/utils/test-utils.h" +#include "fflas-ffpack/utils/fflas_randommatrix.h" +#include "fflas-ffpack/utils/fflas_io.h" +#include "fflas-ffpack/utils/args-parser.h" +#include "fflas-ffpack/ffpack/ffpack_rrrgen.inl" + + + +using namespace std; +using namespace FFLAS; +using namespace FFPACK; + + + + +template +void run_with_field(int q, size_t n, size_t m, size_t s, size_t r, size_t iter, uint64_t seed){ + Field F(q); + typedef typename Field::Element_ptr Element_ptr; + FFLAS::Timer chrono; + size_t lda=n; + size_t ldts = m; + size_t rs = n%s; // Size of the partial block + size_t ls = (rs)? rs: s; // Size of the last block + + double time_gens = 0, time_sssxts =0; + Element_ptr A = FFLAS::fflas_new (F, n, n); + Element_ptr A2 = FFLAS::fflas_new (F, n, n); + Element_ptr B = FFLAS::fflas_new (F, n, n); + Element_ptr B2 = FFLAS::fflas_new (F, n, n); + Element_ptr TSS = FFLAS::fflas_new (F, n, m); + Element_ptr D = fflas_new (F, n, s); + Element_ptr P = fflas_new (F, n - s, s); + Element_ptr Q = fflas_new (F, n - ls, s); + Element_ptr R = fflas_new (F, ((n > (s + ls))? (n - s - ls): 0), s); + Element_ptr U = fflas_new (F, n - ls, s); + Element_ptr V = fflas_new (F, n - ls, s); + Element_ptr W = fflas_new (F, ((n > (s + ls))? (n - s - ls): 0), s); + Element_ptr Res = fflas_new(F, n, m); // Inadequate name + size_t * p = FFLAS::fflas_new (n); + for (size_t i = 0; i < ceil(n/2.); i++) + { + p[i] = n - i - 1; + } + + double time_genb = 0, time_cbxts =0; + Element_ptr H = FFLAS::fflas_new (F, n, 1); + Element_ptr TSB = FFLAS::fflas_new (F, n, m); + size_t * pa = fflas_new (n); + size_t * qa = fflas_new (n); + Element_ptr L = fflas_new(F,n,n); + Element_ptr Ua = fflas_new(F,n,n); + Element_ptr Xu = fflas_new(F, 2*s, n); + size_t * Ku = fflas_new (r+1); + size_t * Mu = fflas_new (n); + size_t * Tu = fflas_new(r); + Element_ptr Xl = fflas_new(F, n, 2*s); + size_t * Kl = fflas_new (r+1); + size_t * Ml = fflas_new (n); + size_t * Tl = fflas_new(r); + size_t * pb = fflas_new (n); + size_t * qb = fflas_new (n); + Element_ptr Lb= fflas_new(F,n,n); + Element_ptr Ub= fflas_new(F,n,n); + Element_ptr Xub= fflas_new(F, 2*s, n); + size_t * Kub= fflas_new (r+1); + size_t * Mub= fflas_new (n); + size_t * Tub= fflas_new(r); + Element_ptr Xlb= fflas_new(F, n, 2*s); + size_t * Klb= fflas_new (r+1); + size_t * Mlb= fflas_new (n); + size_t * Tlb= fflas_new(r); + size_t r2; + size_t r3; + Element_ptr CBruhat = fflas_new(F, n, m); + + double time_genr = 0, time_rrrxts =0; + RRRgen* RRRA; + Element_ptr Result = fflas_new(F, n, m); + + for (size_t i=0;i(F, n, s, A, lda,false,true); + chrono.stop(); + time_genr+=chrono.usertime(); + + // SSS generation + chrono.clear(); + chrono.start(); + DenseToSSS (F, n, s, P, s, Q, s, R, s, U, s, V, s, W, s, D, s, A, n); + chrono.stop(); + time_gens+=chrono.usertime(); + + // CB gen + chrono.clear(); + chrono.start(); + r2 = LTBruhatGen (F, FflasNonUnit, n, A2, lda, pa, qa); + r3 = LTBruhatGen (F, FflasNonUnit, n, B2, lda, pb, qb); + chrono.stop(); + if ((r2!=r)||(r3 != r)) + { + std::cerr<<"ERROR: r != r2 or r3"< >(q, n, m, t, r, iter, seed); + FFLAS::writeCommandString(std::cout, as); + std::cout<s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/benchmarks/benchmark-qscomp.C b/benchmarks/benchmark-qscomp.C index 904f70fa..2d1036e7 100644 --- a/benchmarks/benchmark-qscomp.C +++ b/benchmarks/benchmark-qscomp.C @@ -96,7 +96,7 @@ void run_with_field(int q, size_t n, size_t m, size_t s, size_t r, size_t iter, size_t * Mlb= fflas_new (n); size_t * Tlb= fflas_new(r); size_t r2; - size_t r3; + size_t r3; Element_ptr CBruhat = fflas_new(F, n, m); for (size_t i=0;i= 1000 && r <= 1400){ + size_t * rows1 = FFLAS::fflas_new (r); + size_t * cols1 = FFLAS::fflas_new (r); + + chrono.clear(); + chrono.start(); + RandomLTQSRankProfileMatrix (n, r, t, rows1, cols1); + chrono.stop(); + + time_RPMGen += chrono.usertime(); + FFLAS::fflas_delete(rows1,cols1); + } + + + //std::cout << "rows :"; + //for ( size_t i = 0; i (r); + size_t * cols2 = FFLAS::fflas_new (r); + + chrono.clear(); + chrono.start(); + RandomLTQSRankProfileMatrix_Tom (n, r, t, rows2, cols2); + chrono.stop(); + + time_RPMGen_Tom += chrono.usertime(); + + FFLAS::fflas_delete(rows2,cols2); + A = FFLAS::fflas_new (F, n, n); size_t lda=n; TS = FFLAS::fflas_new (F, n, m); @@ -104,8 +143,9 @@ void run_with_field(int q, size_t n, size_t m, size_t t, size_t r, size_t iter, } // ----------- // Standard output for benchmark - Alexis Breust 2014/11/14 - std::cout << "Time: " << (time_gen + time_cbxts) / double(iter) << " Gfops: Irrelevant (Generator) Specific times: " << time_gen / double(iter)<<" (for construction)" << time_cbxts / double(iter)<<" (for CB x TS)" ; - + // std::cout << "Time: " << (time_gen + time_cbxts) / double(iter) << " Gfops: Irrelevant (Generator) Specific times: " << time_gen / double(iter)<<" (for construction)" << time_cbxts / double(iter)<<" (for CB x TS)" ; + std::cout << "Time_stand : " << time_RPMGen / double(iter) << std::endl; + std::cout << "Time_Tom : " << time_RPMGen_Tom / double(iter) << std::endl; } int main(int argc, char** argv) { diff --git a/benchmarks/benchmark-rrrgen.C b/benchmarks/benchmark-rrrgen.C new file mode 100644 index 00000000..9399dee6 --- /dev/null +++ b/benchmarks/benchmark-rrrgen.C @@ -0,0 +1,247 @@ +/* Copyright (c) FFLAS-FFPACK + * Written by Hippolyte Signargout + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK 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.1 of the License, or (at your option) any later version. + * + * This library 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 this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + */ + +// Template from benchmark-quasisep.C + +#define __FFLASFFPACK_OPENBLAS_NT_ALREADY_SET 1 + +#include "fflas-ffpack/fflas-ffpack-config.h" +#include +#include +#include + +#include "fflas-ffpack/fflas-ffpack.h" +#include "fflas-ffpack/utils/timer.h" +#include "fflas-ffpack/utils/test-utils.h" +#include "fflas-ffpack/utils/fflas_randommatrix.h" +#include "fflas-ffpack/utils/fflas_io.h" +#include "fflas-ffpack/utils/args-parser.h" +#include "fflas-ffpack/ffpack/ffpack_rrrgen.inl" + + + +using namespace std; +using namespace FFLAS; +using namespace FFPACK; + + + + +template +void run_with_field(int q, size_t n, size_t m, size_t t, size_t r, size_t iter, uint64_t seed){ + + Field F(q); + typedef typename Field::Element_ptr Element_ptr; + typename Field::RandIter G (F, seed); + FFLAS::Timer chrono; + + size_t lda = n; + // size_t ldts = m; + + + // Element_ptr A, B, TS; + Element_ptr A; + A = FFLAS::fflas_new (F, n, n); + // B = FFLAS::fflas_new (F, n, n); + // TS = FFLAS::fflas_new (F, n, m); + // Element_ptr Res = fflas_new(F, n, m); // Inadequate name + RRRgen* RRRA; + // RRRgen* RRRB; + RRRgen* RRRres; + RRRgen* RRRL; + RRRgen* RRRU; + + + Element_ptr A2 = fflas_new (F, n, n); + // Element_ptr B2 = fflas_new (F, n, n); + double time_invert = 0, time_LU = 0; + // double time_RRRxTS = 0, time_RRRxRRR = 0,time_gen_qs = 0, time_gen_rrr = 0; + size_t * p = FFLAS::fflas_new (ceil(n/2.)); + for (size_t i = 0; i < ceil(n/2.); i++) + { + p[i] = n - i - 1; + } + + Givaro::GeneralRingNonZeroRandIter nzG (G); + for (size_t i=0; i(F, n, t, A, lda,true,true); + // chrono.stop(); + // time_gen_rrr += chrono.usertime(); + + // create RRR + // chrono.clear(); + // chrono.start(); + // RRRB = new RRRgen(F, n, t, B, lda,true,true); + // chrono.stop(); + // time_gen_rrr += chrono.usertime(); + + // RRRxTS product + // chrono.clear(); + // chrono.start(); + // RRRxTS(F,n,m,RRRA,TS,m, Res,m); + // chrono.stop(); + // time_RRRxTS += chrono.usertime(); + + // RRRxRRR + // chrono.clear(); + // chrono.start(); + // RRRres = RRRxRRR(F,RRRA,RRRB); + // chrono.stop(); + // time_RRRxRRR += chrono.usertime(); + // delete RRRres; + + + + RRRL = nullptr; + RRRU = nullptr; + chrono.clear(); + try { + chrono.start(); + LUfactRRR(F,RRRA,RRRL,RRRU); + chrono.stop(); + } + catch(...){ + delete RRRU; + delete RRRL; + delete RRRA; + FFLAS::fflas_delete(A); + FFLAS::fflas_delete(A2, p); + throw std::runtime_error("RRR Error: LU with non invertible matrix "); + } + time_LU += chrono.usertime(); + + chrono.clear(); + chrono.start(); + RRRres = RRRinvert(F,RRRA); + chrono.stop(); + time_invert += chrono.usertime(); + delete RRRres; + + + + delete RRRA; + delete RRRU; + delete RRRL; + // delete RRRB; + } + + FFLAS::fflas_delete(A); + // FFLAS::fflas_delete(B); + // FFLAS::fflas_delete(TS); + // FFLAS::fflas_delete(Res); + // FFLAS::fflas_delete(A2,B2, p); + FFLAS::fflas_delete(A2, p); + + // double mean_time_RRRxTS = time_RRRxTS / double(iter); + // double mean_time_gen_RRR = time_gen_rrr / double(2*iter); + // double mean_time_RRRxRRR = time_RRRxRRR / double(iter); + double mean_time_invert = time_invert / double(iter); + double mean_time_LU_fact = time_LU / double(iter); + double time = mean_time_invert + mean_time_LU_fact; + #define GFOPS_LU(x) (double(10*n*t*t)*log(double(n))/x) //Gfops = 10*n*s^2*log(n/s) / Time + + std::cout << "Time: " << time + << " Gfops: " << GFOPS_LU(time) + << " | details { Time : " + << "LU_RRR : " << mean_time_LU_fact << "; " + << "Invert_RRR : " << mean_time_invert << "; " + << "Gfops : " + << "Gfops_LU_RRR : " << GFOPS_LU(mean_time_LU_fact) << "; " + << "Gfops_Invert_RRR : " << GFOPS_LU(mean_time_invert) << "} "; + return ; +} + +int main(int argc, char** argv) { + +#ifdef __FFLASFFPACK_OPENBLAS_NUM_THREADS + openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); +#endif + + size_t iter = 10; + int q = 131071; + size_t n = 2167; + size_t m = 455; + size_t t = 236; + size_t r = 1100; + uint64_t seed = FFLAS::getSeed(); + + Argument as[] = { + { 'q', "-q Q", "Set the field characteristic (-1 for the ring ZZ).", TYPE_INT , &q }, + { 'n', "-n N", "Set the order of the square matrix A.", TYPE_INT , &n }, + { 'm', "-m M", "Set the column dimension of n x m RHS matrix B.", TYPE_INT , &m }, + { 't', "-t T", "Set the quasiseparability order of A.", TYPE_INT , &t }, + { 'r', "-r R", "Set the rank of each upper/lower triangular part of A.", TYPE_INT , &r }, + { 'i', "-i R", "Set number of repetitions.", TYPE_INT , &iter }, + { 's', "-s S", "Sets seed.", TYPE_INT , &seed }, + END_OF_ARGUMENTS + }; + + FFLAS::parseArguments(argc,argv,as); + run_with_field >(q, n, m, t, r, iter, seed); + FFLAS::writeCommandString(std::cout, as) << std::endl; + return 0; +} + +/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/benchmarks/benchmark-rrroperations.C b/benchmarks/benchmark-rrroperations.C new file mode 100644 index 00000000..ccd873af --- /dev/null +++ b/benchmarks/benchmark-rrroperations.C @@ -0,0 +1,205 @@ +/* Copyright (c) FFLAS-FFPACK + * Written by Hippolyte Signargout + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK 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.1 of the License, or (at your option) any later version. + * + * This library 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 this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + */ + +// Template from benchmark-quasisep.C + +#define __FFLASFFPACK_OPENBLAS_NT_ALREADY_SET 1 + +#include "fflas-ffpack/fflas-ffpack-config.h" +#include +#include + +#include "fflas-ffpack/fflas-ffpack.h" +#include "fflas-ffpack/utils/timer.h" +#include "fflas-ffpack/utils/test-utils.h" +#include "fflas-ffpack/utils/fflas_randommatrix.h" +#include "fflas-ffpack/utils/fflas_io.h" +#include "fflas-ffpack/utils/args-parser.h" +#include "fflas-ffpack/ffpack/ffpack_rrrgen.inl" + + + +using namespace std; +using namespace FFLAS; +using namespace FFPACK; + +template +void run_with_field(int q, size_t n, size_t m, size_t t, size_t r, size_t iter, uint64_t seed){ + + Field F(q); + typedef typename Field::Element_ptr Element_ptr; + typename Field::RandIter G (F, seed); + FFLAS::Timer chrono; + Element_ptr A, B; + A = FFLAS::fflas_new (F, n, n); + B = FFLAS::fflas_new (F, n, n); + + RRgen* RRA; + RRgen* RRB; + RRgen* RRres; + + RRRgen* RRRA; + RRRgen* RRRB; + RRRgen* RRRres; + + Element_ptr A2 = fflas_new (F, n, n); + Element_ptr B2 = fflas_new (F, n, n); + + size_t lda=n; + size_t ldb = n; + + double time_rrxrr = 0, time_rraddrr = 0, time_rrraddrr = 0, time_rrrxrr = 0, time_rrrxrrr = 0, time_invert = 0; + for (size_t i=0;i (ceil(n/2.)); + for (size_t i = 0; i < ceil(n/2.); i++) + { + p[i] = n - i - 1; + } + applyP (F, FFLAS::FflasLeft, FFLAS::FflasNoTrans, n, 0, ceil(n/2.), A, n, p); + applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, n, 0, ceil(n/2.), A2, n, p); + faddin (F, n, n, A2, n, A, n); + FFLAS::fflas_delete(A2, p); + + // generate B a t qsmatrix + RandomLTQSMatrixWithRankandQSorder (F,n,r,t,B, n,G); + RandomLTQSMatrixWithRankandQSorder (F,n,r,t, B2, n,G); + p = FFLAS::fflas_new (ceil(n/2.)); + for (size_t i = 0; i < ceil(n/2.); i++) + { + p[i] = n - i - 1; + } + applyP (F, FFLAS::FflasLeft, FFLAS::FflasNoTrans, n, 0, ceil(n/2.), B, n, p); + applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, n, 0, ceil(n/2.), B2, n, p); + faddin (F, n, n, B2, n, B, n); + FFLAS::fflas_delete(B2, p); + + // Computes the RRgen of A and B + RRA = new RRgen(F, n, n, A, lda); + RRB = new RRgen(F, n, n, B, ldb); + + // Computes the RRRgen of A and B + RRRA = new RRRgen(F, n, t, A, lda,true,true); + RRRB = new RRRgen(F, n, t, B, lda,true,true); + + + // RRxRR + chrono.clear(); + chrono.start(); + RRres = RRxRR(F,RRA,RRB); + chrono.stop(); + delete RRres; + time_rrxrr+=chrono.usertime(); + + // RRaddRR + chrono.clear(); + chrono.start(); + RRres = RRaddRR(F,RRA,RRB); + chrono.stop(); + delete RRres; + time_rraddrr+=chrono.usertime(); + + // RRRaddRR + chrono.clear(); + chrono.start(); + RRRres = RRRaddRR(F,RRRA,RRB); + chrono.stop(); + delete RRRres; + time_rrraddrr+=chrono.usertime(); + + // RRRxRR + chrono.clear(); + chrono.start(); + RRres = RRRxRR(F,RRRA,RRB); + chrono.stop(); + delete RRres; + time_rrrxrr+=chrono.usertime(); + + // RRRxRRR + chrono.clear(); + chrono.start(); + RRRres = RRRxRRR(F,RRRA,RRRB); + chrono.stop(); + delete RRRres; + time_rrrxrrr+=chrono.usertime(); + + // RRR invert + chrono.clear(); + chrono.start(); + RRRres = RRRinvert(F,RRRA); + chrono.stop(); + delete RRRres; + time_invert+=chrono.usertime(); + + + delete RRA; + delete RRB; + delete RRRA; + delete RRRB; + } + FFLAS::fflas_delete(A, B); + std::cout << "Total Time: " << (time_rrxrr + time_rraddrr + time_rrraddrr + time_rrrxrr + time_rrrxrrr + time_invert) / double(iter) << std::endl + << " RRxRR : " << time_rrxrr / double(iter) << std::endl + << " RR+RR : " << time_rraddrr / double(iter) << std::endl + << " RRR+RR : " << time_rrraddrr / double(iter) << std::endl + << " RRRxRR : " << time_rrrxrr / double(iter) << std::endl + << " RRRxRRR : " << time_rrrxrrr / double(iter) << std::endl + << " RRR invert : " << time_invert / double(iter) << std::endl; +} + +int main(int argc, char** argv) { + +#ifdef __FFLASFFPACK_OPENBLAS_NUM_THREADS + openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); +#endif + + size_t iter = 10; + int q = 131071; + size_t n = 2167; + size_t m = 455; + size_t t = 236; + size_t r = 1100; + uint64_t seed = FFLAS::getSeed(); + + Argument as[] = { + { 'q', "-q Q", "Set the field characteristic (-1 for the ring ZZ).", TYPE_INT , &q }, + { 'n', "-n N", "Set the order of the square matrix A.", TYPE_INT , &n }, + { 'm', "-m M", "Set the column dimension of n x m RHS matrix B.", TYPE_INT , &m }, + { 't', "-t T", "Set the quasiseparability order of A.", TYPE_INT , &t }, + { 'r', "-r R", "Set the rank of each upper/lower triangular part of A.", TYPE_INT , &r }, + { 'i', "-i R", "Set number of repetitions.", TYPE_INT , &iter }, + { 's', "-s S", "Sets seed.", TYPE_INT , &seed }, + END_OF_ARGUMENTS + }; + + FFLAS::parseArguments(argc,argv,as); + + run_with_field >(q, n, m, t, r, iter, seed); + + std::cout << "( "; + FFLAS::writeCommandString(std::cout, as) << ")" << std::endl; + return 0; +} + +/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/fflas-ffpack/ffpack/ffpack.h b/fflas-ffpack/ffpack/ffpack.h index 4f665875..77b697b0 100755 --- a/fflas-ffpack/ffpack/ffpack.h +++ b/fflas-ffpack/ffpack/ffpack.h @@ -2146,6 +2146,7 @@ namespace FFPACK { /* not used */ #include "ffpack_bruhatgen.inl" #include "ffpack_sss.inl" #include "ffpack.inl" +#include "ffpack_rrrgen.inl" #endif // __FFLASFFPACK_ffpack_H diff --git a/fflas-ffpack/ffpack/ffpack_invert.inl b/fflas-ffpack/ffpack/ffpack_invert.inl index 96058c0e..2e09a640 100644 --- a/fflas-ffpack/ffpack/ffpack_invert.inl +++ b/fflas-ffpack/ffpack/ffpack_invert.inl @@ -46,7 +46,7 @@ namespace FFPACK { } size_t * P = FFLAS::fflas_new(M); size_t * Q = FFLAS::fflas_new(M); - size_t R = ReducedRowEchelonForm (F, M, M, A, lda, P, Q, true, FfpackGaussJordanTile); + size_t R = ReducedRowEchelonForm (F, M, M, A, lda, P, Q, true, FfpackGaussJordanSlab); nullity = (int)(M - R); applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, M, 0, (int)R, A, lda, P); FFLAS::fflas_delete(P); diff --git a/fflas-ffpack/ffpack/ffpack_rrrgen.inl b/fflas-ffpack/ffpack/ffpack_rrrgen.inl new file mode 100644 index 00000000..024ca950 --- /dev/null +++ b/fflas-ffpack/ffpack/ffpack_rrrgen.inl @@ -0,0 +1,1287 @@ +#ifndef __FFLASFFPACK_ffpack_rrrgen_inl +#define __FFLASFFPACK_ffpack_rrrgen_inl + +#include + +namespace FFPACK{ + +/// @brief RRgen Class for easier representation and use. +/// A = PxLxUxQ. A (nxm) rank r +template +class RRgen { +public: + size_t n; // size of lines of original matrix + size_t m; // size of columns of original matrix + size_t r; // rank of the original matrix + typename Field::Element_ptr PL; // PL from PLUQ (n,r) + size_t ldPL; + typename Field::Element_ptr UQ; // UQ from PLUQ (r*m) + size_t ldUQ; + bool memory_owner; // whether the structure owns the memory or whether it is a view on memory area. Used in destructors + + RRgen(const Field& Fi, size_t n, size_t m, size_t r, + typename Field::Element_ptr PL, size_t ldPL, + typename Field::Element_ptr UQ,size_t ldUQ, bool memory_owner = true) + : n(n), m(m), r(r), PL(PL), ldPL(ldPL), UQ(UQ), ldUQ(ldUQ), memory_owner(memory_owner) {} + + + RRgen(const Field& Fi, size_t n_A, size_t m_A, size_t r_A, + typename Field::Element_ptr PL_A, size_t ldPL_A, + typename Field::Element_ptr UQ_A,size_t ldUQ_A, bool mem_owner, bool copy):n(n_A),m(m_A),r(r_A){ + if (copy){ + PL = FFLAS::fflas_new(Fi, n_A, r_A); + FFLAS::fassign(Fi,n_A,r_A,PL_A,ldPL_A,PL,r_A); + ldPL = r_A; + memory_owner = true; + if (UQ_A){ + UQ = FFLAS::fflas_new(Fi, r_A, m_A); + FFLAS::fassign(Fi,r_A,m_A,UQ_A,ldUQ_A,UQ,m_A); + ldUQ = m_A; + } + else { + UQ = UQ_A; + ldUQ = 0; + } + } + else { + ldPL = ldPL_A; + ldUQ = ldUQ_A; + PL = PL_A; + UQ = UQ_A; + memory_owner = mem_owner; + } + } + + + RRgen(const Field& Fi, size_t n_A, size_t m_A, typename Field::ConstElement_ptr A, size_t ldA):n(n_A),m(m_A){ + size_t* P = FFLAS::fflas_new(n_A); + size_t* Q = FFLAS::fflas_new(m_A); + + + typename Field::Element_ptr A_copy = FFLAS::fflas_new(Fi, n_A, m_A); + FFLAS::fassign(Fi, n_A, m_A, A, ldA, A_copy, m_A); + + + r = PLUQ(Fi, FFLAS::FflasNonUnit, n_A, m_A, A_copy, m_A, P, Q); + + + PL = FFLAS::fflas_new(Fi, n_A, r); + UQ = FFLAS::fflas_new(Fi, r, m_A); + ldPL = r; + ldUQ = m_A; + + // extraction of U + getTriangular(Fi, FFLAS::FflasUpper, FFLAS::FflasNonUnit, n_A, m_A, r, A_copy, m_A, UQ, m_A, true); + + // extraction of L + getTriangular(Fi, FFLAS::FflasLower, FFLAS::FflasUnit, n_A, m_A, r, A_copy, m_A, PL, r, true); + + // apply the permutations (P on L and Q on U) + applyP(Fi, FFLAS::FflasLeft, FFLAS::FflasTrans, r, 0, n_A, PL, r, P); + + + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasNoTrans, r, 0, m_A, UQ, m_A, Q); + + FFLAS::fflas_delete(A_copy); + FFLAS::fflas_delete(P,Q); + memory_owner = true; + } + + RRgen(const Field& Fi, size_t n_A, size_t m_A, typename Field::Element_ptr A, size_t ldA):n(n_A),m(m_A){ + + size_t* P = FFLAS::fflas_new(n_A); + size_t* Q = FFLAS::fflas_new(m_A); + + + + r = PLUQ(Fi, FFLAS::FflasNonUnit, n_A, m_A, A, m_A, P, Q); + + + PL = FFLAS::fflas_new(Fi, n_A, r); + UQ = FFLAS::fflas_new(Fi, r, m_A); + ldPL = r; + ldUQ = m_A; + + // extraction of U + getTriangular(Fi, FFLAS::FflasUpper, FFLAS::FflasNonUnit, n_A, m_A, r, A, m_A, UQ, m_A, true); + + // extraction of L + getTriangular(Fi, FFLAS::FflasLower, FFLAS::FflasUnit, n_A, m_A, r, A, m_A, PL, r, true); + + // apply the permutations (P on L and Q on U) + applyP(Fi, FFLAS::FflasLeft, FFLAS::FflasTrans, r, 0, n_A, PL, r, P); + + + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasNoTrans, r, 0, m_A, UQ, m_A, Q); + + FFLAS::fflas_delete(P,Q); + memory_owner = true; + } + + void RRExpand( const Field& Fi, typename Field::Element_ptr A, size_t ldA){ + + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + n, m, r, Fi.one, + PL, ldPL, + UQ, ldUQ, + Fi.zero, A, ldA); + } + + RRgen* RRcopy(const Field& Fi){ + if(UQ){ + return new RRgen(Fi,n,m,r,PL,ldPL,UQ,ldUQ,true,true); + } + else { + return new RRgen(Fi,n,m,r,PL,ldPL,nullptr,0,true,true); + } + } + + ~RRgen() { + if (memory_owner){ + + FFLAS::fflas_delete(PL); + if (UQ){ + FFLAS::fflas_delete(UQ); + } + } + } +}; + +/// @brief RRRgen Class for the tree representation. A = [A11 A12] +/// If the node is a leaf the matrix is stored in LU_right->PL and LU_left, left and right are None. [A21 A22] +template +class RRRgen { + +public: + RRgen* LU_right; // PLUQ representation of the top right submatrix (N2*rr)*(rr*N1) = A12 + RRgen* LU_left; // PLUQ representation of the down left submatrix (N1*rl)*(rl*N2) = A21 + size_t size_N1; // size of N/2 (N is the size of the matrix represented by the RRRgen) + size_t size_N2; // size of N-N1 + size_t t; // threshold for recursive representation : when N1+N2 < t, the matrix is stored completly + RRRgen* left; // recursively pointing on the same representation of the up left submatrix = A11 + RRRgen* right; // recursively pointing on the same representation of the down right submatrix = A22 + bool memory_owner; // whether the structure owns the memory or whether it is a view on memory area. Used in destructors + + + RRRgen(RRgen* LU_right, RRgen* LU_left, size_t size_N1, size_t size_N2,size_t t, RRRgen* left, RRRgen* right,bool memory_owner) + : LU_right(LU_right), LU_left(LU_left), size_N1(size_N1), size_N2(size_N2),t(t), left(left), right(right),memory_owner(memory_owner){} + + + RRRgen( const Field& F,typename Field::Element_ptr leaf, size_t n,size_t threshold,bool mem_owner,bool copy){ + + LU_right = new RRgen(F,n,n,n,leaf,n,nullptr,0,mem_owner,copy); + LU_left = nullptr; + size_N1 = n; + size_N2 = 0; + t = threshold; + left = nullptr; + right = nullptr; + memory_owner = true; + } + + RRRgen( const Field& Fi,const size_t N, const size_t threshold,typename Field::Element_ptr A, const size_t lda, bool mem_owner, bool copy = false): t (threshold) + { + if (N <= threshold) { + // leaf + LU_right = new RRgen(Fi,N,N,N,A,lda,nullptr,0,mem_owner,copy); + LU_left = nullptr; + size_N1 = N; + size_N2 = 0; + left = nullptr; + right = nullptr; + memory_owner = true; + return; + } + + size_N1 = N/2; + size_N2 = N - size_N1; + // cut the matrix A in four parts + // in case of odd N, A12 and A21 are rectangles and not squares + typename Field::Element_ptr A11 = A; + typename Field::ConstElement_ptr A12 = A + size_N1; + typename Field::ConstElement_ptr A21 = A + lda*size_N1; + typename Field::Element_ptr A22 = (typename Field::Element_ptr)A21 + size_N1; + + ///////// PLUQ Factorisation for A12 and A21 + LU_right = new RRgen(Fi, size_N1, size_N2, A12, lda); + LU_left = new RRgen(Fi, size_N2, size_N1, A21, lda); + + + // recursion on A11 and A22 + left = new RRRgen( Fi, size_N1, threshold, A11, lda, mem_owner,copy); + + right = new RRRgen( Fi, size_N2, threshold, A22, lda, mem_owner,copy); + memory_owner = true; + return; + } + + + RRRgen* RRRcopy(const Field& Fi){ + RRRgen* new_left; + RRRgen* new_right; + RRgen* new_LU_left; + RRgen* new_LU_right; + + if (left){ + new_left = left->RRRcopy(Fi); + } + else { + new_left = nullptr; + } + + if (right){ + new_right = right->RRRcopy(Fi); + } + else { + new_right = nullptr; + } + + if (LU_left){ + new_LU_left = LU_left.RRcopy(Fi); + } + else { + new_LU_left = nullptr; + } + + if (LU_right){ + new_LU_right = LU_right.RRcopy(Fi); + } + else { + new_LU_right = nullptr; + } + + return new RRRgen(new_LU_right,LU_left,size_N1,size_N2,t,new_left,new_right, true); + } + + ~RRRgen() { + if (memory_owner){ + delete LU_right; + + if (LU_left) { + delete (LU_left); + } + + if (left) { + delete left; + } + + if (right) { + delete right; + } + } + } +}; + + +/// @brief Computes the dense matrix of RRR(A) in B. B needs to be at least (A->sizeN1 + A->sizeN2) x (A->sizeN1 + A->sizeN2) preallocated. +/// @tparam Field +/// @param Fi +/// @param A +/// @param B +/// @param ldb +template +inline void RRRExpand (const Field& Fi, + const RRRgen* A, + typename Field::Element_ptr B, const size_t ldb) + { + // leaf + if (A->left == nullptr){ + size_t N = A->size_N1; + FFLAS::fassign(Fi, N, N, A->LU_right->PL, A->LU_right->ldPL, B, ldb); + return; + } + + size_t N1 = A->size_N1; + size_t N2 = A->size_N2; + + // cut the matrix A in four parts + // in case of odd N, A12 and A21 are rectangles and not squares + typename Field::Element_ptr B11 = B; + typename Field::Element_ptr B12 = B + N1; + typename Field::Element_ptr B21 = B + ldb*N1; + typename Field::Element_ptr B22 = B21 + N1; + + + // B11 < RRRExpand(A11) + RRRExpand(Fi, A->left, B11, ldb); + + // B22 < RRRExpand(A22) + RRRExpand(Fi, A->right, B22, ldb); + + // B12 < L_u * U_u + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + N1, N2, A->LU_right->r, Fi.one, + A->LU_right->PL, A->LU_right->ldPL, + A->LU_right->UQ, A->LU_right->ldUQ, + Fi.zero, B12, ldb); + + + // B21 < L_l * U_l + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + N2, N1, A->LU_left->r, Fi.one, + A->LU_left->PL, A->LU_left->ldPL, + A->LU_left->UQ, A->LU_left->ldUQ, + Fi.zero, B21, ldb); + + + } + + + +/// @brief Multiplies two matrices stored as rank revealing factorization. C = minus ? -1 : 1 A*B +/// @tparam Field +/// @param Fi +/// @param A stored with an RRgen +/// @param B stored with an RRgen +template +inline RRgen* RRxRR (const Field& Fi,RRgen* A, RRgen*B, bool minus) + { + size_t m = B->m; + size_t n = A->n; + size_t k = A->m; + + // X < UA * LB + typename Field::Element_ptr X = FFLAS::fflas_new(Fi, A->r, B->r); + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + A->r, B->r, k, Fi.one, + A->UQ, A->ldUQ, + B->PL, B->ldPL, + Fi.zero, X, B->r); + + RRgen* RR_X = new RRgen(Fi,A->r,B->r,X,B->r); // X will be modified + FFLAS::fflas_delete(X); + + // LC < LA*LX + typename Field::Element_ptr L_C = FFLAS::fflas_new(Fi, n, RR_X->r); + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + n, RR_X->r, A->r, + minus ? Fi.mOne : Fi.one, A->PL, A->ldPL, + RR_X->PL,RR_X->ldPL, + Fi.zero, L_C, RR_X->r); + + + // RC < RX*RB + typename Field::Element_ptr R_C = FFLAS::fflas_new(Fi, RR_X->r, m); + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + RR_X->r, m, B->r, Fi.one, + RR_X->UQ, RR_X->ldUQ, + B->UQ, B->ldUQ, + Fi.zero, R_C, m); + + + RRgen* C = new RRgen(Fi,n,m,RR_X->r,L_C,RR_X->r,R_C,m,true); // can be modifed + delete(RR_X); + + return C; + } + +/// @brief multiplies two matrices stored as rank revealing factorization. +/// @tparam Field +/// @param Fi +/// @param A stored with an RRgen +/// @param B stored with an RRgen +template +inline RRgen* RRxRR (const Field& Fi,RRgen* A, RRgen*B){ + return RRxRR (Fi, A, B, false); +} + + +/// @brief add two matrices stored as rank revealing factorization +/// @tparam Field +/// @param Fi +/// @param A stored with an RRgen +/// @param B stored with an RRgen +template +inline RRgen* RRaddRR (const Field& Fi, RRgen* A, RRgen*B) + { + // X < [LA LB] + typename Field::Element_ptr X = FFLAS::fflas_new(Fi, A->n, A->r + B->r); + FFLAS::fassign(Fi,A->n,A->r,A->PL,A->ldPL,X,A->r+ B->r); + typename Field::Element_ptr X2 = X + A->r; + FFLAS::fassign(Fi,B->n,B->r,B->PL,B->ldPL,X2,A->r+B->r); + + // Y < [RA] + // [RB] + + typename Field::Element_ptr Y = FFLAS::fflas_new(Fi, A->r + B->r, A->m); + + FFLAS::fassign(Fi,A->r,A->m,A->UQ,A->ldUQ,Y,A->m); + typename Field::Element_ptr Y2 = Y+(A->r)*(A->m); + FFLAS::fassign(Fi,B->r,B->m,B->UQ,B->ldUQ,Y2,B->m); + + + // LX,RX < facto(X) + RRgen* X_fact = new RRgen(Fi, A->n, A->r + B->r, X, A->r + B->r); + + // LY,RY < facto(Y) + RRgen* Y_fact = new RRgen(Fi, A->r + B->r, B->m, Y, B->m); + + + // D < RRxRR(X,Y) + RRgen* D = RRxRR(Fi,X_fact,Y_fact); + + FFLAS::fflas_delete(Y); + FFLAS::fflas_delete(X); + delete(X_fact); + delete(Y_fact); + + return D; + } + +/// @brief Adds a quasiseparable matrix in RRR representation and a rank revealing factorization. +/// @tparam Field +/// @param Fi +/// @param A size n*n in RRR representation +/// @param B size n*n in RR representation +template +inline RRRgen* RRRaddRR (const Field& Fi, RRRgen* A, RRgen* B, bool same_threshold = false) +{ + size_t threshold = A->t + B->r; + if (same_threshold){ + threshold = A->t; + } + // std::cout << "n et rank " <size_N1+A->size_N2 << " " << A->t + B->r << std::endl; + if (A->size_N1+A->size_N2 <= threshold ){ + + // C = RRRexpand(A) + typename Field::Element_ptr C = FFLAS::fflas_new(Fi, A->size_N1+A->size_N2, A->size_N1+A->size_N2); + RRRExpand(Fi, A, C, A->size_N1 + A->size_N2); + + // B_expanded = RRexpand(B) + + typename Field::Element_ptr B_expanded = FFLAS::fflas_new(Fi,B->n,B->m); + B->RRExpand(Fi,B_expanded,B->m); + + // C = B+C + FFLAS::faddin(Fi,A->size_N1+ A->size_N2,A->size_N1 +A->size_N2,B_expanded,B->n,C,A->size_N1 +A->size_N2); + FFLAS::fflas_delete(B_expanded); + RRRgen* D = new RRRgen(Fi,C,A->size_N1 + A->size_N2,threshold,true,true); + FFLAS::fflas_delete(C); + + return D; + } + else { + // B_expanded = [RR_B11 RR_B12] + // [RR_B21 RR_B22] + + RRgen* RR_B11 = new RRgen(Fi, A->size_N1, A->size_N1, B->r, B->PL, B->ldPL, B->UQ,B->ldUQ,false); + RRgen* RR_B12 = new RRgen(Fi, A->size_N1, A->size_N2, B->r, B->PL, B->ldPL, B->UQ + A->size_N1,B->ldUQ,false); + RRgen* RR_B21 = new RRgen(Fi, A->size_N2, A->size_N1, B->r, B->PL + (A->size_N1 * B->ldPL), B->ldPL, B->UQ ,B->ldUQ,false); + RRgen* RR_B22 = new RRgen(Fi, A->size_N2, A->size_N2, B->r, B->PL + (A->size_N1 * B->ldPL), B->ldPL, B->UQ + A->size_N1,B->ldUQ,false); + + + + // C11 = RRR+RR(A11,B11) + RRRgen* C11 = RRRaddRR(Fi, A->left,RR_B11,same_threshold); + + // C22 = RRR+RR(A22,B22) + RRRgen* C22 = RRRaddRR(Fi, A->right,RR_B22,same_threshold); + + // C12 = RR+RR(A12,B12) + RRgen* C12 = RRaddRR(Fi, A->LU_right,RR_B12); + + // C21 = RR+RR(A21,B21) + RRgen* C21 = RRaddRR(Fi, A->LU_left,RR_B21); + + delete RR_B11; + delete RR_B12; + delete RR_B21; + delete RR_B22; + + return new RRRgen(C12,C21,A->size_N1,A->size_N2,threshold,C11,C22,true); + } +} + + + +/// @brief Multiplies a quasiseparable matric in RRR representation with a tall and skinny matrix +/// @tparam Field +/// @param Fi +/// @param n +/// @param t +/// @param A size n*n in RRR representation +/// @param B size n*t +/// @param ldB leading dimension of B +/// @param C size n*t +/// @param ldC leading dimension of C +template +inline void RRRxTS (const Field& Fi, size_t n, size_t m_B, + const RRRgen* A, typename Field::ConstElement_ptr B, size_t ldB, + typename Field::Element_ptr C, size_t ldC) + { + if (n<=A->t+m_B){ + typename Field::Element_ptr Adense =FFLAS::fflas_new (Fi, n, n); + RRRExpand(Fi, A, Adense, n); + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + n, m_B, n, + Fi.one, Adense, n, + B, ldB, + 0, C, ldC); + FFLAS::fflas_delete(Adense); + return; + } + + else { + + size_t N1 = A->size_N1; + size_t N2 = A->size_N2; + + // split the matrices as [C1] = [A11 A12] [B1] + // [C2] = [A12 A22] [B2] + typename Field::Element_ptr C1 = C; + typename Field::Element_ptr C2 = C1 + N1*ldC; + + typename Field::ConstElement_ptr B1 = B; + typename Field::ConstElement_ptr B2 = B1 + N1*ldB; + + // C1 < RRRxTS(A11,B1) + RRRxTS(Fi, N1, m_B, A->left, B1, ldB, C1, ldC); + + // C2 < RRRxTS(A22,B2) + RRRxTS(Fi, N2, m_B, A->right, B2, ldB, C2, ldC); + + // X < RA12 x B2 + // X of size ru*m_B + typename Field::Element_ptr X = FFLAS::fflas_new (Fi, A->LU_right->r, m_B); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + A->LU_right->r, m_B, N2, + Fi.one, A->LU_right->UQ, A->LU_right->ldUQ, + B2, ldB, + 0, X, m_B); + + // C1 < C1 + LA12 x X + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + N1, m_B , A->LU_right->r , + Fi.one, A->LU_right->PL, A->LU_right->ldPL, + X, m_B, + Fi.one, C1, ldC); + + FFLAS::fflas_delete(X); + + // Y < RA21 x B1 + typename Field::Element_ptr Y = FFLAS::fflas_new (Fi, A->LU_left->r, m_B); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + A->LU_left->r, m_B, N1, + Fi.one,A->LU_left->UQ, A->LU_left->ldUQ, + B1, ldB, + 0, Y, m_B); + + // C2 < C2 + LA21 x Y + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + N2, m_B, A->LU_left->r, + Fi.one, A->LU_left->PL, A->LU_left->ldPL, + Y, m_B, + Fi.one, C2, ldC); + + FFLAS::fflas_delete(Y); + // C = [C1] + // [C2] + + } + } + + +/// @brief Multiplies a tall and skinny matrix with a quasiseparable matric in RRR representation +/// @tparam Field +/// @param Fi +/// @param n +/// @param t +/// @param B size t*n +/// @param ldB leading dimension of B +/// @param A size n*n in RRR representation +/// @param C size t*n +/// @param ldC leading dimension of C +template +inline void TSxRRR (const Field& Fi, size_t n, size_t t, + typename Field::ConstElement_ptr B, size_t ldB, const RRRgen* A, + typename Field::Element_ptr C, size_t ldC) + { + if (n<=A->t+t){ + typename Field::Element_ptr Adense = FFLAS::fflas_new (Fi, n, n); + RRRExpand(Fi, A, Adense, n); + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t, n, n, + Fi.one,B, ldB, + Adense, n, + 0, C, ldC); + FFLAS::fflas_delete(Adense); + return; + } + + else { + + size_t N1 = A->size_N1; + size_t N2 = A->size_N2; + + // split the matrices as [A11 A12] + // [C1 C2] = [B1 B2][A12 A22] + typename Field::Element_ptr C1 = C; + typename Field::Element_ptr C2 = C1 + N1; + + typename Field::Element_ptr B1 = (typename Field::Element_ptr)B; + typename Field::Element_ptr B2 = B1 + N1; + + // C1 < TSxRRR(A11,B1) + TSxRRR(Fi, N1, t, B1, ldB, A->left, C1, ldC); + + // C2 < RRRxTS(A22,B2) + TSxRRR(Fi, N2, t, B2, ldB, A->right, C2, ldC); + + + // X < B2*A21->PL + // X of size t*r + typename Field::Element_ptr X = FFLAS::fflas_new (Fi, t, A->LU_left->r); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t, A->LU_left->r, N2, + Fi.one,B2, ldB , + A->LU_left->PL, A->LU_left->ldPL, + 0, X,A->LU_left->r); + + // C1 < C1 + X*A21->UQ + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t, N1 , A->LU_left->r , + Fi.one, X, A->LU_left->r, + A->LU_left->UQ, A->LU_left->ldUQ, + Fi.one, C1, ldC); + + FFLAS::fflas_delete(X); + + // Y < B1*A12->PL + // Y of size t*r + typename Field::Element_ptr Y = FFLAS::fflas_new (Fi, t, A->LU_right->r); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t, A->LU_right->r, N1, + Fi.one,B1, ldB , + A->LU_right->PL, A->LU_right->ldPL, + 0, Y,A->LU_right->r); + + // C2 < C2 + Y*A12->UQ + fgemm(Fi, + FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t, N2 , A->LU_right->r , + Fi.one, Y, A->LU_right->r, + A->LU_right->UQ, A->LU_right->ldUQ, + Fi.one, C2, ldC); + + FFLAS::fflas_delete(Y); + + + // C = [C1 C2] + + } + } + + + +/// @brief Multiplies a QS matrix in RRR representation with a rank revealing factorization +/// Computes C = A*B if minus == false , else computes C = -A*B +/// @tparam Field +/// @param Fi +/// @param A size n*n in a RRR representation +/// @param B size n*m in RR representation +template +inline RRgen* RRxRRR (const Field& Fi,const RRRgen* A,const RRgen* B,bool minus) + { + + // X < TSxRRR(RB,A) + size_t n = A->size_N1+A->size_N2; + if(B->m != n){ + std::cout << "PAS LE BON FORMAT " << std::endl; + return nullptr; + } + + if (B->UQ == nullptr){ + // if B is a leaf of the RRRgen then UQ is empty ! + // and the dense matrix is stored in PL + typename Field::Element_ptr X = FFLAS::fflas_new (Fi, B->n, n); + TSxRRR(Fi,n,B->n,B->PL,B->ldPL,A,X,n); + RRgen* RR_X = new RRgen(Fi, B->r, n, X, n); + FFLAS::fflas_delete(X); + return RR_X; + } + + typename Field::Element_ptr X = FFLAS::fflas_new (Fi, B->r, n); + TSxRRR(Fi,n,B->r,B->UQ,B->ldUQ,A,X,n); + + // (LX, RX) < RRF(X) + RRgen* RR_X = new RRgen(Fi, B->r, n, X, n); + FFLAS::fflas_delete(X); + + // RD < RX + typename Field::Element_ptr UQ_D = FFLAS::fflas_new (Fi, RR_X->r, n); + FFLAS::fassign(Fi,RR_X->r,n,RR_X->UQ,RR_X->ldUQ,UQ_D,n); + + // LD < (minus ? (-1) : 1) LBxLX + typename Field::Element_ptr PL_D = FFLAS::fflas_new (Fi, B->n, RR_X->r); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + B->n, RR_X->r, B->r, + minus ? Fi.mOne : Fi.one , B->PL, B->ldPL, + RR_X->PL, RR_X->ldPL, + Fi.zero, PL_D, RR_X->r); + + + RRgen* D = new RRgen(Fi,B->n,n,RR_X->r,PL_D,RR_X->ldPL,UQ_D,B->m,true); + delete RR_X; + return D; + + } + +/// @brief (algo 9) Multiplies a QS matrix in RRR representation with a rank revealing factorization +/// Computes C = A*B if minus == false , else computes C = -A*B +/// @tparam Field +/// @param Fi +/// @param A size n*n in a RRR representation +/// @param B size n*m in RR representation +/// @param minus boolean in order to compute C = -A*B if true. +template +inline RRgen* RRRxRR (const Field& Fi,const RRRgen* A,const RRgen* B,bool minus) + { + // X < RRRxTS(A,LB) + size_t n = A->size_N1+A->size_N2; + typename Field::Element_ptr X = FFLAS::fflas_new (Fi, n, B->r); + RRRxTS(Fi,n,B->r,A,B->PL,B->ldPL,X,B->r); + + // (LX, RX) < RRF(X) + RRgen* RR_X = new RRgen(Fi, n, B->r, X, B->r); + FFLAS::fflas_delete(X); + if (B->UQ == nullptr){ + return RR_X; + } + // LD < LX + typename Field::Element_ptr PL_D = FFLAS::fflas_new (Fi, n, RR_X->r); + FFLAS::fassign(Fi,n,RR_X->r,RR_X->PL,RR_X->ldPL,PL_D,RR_X->r); + + // RD < (minus ? (-1) : 1) RX x RB + typename Field::Element_ptr UQ_D = FFLAS::fflas_new (Fi, RR_X->r, B->m); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + RR_X->r, B->m, B->r, + minus ? Fi.mOne : Fi.one , RR_X->UQ, RR_X->ldUQ, + B->UQ, B->ldUQ, + Fi.zero, UQ_D, B->m); + + + RRgen* D = new RRgen(Fi,n,B->m,RR_X->r,PL_D,RR_X->ldPL,UQ_D,B->m,true); + delete RR_X; + return D; + } + + +/// @brief Multiplies a QS matrix in RRR representation with a rank revealing factorization +/// Computes C = A*B +/// @tparam Field +/// @param Fi +/// @param A size n*n in a RRR representation +/// @param B size n*m in RR representation +template +inline RRgen* RRRxRR (const Field& Fi,const RRRgen* A,const RRgen* B){ + return RRRxRR(Fi,A,B,false); +} +/// @brief (algo 9) Multiplies a QS matrix in RRR representation with a rank revealing factorization +/// Computes C = A*B +/// @tparam Field +/// @param Fi +/// @param A size n*n in a RRR representation +/// @param B size n*m in RR representation +template +inline RRgen* RRxRRR (const Field& Fi,const RRRgen* A,const RRgen* B){ + return RRxRRR(Fi,A,B,false); +} + +/// @brief multiplies two QS matrices +/// @tparam Field +/// @param Fi +/// @param A size n*n +/// @param B size n*n +template +inline RRRgen* RRRxRRR (const Field& Fi, const RRRgen* A, const RRRgen* B){ + size_t n = A->size_N1+A->size_N2; + if (n<= A->t+B->t){ + //return RRRExpand(A) x RRRExpand(B) + typename Field::Element_ptr A_expanded =FFLAS::fflas_new (Fi, n, n); + typename Field::Element_ptr B_expanded =FFLAS::fflas_new (Fi, n, n); + typename Field::Element_ptr C = FFLAS::fflas_new (Fi, n, n); + RRRExpand(Fi,A,A_expanded,n); + RRRExpand(Fi,B,B_expanded,n); + fgemm(Fi, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + n, n, n, + Fi.one, A_expanded, n, + B_expanded, n, + 0, C, n); + FFLAS::fflas_delete(A_expanded); + FFLAS::fflas_delete(B_expanded); + RRRgen* D = new RRRgen(Fi,n,A->t+B->t,C,n,true,true); + FFLAS::fflas_delete(C); + return D; + } + + else { + // C11 < RRRxRRR(A11,B11) + RRRgen* C11_ = RRRxRRR(Fi,A->left,B->left); + + // C22 < RRRxRRR(A22,B22) + RRRgen* C22_ = RRRxRRR(Fi,A->right,B->right); + + // X < RRxRR(A12,B21) + RRgen* X = RRxRR(Fi,A->LU_right,B->LU_left); + + // Y < RRxRR(A21,B12) + RRgen* Y = RRxRR(Fi,A->LU_left,B->LU_right); + + // C11 < RRR+RR(C11,X) + RRRgen* C11 = RRRaddRR(Fi,C11_,X); + + // C22 < RRR+RR(C22,Y) + RRRgen* C22 = RRRaddRR(Fi,C22_,Y); + + delete C11_; + delete C22_; + delete X; + delete Y; + + // LX < RRRxTS(A11,LB12) ; RX < RB12 + typename Field::Element_ptr LX =FFLAS::fflas_new (Fi, A->size_N1, B->LU_right->r); + RRRxTS(Fi,A->size_N1,B->LU_right->r,A->left,B->LU_right->PL,B->LU_right->ldPL,LX,B->LU_right->r); + + typename Field::Element_ptr RX =FFLAS::fflas_new (Fi, B->LU_right->r, B->size_N2); + FFLAS::fassign(Fi,B->LU_right->r,B->size_N2, + B->LU_right->UQ,B->LU_right->ldUQ, + RX,B->size_N2); + + // LY < LA12 ; RY < TSxRRR(RA12,B22) + typename Field::Element_ptr LY =FFLAS::fflas_new (Fi, A->size_N1, A->LU_right->r); + FFLAS::fassign(Fi,A->size_N1,A->LU_right->r, + A->LU_right->PL,A->LU_right->ldPL, + LY,A->LU_right->r); + + typename Field::Element_ptr RY =FFLAS::fflas_new (Fi, A->LU_right->r, B->size_N2); + TSxRRR(Fi,B->size_N2,A->LU_right->r,A->LU_right->UQ,A->LU_right->ldUQ,B->right,RY,B->size_N2); + + X = new RRgen(Fi, A->size_N1, B->size_N2,B->LU_right->r,LX,B->LU_right->r,RX,B->size_N2,true); + Y = new RRgen(Fi, A->size_N1, B->size_N2,A->LU_right->r,LY,A->LU_right->r,RY,B->size_N2,true); + + // C12 < RR+RR(X,Y) + RRgen* C12 = RRaddRR(Fi,X,Y); + delete X; + delete Y; + + // LX < RRRxTS(A22,LB21) ; RX < RB21 + LX = FFLAS::fflas_new (Fi, A->size_N2, B->LU_left->r); + RRRxTS(Fi,A->size_N2,B->LU_left->r,A->right,B->LU_left->PL,B->LU_left->ldPL,LX,B->LU_left->r); + + RX = FFLAS::fflas_new (Fi, B->LU_left->r, B->size_N1); + FFLAS::fassign(Fi,B->LU_left->r,B->size_N1, + B->LU_left->UQ,B->LU_left->ldUQ, + RX,B->size_N1); + + // LY < LA21 ; RY < TSxRRR(RA21,B11) + LY =FFLAS::fflas_new (Fi, A->size_N2, A->LU_left->r); + FFLAS::fassign(Fi,A->size_N2,A->LU_left->r, + A->LU_left->PL,A->LU_left->ldPL, + LY,A->LU_left->r); + + RY = FFLAS::fflas_new (Fi, A->LU_left->r, B->size_N1); + TSxRRR(Fi,B->size_N1,A->LU_left->r,A->LU_left->UQ,A->LU_left->ldUQ,B->left,RY,B->size_N1); + + X = new RRgen(Fi, A->size_N2, B->size_N1,B->LU_left->r,LX,B->LU_left->r,RX,B->size_N1,true); + Y = new RRgen(Fi, A->size_N2, B->size_N1,A->LU_left->r,LY,A->LU_left->r,RY,B->size_N1,true); + + // C21 < RR+RR(X,Y) + RRgen* C21 = RRaddRR(Fi,X,Y); + delete X; + delete Y; + + // RETURN C = [C11 C12] + // [C12 C21] + return new RRRgen(C12,C21,A->size_N1,B->size_N2,A->t+B->t,C11,C22,true); + } + +} + + +/// @brief Computes the inverse in RRR representation adn returns it in RRRgen. +/// @tparam Field +/// @param Fi +/// @param A in RRR representation +template +inline RRRgen* RRRinvert (const Field& Fi,const RRRgen* A){ + if (!A->left){ + // leaf + // Y < RRRExpand(A) + size_t N1 = A->size_N1; + typename Field::Element_ptr Y =FFLAS::fflas_new (Fi, N1, N1); + RRRExpand(Fi,A,Y,N1); + // return Invert(Y) + int nullity; + FFPACK::Invert (Fi, N1,Y, N1, nullity); + if (nullity != 0) { + FFLAS::fflas_delete(Y); + throw std::runtime_error("RRR Error: non invertible matrix "); + } + RRRgen* Y_RRR = new RRRgen(Fi,N1,A->t,Y,N1,true,true); + FFLAS::fflas_delete(Y); + + return Y_RRR; + } + + // split the matrix as A = [A11 A12] and X = [X11 X12] + // [A21 A22] [X21 X22] + + // Y11 < RRRinvertrec(A11) + RRRgen* Y11 = nullptr; + try { + Y11 = RRRinvert(Fi,A->left); + } + catch(...){ + delete Y11; + throw std::runtime_error("RRR Error: non invertible matrix "); + } + + // Y12 < RRRxRR(Y11,A12) + RRgen* Y12 = RRRxRR(Fi,Y11,A->LU_right); + + // Y21 < RRxRRR(A21,Y11) + RRgen* Y21 = RRxRRR(Fi,Y11,A->LU_left); + + // Z < -RRxRR(A21,Y12) + RRgen* Z = RRxRR(Fi,A->LU_left,Y12,true); + + // D < RRRaddRR(A22,Z) + RRRgen* D = RRRaddRR(Fi,A->right,Z); + delete Z; + RRRgen* X22 = nullptr; + // X22 < RRRinvert(D) + try { + X22 = RRRinvert(Fi,D); + } + catch(...){ + delete Y11; + delete Y12; + delete Y21; + delete D; + delete X22; + throw std::runtime_error("RRR Error: non invertible matrix "); + } + delete D; + + // X21 < -RRRxRR(X22,Y21) + RRgen* X21 = RRRxRR(Fi,X22,Y21,true); + + // W < -RRxRR(Y12,X21) + RRgen* W = RRxRR(Fi,Y12,X21,true); + + // X12 < -RRxRRR(Y12,X22) + RRgen* X12 = RRxRRR(Fi,X22,Y12,true); + + // X11 < RRRaddRR(Y11,W) + RRRgen* X11 = RRRaddRR(Fi, Y11, W); + delete W; + delete Y11; + delete Y12; + delete Y21; + + // return X + return new RRRgen(X12,X21,A->size_N1,A->size_N2,A->t,X11,X22,true); +} + + +/// @brief Computes B = U^-1*B. +/// @tparam Field +/// @param Fi +/// @param U in RRR representation, must be inversible +/// @param B TS n*t +/// @param ldb +/// @param t number of columns of B +template +inline void TRSM_RRR_TS (const Field& Fi,const FFLAS::FFLAS_SIDE Side, const FFLAS::FFLAS_UPLO Uplo, const RRRgen* M, typename Field::Element_ptr B, size_t ldb, size_t t) +{ + size_t size_N1 = M->size_N1; + if (!(M->left)){ + if (Side == FFLAS::FflasRight){ + if (M->LU_right->UQ){ + size_t* Q = reinterpret_cast(M->LU_right->UQ); + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasTrans, size_N1, 0, size_N1, M->LU_right->PL, M->LU_right->ldPL, Q); + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasTrans, t, 0, size_N1, B, ldb, Q); + FFLAS::ftrsm(Fi,Side, Uplo,FFLAS::FflasNoTrans,FFLAS::FflasNonUnit, t, M->size_N1, Fi.one,M->LU_right->PL,M->LU_right->ldPL,B,ldb); + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasNoTrans, size_N1, 0, size_N1, M->LU_right->PL, M->LU_right->ldPL, Q); + } + else{ + FFLAS::ftrsm(Fi,Side, Uplo,FFLAS::FflasNoTrans,FFLAS::FflasNonUnit, t, M->size_N1, Fi.one,M->LU_right->PL,M->LU_right->ldPL,B,ldb); + } + } + else { + if (M->LU_right->UQ){ + size_t* Q = reinterpret_cast(M->LU_right->UQ); + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasTrans, size_N1, 0, size_N1, M->LU_right->PL, M->LU_right->ldPL, Q); + FFLAS::ftrsm(Fi,Side, Uplo,FFLAS::FflasNoTrans,FFLAS::FflasNonUnit, M->size_N1,t,Fi.one,M->LU_right->PL,M->LU_right->ldPL,B,ldb); + applyP(Fi, FFLAS::FflasLeft, FFLAS::FflasTrans, size_N1, 0, size_N1, B, ldb, Q); + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasNoTrans, size_N1, 0, size_N1, M->LU_right->PL, M->LU_right->ldPL, Q); + } + else { + FFLAS::ftrsm(Fi,Side, Uplo,FFLAS::FflasNoTrans,FFLAS::FflasNonUnit, M->size_N1,t,Fi.one,M->LU_right->PL,M->LU_right->ldPL,B,ldb); + } + } + return; + } + + + size_t size_N2 = M->size_N2; + + + // B = M^-1*B. + if ((Side == FFLAS::FflasLeft && Uplo == FFLAS::FflasUpper) || (Side == FFLAS::FflasRight && Uplo == FFLAS::FflasLower)) + { + if (Uplo == FFLAS::FflasUpper){ + // B2 = TRSM_RRR_TS(M2,B2) + typename Field::Element_ptr B2 = B + size_N1 * ldb; + TRSM_RRR_TS(Fi, Side, Uplo, M->right,B2,ldb,t); + size_t r = M->LU_right->r; + // B1 = B1 - M12*B2 + typename Field::Element_ptr M_int = FFLAS::fflas_new(Fi,r,t); + + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + r,t,size_N2,Fi.one, + M->LU_right->UQ,M->LU_right->ldUQ, + B2,ldb, + Fi.zero, M_int, t); + + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + size_N1,t,r,Fi.mOne, + M->LU_right->PL,M->LU_right->ldPL, + M_int,t, + Fi.one, B, ldb); + FFLAS::fflas_delete(M_int); + // B1 = TRSM_RRR_TS(M1,B1) + TRSM_RRR_TS(Fi, Side, Uplo, M->left,B,ldb,t); + } + + else { + // B2 = TRSM_RRR_TS(M2,B2) + size_t r = M->LU_left->r; + typename Field::Element_ptr B2 = B + size_N1; + TRSM_RRR_TS(Fi, Side, Uplo, M->right,B2,ldb,t); + + // B1 = B1 - B2*M21 + typename Field::Element_ptr M_int = FFLAS::fflas_new(Fi,t,r); + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t,r,size_N2,Fi.one, + B2,ldb, + M->LU_left->PL,M->LU_left->ldPL, + Fi.zero, M_int, r); + + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t,size_N1,r,Fi.mOne, + M_int,r, + M->LU_left->UQ,M->LU_left->ldUQ, + Fi.one, B, ldb); + FFLAS::fflas_delete(M_int); + TRSM_RRR_TS(Fi, Side, Uplo, M->left,B,ldb,t); + } + } + + + else{ + // B1 = TRSM_RRR_TS(M1,B1) + typename Field::Element_ptr B2; + TRSM_RRR_TS(Fi, Side, Uplo, M->left,B,ldb,t); + // B2 = B2 - (B1*M12 / M21*B1) + + + // B2 = B2 - B1*M12 + if (Uplo == FFLAS::FflasUpper){ + size_t r = M->LU_right->r; + B2 = B + size_N1; + typename Field::Element_ptr M_int = FFLAS::fflas_new(Fi,t,r); + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t,r,size_N1,Fi.one, + B,ldb, + M->LU_right->PL,M->LU_right->ldPL, + Fi.zero, M_int, r); + + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + t,size_N2,r,Fi.mOne, + M_int,r, + M->LU_right->UQ,M->LU_right->ldUQ, + Fi.one, B2, ldb); + FFLAS::fflas_delete(M_int); + } + + + + // B2 = B2 - M21*B1 + else { + size_t r = M->LU_left->r; + B2 = B + size_N1*ldb; + typename Field::Element_ptr M_int = FFLAS::fflas_new(Fi,r,t); + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + r,t,size_N1,Fi.one, + M->LU_left->UQ,M->LU_left->ldUQ, + B,ldb, + Fi.zero, M_int, t); + + fgemm(Fi,FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + size_N2,t,r,Fi.mOne, + M->LU_left->PL,M->LU_left->ldPL, + M_int,t, + Fi.one, B2, ldb); + FFLAS::fflas_delete(M_int); + } + // B2 = TRSM_RRR_TS(M2,B2) + TRSM_RRR_TS(Fi, Side, Uplo, M->right,B2,ldb,t); + } + +} + +/// @brief Computes X = M^-1*B or X = B*M^-1 +/// @tparam Field +/// @param Fi +/// @param Side Side to apply TRSM +/// @param Uplo Whether the matrix M is lower triangular or upper triangular +/// @param U in RRR representation and triangular, must be inversible +/// @param B RRgen +template +inline void TRSM_RRR_RR (const Field& Fi,const FFLAS::FFLAS_SIDE Side, const FFLAS::FFLAS_UPLO Uplo, const RRRgen* M, RRgen* B) +{ + // if X = M^-1*B + if (Side == FFLAS::FflasLeft){ + return TRSM_RRR_TS(Fi,Side,Uplo,M,B->PL,B->ldPL,B->r); + } + // if X = B*M^-1 + return TRSM_RRR_TS(Fi,Side,Uplo,M,B->UQ,B->ldUQ,B->r); +} + +/// @brief Computes the L U factorization in RRR representation and returns it in RRRgen. +/// @tparam Field +/// @param Fi +/// @param A in RRR representation +/// @param L RRRgen uninitialized +/// @param U RRRgen uninitialized +template +inline void LUfactRRR (const Field& Fi, const RRRgen* A, RRRgen*& L, RRRgen*& U) +{ + size_t size_N1 = A->size_N1; + size_t size_N2 = A->size_N2; + size_t N = size_N1 + size_N2; + if (N <= A->t){ + size_t* P = FFLAS::fflas_new(size_N1); + size_t* Q = FFLAS::fflas_new(size_N1); + typename Field::Element_ptr A_copy = FFLAS::fflas_new(Fi, size_N1, size_N1); + FFLAS::fassign(Fi, size_N1, size_N1, A->LU_right->PL, A->LU_right->ldPL, A_copy, size_N1); + + size_t r = PLUQ(Fi, FFLAS::FflasNonUnit, size_N1, size_N1, A_copy, size_N1, P, Q); + if (r(Fi, FFLAS::FflasUpper, FFLAS::FflasNonUnit, size_N1, size_N1, r, A_copy, size_N1, UQ, size_N1, true); + + // extraction of L + getTriangular(Fi, FFLAS::FflasLower, FFLAS::FflasUnit, size_N1, size_N1, r, A_copy, size_N1, PL, r, true); + + // apply the permutations (P on L and Q on U) + bool flag = false; + for (size_t i = 0; it,true,false); + U = new RRRgen(Fi,UQ,N,A->t,true,false); + if (flag){ + applyP(Fi, FFLAS::FflasRight, FFLAS::FflasNoTrans, r, 0, size_N1, UQ, size_N1, Q); + U->LU_right->UQ = reinterpret_cast(Q); + } + else { + FFLAS::fflas_delete(Q); + } + FFLAS::fflas_delete(A_copy); + FFLAS::fflas_delete(P); + return; + } + + // L11/U11 = LUfactRRR(A11) + RRRgen* L11 = nullptr; + RRRgen* U11 = nullptr; + try { + LUfactRRR(Fi,A->left,L11,U11); + } + catch(...){ + delete L11; + delete U11; + throw std::runtime_error("RRR Error: LU with non invertible matrix "); + } + + // D2 = L11^{-1}*L12 / U12 + RRgen* D2 = A->LU_right->RRcopy(Fi); + TRSM_RRR_RR(Fi,FFLAS::FflasLeft, FFLAS::FflasLower,L11,D2); + + // D1 = L21 / U21xU11^{-1} + RRgen* D1 = A->LU_left->RRcopy(Fi); + TRSM_RRR_RR(Fi,FFLAS::FflasRight, FFLAS::FflasUpper,U11,D1); + + + // X22 = A22 - D1*D2 + RRgen* D1xD2 = RRxRR(Fi,D1,D2,true); + RRRgen* X22 = RRRaddRR(Fi,A->right,D1xD2,true); + delete D1xD2; + + + // L2/U2 = LUfactRRR(X22) + RRRgen* L22 = nullptr; + RRRgen* U22 = nullptr; + try { + LUfactRRR(Fi,X22,L22,U22); + } + catch(...){ + delete L11; + delete U11; + delete L22; + delete U22; + delete X22; + delete D1; + delete D2; + throw std::runtime_error("RRR Error: LU with non invertible matrix "); + } + delete X22; + + // L = [L11 0] U = [U11 D2] + // [D1 L22] [0 U22] + typename Field::Element_ptr L12_PL = FFLAS::fflas_new(Fi, size_N1, 1); + FFLAS::fzero(Fi, size_N1, 1,L12_PL,1); + typename Field::Element_ptr L12_UQ = FFLAS::fflas_new(Fi, 1, size_N2); + FFLAS::fzero(Fi,1, size_N2,L12_UQ,size_N2); + RRgen* L12 = new RRgen(Fi,size_N1,size_N2,1, L12_PL,1,L12_UQ,size_N2,true); + L = new RRRgen(L12,D1,size_N1,size_N2,A->t,L11,L22,true); + + + + typename Field::Element_ptr U21_PL = FFLAS::fflas_new(Fi, size_N2, 1); + FFLAS::fzero(Fi,size_N2, 1,U21_PL,1); + typename Field::Element_ptr U21_UQ = FFLAS::fflas_new(Fi, 1, size_N1); + FFLAS::fzero(Fi,1, size_N1,U21_UQ,size_N1); + RRgen* U21 = new RRgen(Fi,size_N2,size_N1,1, U21_PL,1,U21_UQ,size_N1,true); + U = new RRRgen(D2,U21,size_N1,size_N2,A->t,U11,U22,true); +} + + + + + +} +#endif //_FFPACK_ffpack_rrrgen_inl \ No newline at end of file diff --git a/fflas-ffpack/utils/fflas_randommatrix.h b/fflas-ffpack/utils/fflas_randommatrix.h index a5bce3f3..57b56c5b 100644 --- a/fflas-ffpack/utils/fflas_randommatrix.h +++ b/fflas-ffpack/utils/fflas_randommatrix.h @@ -467,6 +467,264 @@ namespace FFPACK{ FFLAS::fflas_delete (randcols,randrows); } + // O(l*len(L)+M*len(M)) + // compute all the valid moves in the array valid_moves depending on the free areas described in L and M + inline void compute_valid_moves(std::vector* valid_moves, size_t* start_up, size_t l, size_t* L, size_t lengthL, int llim, size_t m, size_t* M, size_t lengthM, int mlim, std::vector availrows, std::vector availcols) { + valid_moves->clear(); + + // Base case: Identity + valid_moves->push_back(l); + + // Compute all the empty columns/lines + for (size_t j = (size_t) llim+1; j < l; ++j) { + if (availcols[j]) { + valid_moves->push_back(j); + } + } + + // start_up point where the valid moves change direction + *start_up = valid_moves->size(); + + for (size_t i = (size_t) mlim+1; i < m; ++i) { + if (availrows[i]) { + valid_moves->push_back(i); + } + } + } + + // O(1) + // estimate what rank a submatrix is able to reach given the current configuration of the algorithm + size_t valuate(int ti, int r, int n, int i, int Is_pivot_not_moved,int nb_pivot_after_i, int nb_pivot_bf_i, int nb_pivot_after_i_not_moved, int nb_pivot_before_i_not_moved) { + // Compute the valuation of a given leading submatrix ti + return ti + std::min(i + 1 - nb_pivot_bf_i, nb_pivot_after_i_not_moved - Is_pivot_not_moved) + + std::min(n - (i + 1) - nb_pivot_after_i, nb_pivot_before_i_not_moved - Is_pivot_not_moved); + } + + /** @brief genenration of an Rank Profile Matrix randomly + * @param n Dimension + * @param r rank of the n*n matrix + * @param t quasi-separability order + * @param rows + * @param cols + */ + void RandomLTQSRankProfileMatrix_Tom(size_t n, size_t r, size_t t, size_t* rows, size_t* cols) { + + if (r <= 0 || n <= 0 || t <= 0) { + std::cout << "care, all the arguments must be positive" << std::endl; + return; + } + + if (r >= n) { + std::cout << "care, you called the function with invalid arguments : r >= n IMPOSSIBLE" << std::endl; + return; + } + + if (t > r) { + std::cout << "care, you called the function with invalid arguments : t >= r IMPOSSIBLE" << std::endl; + return; + } + + if (r + t > n) { + std::cout << "care, you called the function with invalid arguments : r+t > n IMPOSSIBLE" << std::endl; + return; + } + + /// Initialisation + + size_t* H = FFLAS::fflas_new(n-1); + size_t* T = FFLAS::fflas_new(n-1); + size_t* nb_pivot_before_index = FFLAS::fflas_new(n-1); // pivots that are on lines before line i + size_t* nb_pivot_after_index = FFLAS::fflas_new(n-1); // pivots that are on columns before the last column of the line i + size_t* nb_pivot_before_index_not_moved = FFLAS::fflas_new(n-1); // same idea with pivots that didn't move yet + size_t* nb_pivot_after_index_not_moved = FFLAS::fflas_new(n-1); + + std::vector availablerows (n-1, true); // indicate wether the row i is available or not + std::vector availablecols (n-1, true); // same idea with col j + + for (size_t i = 0; i valid_moves; + size_t start_up; + int ilim, jlim; + int h = 0; + int prev_h = -1; + + while ((size_t) h < r) { + size_t prev_i = rows[h]; + size_t prev_j = cols[h]; + + size_t new_i = prev_i; + size_t new_j = prev_j; + + + /// computating valid_moves + // complex : 2n + if (prev_h != h) { + // compute impossible ones (those which make the configuration having a greater t than asked) + ilim = -1; + jlim = -1; + int k = rows[h]+1; + while ((size_t) k <= n-2) { + if (T[k] == t) { + jlim = n-k-2; + break; + } + ++k; + } + k = rows[h]-1; + while (k >= 0) { + if (T[k] == t) { + ilim = k; + break; + } + --k; + } + compute_valid_moves(& valid_moves, & start_up, prev_j, cols, r, jlim, prev_i, rows, r, ilim, availablerows, availablecols); + } + + /// choose and make the chosen move + size_t rand_index; + if (valid_moves.size()>1){ + rand_index = RandInt(0, valid_moves.size()); + } + else { + rand_index = 0; + } + + if (rand_index < start_up) { + new_j = valid_moves[rand_index]; + cols[h] = new_j; + // updates of the arrays for the valuation function + for (size_t k = prev_i + 1; k <= prev_i + (prev_j - new_j); ++k) { + T[k] += 1; + nb_pivot_after_index[k] += 1; + } + } else { + new_i = valid_moves[rand_index]; + rows[h] = new_i; + // updates of the arrays for the valuation function + for (size_t k = new_i; k < prev_i; ++k) { + T[k] += 1; + nb_pivot_before_index[k] += 1; + } + } + + // updates of the arrays for the valuation function + // O(n) + for (size_t k = 0; k= prev_i) { + nb_pivot_before_index_not_moved[k]--; + } + } + + H[prev_i] = 0; + + /// valutation + // O(n) + bool undo = true; + for (size_t k = 0; k < n - 1; ++k) { + if (valuate(T[k], r, n, k, H[k], nb_pivot_after_index[k], nb_pivot_before_index[k], nb_pivot_after_index_not_moved[k], nb_pivot_before_index_not_moved[k]) >= t) { + undo = false; + break; + } + } + + /// undo if needed + if (undo) { + if (rand_index >= start_up) { + // undo the move and erase it from the valid moves + rows[h] = prev_i; + valid_moves.erase(valid_moves.begin()+rand_index); + // undo of the arrays for valuation func + for (size_t k = new_i; k < prev_i; ++k) { + T[k] -= 1; + nb_pivot_before_index[k] -= 1; + } + } else { + // undo the move and erase it from the valid moves + cols[h] = prev_j; + valid_moves.erase(valid_moves.begin()+rand_index); + // undo of the arrays for valuation func + start_up--; + for (size_t k = prev_i + 1; k <= prev_i + (prev_j - new_j); ++k) { + T[k] -= 1; + nb_pivot_after_index[k] -= 1; + } + } + + // undo of the arrays for valuation func + // O(n) + for (size_t k = 0; k= prev_i) { + nb_pivot_before_index_not_moved[k]--; + } + } + + H[prev_i] = 0; + + ++loop; + if (prev_h != h) { + prev_h = h; + } + // let's try another move + continue; + } + + /// set up for next loop turn + if (prev_j != new_j){ + availablecols[prev_j] = true; + availablecols[new_j] = false; + } + else if (prev_i != new_i){ + availablerows[prev_i] = true; + availablerows[new_i] = false; + } + ++loop; + ++h; + } + FFLAS::fflas_delete (H); + FFLAS::fflas_delete (T); + FFLAS::fflas_delete(nb_pivot_before_index); + FFLAS::fflas_delete(nb_pivot_after_index); + FFLAS::fflas_delete(nb_pivot_before_index_not_moved); + FFLAS::fflas_delete(nb_pivot_after_index_not_moved); + } + + /** @brief Random Matrix with prescribed rank and rank profile matrix * Creates an \c m x \c n matrix with random entries and rank \c r. * @param F field @@ -775,7 +1033,7 @@ namespace FFPACK{ size_t * pivot_r = FFLAS::fflas_new (r); size_t * pivot_c = FFLAS::fflas_new (r); - RandomLTQSRankProfileMatrix (n, r, t, pivot_r, pivot_c); + RandomLTQSRankProfileMatrix_Tom (n, r, t, pivot_r, pivot_c); // typename Field::Element_ptr R =FFLAS::fflas_new(F,n,n); // getLTBruhatGen(F, n, r, pivot_r, pivot_c, R, n); // FFLAS:: WriteMatrix (std::cerr<<"R = "<(r); size_t * cols = FFLAS::fflas_new(r); - RandomLTQSRankProfileMatrix (n, r, t, rows, cols); + RandomLTQSRankProfileMatrix_Tom (n, r, t, rows, cols); typename Field::Element_ptr A = fflas_new(F,n,n); getLTBruhatGen(F, n, r, rows, cols, A, n); diff --git a/tests/test-rrr.C b/tests/test-rrr.C new file mode 100644 index 00000000..bbab7c66 --- /dev/null +++ b/tests/test-rrr.C @@ -0,0 +1,747 @@ +//------------------------------------------------------------------------- +// Test suite for the Quasi-Separable matrices in RRR format +//------------------------------------------------------------------------- + +/* Structure taken from test-quasisep.C */ +#define __FFLASFFPACK_PLUQ_THRESHOLD 5 // Recursive vs iterative PLUQ threshold (default 256) + +#include "fflas-ffpack/fflas-ffpack-config.h" +#include +#include +#include + +#include "fflas-ffpack/fflas/fflas.h" +#include "fflas-ffpack/ffpack/ffpack.h" +#include "fflas-ffpack/utils/test-utils.h" +#include "fflas-ffpack/ffpack/ffpack_rrrgen.inl" +#include +#include "fflas-ffpack/utils/args-parser.h" + +#include + +using namespace std; +using namespace FFPACK; +using namespace FFLAS; + +/** + * \brief test equality between a dense matrix and the result of compressing with RR and reconstructing it. + */ +template +bool test_compression_RR (const Field & F, size_t n, size_t m, + typename Field::Element_ptr A, size_t lda) +{ + typename Field::Element_ptr Acheck = fflas_new (F, n, m); + + + RRgen* RRA = new RRgen(F, n, m, (typename Field::ConstElement_ptr)A, lda); + + RRA->RRExpand(F, Acheck, m); + + bool ok = fequal (F, n, m, A, lda, Acheck, m); + if ( !ok ) + { + std::cout << "ERROR: different results for dense to RR and RR to dense (RRGen and RRExpand)"<t = "<< t<* RRRA = new RRRgen(F, n_A, t, A, lda,true,true); + RRgen* RRB = new RRgen(F, n_A, m_A, (typename Field::ConstElement_ptr)B, ldb); + // std::cout << "B->r = "<< RRB->r<r, RRB->PL, RRB->ldPL); + RRRgen* RRRC = RRRaddRR(F,RRRA,RRB); + + RRRExpand(F, RRRC, C_check, n_A); + + bool ok = fequal (F, n_A, m_A, C_init, m_A, C_check, m_A); + if ( !ok ) + { + std::cout << "ERROR: different results for RRRaddRR and fadd "<* LU_A = new RRgen(F,n,n,(typename Field::ConstElement_ptr)A,n); + + + // L/U check = L/U(A) with RRR_LU + typename Field::Element_ptr Lcheck = fflas_new (F, n, n); + typename Field::Element_ptr Ucheck = fflas_new (F, n, n); + typename Field::Element_ptr Acheck = fflas_new (F, n, n); + + RRRExpand(F, RRRL, Lcheck, n); + RRRExpand(F, RRRU, Ucheck, n); + + fgemm(F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, + n, n,n, + F.one, Lcheck, n, + Ucheck, n, + 0, Acheck, n); + + bool ok = fequal (F, n, n, Acheck, n, A, n); + if ( !ok ) + { + std::cout << "ERROR: different results for dense LU and RRR LU"< > (q,b,n,t,m, r,iters,seed); + ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); + ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); + ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); + ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); + ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); + // ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); // Valgrind does not like this one + // ok = ok &&run_with_field > (q,b,n,t,m, r, iters,seed); + // ok = ok &&run_with_field > (q,9, ceil(n/4.), ceil(t / 4.), ceil(m / 4.), ceil(r / 4.), iters,seed); + // ok = ok &&run_with_field > (q,(b?b:224), ceil(n/4.), ceil(t / 4.), ceil(m / 4.), ceil(r / 4.),iters,seed); + seed++; + } while (loop && ok); + + if (!ok) std::cerr<<"with seed = "<