lattice_extraction.c

Go to the documentation of this file.

/*
 
  $Id: lattice_extraction.c 23065 2016-03-02 09:05:50Z coelho $
 
  Copyright 1989-2016 MINES ParisTech
 
  This file is part of PIPS.
 
  PIPS is free software: you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  any later version.
 
  PIPS 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 General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with PIPS.  If not, see <http://www.gnu.org/licenses/>.
 
*/
#ifdef HAVE_CONFIG_H
    #include "pips_config.h"
#endif
/* HPFC module by Fabien COELHO
 */
 
#include "defines-local.h"
#include "workspace-util.h" // for list_to_base()
#include "prettyprint.h" // for sc_syst_debug()
 
#include "matrix.h"
#include "matrice.h"
#include "sparse_sc.h"
 
/*
void normalize_system(Psysteme * ps)
{
    Ptsg sg;
    Psysteme s;
 
    DEBUG_SYST(1, "entry system", *ps);
 
    sg = sc_to_sg_chernikova(*ps);
    sc_rm(*ps);
 
    s = sg_to_sc_chernikova(sg);
    sg_rm(sg);
 
    DEBUG_SYST(1, "exit system", s);
 
    *ps = s;
}
*/
 
/* blindly appends b2 after b1
 */
static Pbase append_to(Pbase b1, Pbase b2)
{
  Pbase b;
 
  if (BASE_NULLE_P(b1)) return b2;
 
  for (b=b1; !BASE_NULLE_P(b->succ); b=b->succ);
  b->succ = b2;
  return b1;
}
 
/* returns a newly allocated base with the scanners ahead
 */
static Pbase scanners_then_others(
  Pbase initial,           /* full base */
  list /* of entity */ ls) /* scanners */
{
    Pbase sb = BASE_NULLE, ob = BASE_NULLE;
    Pbase b;
 
    for (b=initial; !BASE_NULLE_P(b); b=b->succ) {
      Variable v = var_of(b);
      if (gen_in_list_p((entity) v, ls))
        base_add_dimension(&sb, v);
      else
        base_add_dimension(&ob, v);
    }
 
    return append_to(sb, ob);
}
 
/* void extract_lattice
 *
 * what: extracts a lattice from a set of equalities
 * how: Hermite form computation on equalities implying scanners
 *  - equalities translated in F and M;
 *  - variables = scanners (s) + others (o) + cst (1)
 *  - F.s + M.o + V == 0
 *  - H = PFQ; // I assert P==I, maybe optimistic for this implementation...
 *  - (1) s = Q.y
 *  - (2) H.y + M.o + V == 0
 *   then
 *  - (2) => new equalities that define some y(r)
 *  - (1) => ddc in some order..., plus replacement in inequalities
 *  - (1) => newscs, but what about the order?
 *
 * input: Psysteme and scanners
 * output: modified system, new scanners and deducibles
 * side effects:
 *  - may create some new variables
 * bugs or features:
 */
void
extract_lattice(
  Psysteme s,                      /* the system is modified */
  list /* of entity */ scanners,   /* variables to be scanned */
  list /* of entity */ *newscs,    /* returned new scanners */
  list /* of expression */ *ddc)   /* old deduction */
{
    /* - should try to remove deducibles before hand?
     */
    int nscanners, nothers, ntotal, neq;
    Pbase b, bsorted, byr;
    Pmatrix FM, F, M, V, P, H, Q, Hl, Hli, Ql, Qr, QlHli,
      QlHliM, QlHliV, mQr, I, Fnew;
    Value det_P, det_Q;
    int i;
    list /* of entity */ lns = NIL, ltmp = NIL;
    Pcontrainte eq;
 
    neq = sc_nbre_egalites(s);
 
    if (neq==0 || !get_bool_property("HPFC_EXTRACT_LATTICE")) {
      /* void implementation: nothing done!
       */
      *newscs = gen_copy_seq(scanners);
      *ddc = NIL;
      return;
    }
    /* else do the job */
 
    DEBUG_SYST(3, "initial system", s);
    DEBUG_ELST(3, "scanners", scanners);
 
    b = sc_base(s);
    nscanners = gen_length(scanners);
    ntotal = base_dimension(b);
    nothers = ntotal - nscanners;
 
    message_assert("more scanners than equalities", nscanners>=neq);
 
    bsorted = scanners_then_others(b, scanners);
 
    DEBUG_BASE(3, "sorted base", bsorted);
    pips_debug(3, "%d scanners, %d others, %d eqs\n", nscanners, nothers, neq);
 
    /* FM (so) + V == 0 */
    FM = matrix_new(neq, ntotal);
    V  = matrix_new(neq, 1);
 
    constraints_to_matrices(sc_egalites(s), bsorted, FM, V);
 
    DEBUG_MTRX(4, "FM", FM);
    DEBUG_MTRX(4, "V", V);
 
    /* Fs + Mo + V == 0 */
    F = matrix_new(neq, nscanners);
    M = matrix_new(neq, nothers);
 
    ordinary_sub_matrix(FM, F, 1, neq, 1, nscanners);
    ordinary_sub_matrix(FM, M, 1, neq, nscanners+1, ntotal);
 
    matrix_free(FM);
 
    DEBUG_MTRX(4, "F", F);
    DEBUG_MTRX(4, "M", M);
 
    /* H = P * F * Q */
    H = matrix_new(neq, nscanners);
    P = matrix_new(neq, neq);
    Q = matrix_new(nscanners, nscanners);
 
    matrix_hermite(F, P, H, Q, &det_P, &det_Q);
 
    DEBUG_MTRX(4, "H", H);
    DEBUG_MTRX(4, "P", P);
    DEBUG_MTRX(4, "Q", Q);
 
    message_assert("P == I", matrix_diagonal_p(P) && det_P==1);
 
    /* H = (Hl 0) */
    Hl = matrix_new(neq, neq);
    ordinary_sub_matrix(H, Hl, 1, neq, 1, neq);
    matrix_free(H);
 
    DEBUG_MTRX(4, "Hl", Hl);
 
    if (!matrix_triangular_unimodular_p(Hl, true)) {
      pips_user_warning("fast exit, some yes/no lattice skipped\n");
      /* and memory leak, by the way */
      *newscs = gen_copy_seq(scanners);
      *ddc = NIL;
      return;
    }
 
    message_assert("Hl is lower triangular unimodular",
                   matrix_triangular_unimodular_p(Hl, true));
 
    /* Hli = Hl^-1 */
    Hli = matrix_new(neq, neq);
    matrix_unimodular_triangular_inversion(Hl, Hli, true);
    matrix_free(Hl);
 
    DEBUG_MTRX(4, "Hli", Hli);
 
    /* Q = (Ql Qr) */
    Ql = matrix_new(nscanners, neq);
    Qr = matrix_new(nscanners, nscanners-neq);
 
    ordinary_sub_matrix(Q, Ql, 1, nscanners, 1, neq);
    ordinary_sub_matrix(Q, Qr, 1, nscanners, neq+1, nscanners);
 
    DEBUG_MTRX(4, "Ql", Ql);
    DEBUG_MTRX(4, "Qr", Qr);
 
    matrix_free(Q);
 
    /* QlHli = Ql * Hl^-1 */
    QlHli = matrix_new(nscanners, neq);
    matrix_multiply(Ql, Hli, QlHli);
 
    matrix_free(Ql);
    matrix_free(Hli);
 
    /* QlHliM = QlHli * M */
    QlHliM = matrix_new(nscanners, nothers);
    matrix_multiply(QlHli, M, QlHliM);
 
    matrix_free(M);
 
    /* QlHliV = QlHli * V */
    QlHliV = matrix_new(nscanners, 1);
    matrix_multiply(QlHli, V, QlHliV);
 
    matrix_free(V);
    matrix_free(QlHli);
 
    /* I */
    I = matrix_new(nscanners, nscanners);
    matrix_identity(I, 0);
 
    /* mQr = - Qr */
    mQr = matrix_new(nscanners, nscanners-neq);
    matrix_uminus(Qr, mQr);
    matrix_free(Qr);
 
    /* create nscanners-neq new scanning variables... they are the yr's.
     */
    for (i=0; i<nscanners-neq; i++)
      lns = CONS(ENTITY,
                 hpfc_new_variable(node_module, MakeBasic(is_basic_int)), lns);
 
    byr = list_to_base(lns);
    bsorted = append_to(byr, bsorted); byr = BASE_NULLE;
 
    /* We have: mQr yr + I s + QlHliM o + QlHliV == 0
     * yr are the new scanners, s the old ones, deducable from the new ones.
     * the equation must also be used to remove s from the inequalities.
     *
     * Fnew = ( mQr I QlHliM )
     */
    Fnew = matrix_new(nscanners, 2*nscanners-neq+nothers);
 
    insert_sub_matrix(Fnew, mQr, 1, nscanners, 1, nscanners-neq);
    insert_sub_matrix(Fnew, I, 1, nscanners, nscanners-neq+1, 2*nscanners-neq);
    insert_sub_matrix(Fnew, QlHliM, 1, nscanners,
                      2*nscanners-neq+1, 2*nscanners-neq+nothers);
 
    matrix_free(I);
    matrix_free(mQr);
    matrix_free(QlHliM);
 
    /* Now we have:
     *   (a) Fnew (yr s o)^t + QlHliV == 0
     *   (b) lns -- the new scanners
     *
     * we must
     *  (1) generate deducables from (a),
     *  (2) regenerate inequalities on yr's.
     */
 
    matrices_to_constraints(&eq, bsorted, Fnew, QlHliV);
    matrix_free(Fnew);
    matrix_free(QlHliV);
 
    /* clean the new system
     */
    contraintes_free(sc_egalites(s));
    sc_egalites(s) = eq;
    base_rm(sc_base(s));
    sc_creer_base(s);
 
    /* old scanners are deduced now:
     */
    *ddc = gen_append(*ddc, simplify_deducable_variables(s, scanners, &ltmp));
    pips_assert("no vars left", ENDP(ltmp));
 
    *newscs = lns;
 
    base_rm(bsorted);
 
    DEBUG_SYST(3, "resulting system", s);
    DEBUG_ELST(3, "new scanners", lns);
}
 
/* that is all
 */