source/matrix.c

   1
   2    /**-------------------------------------------------------------------**
   3     **                               CLooG                               **
   4     **-------------------------------------------------------------------**
   5     **                             matrix.c                              **
   6     **-------------------------------------------------------------------**
   7     **                    First version: april 17th 2005                 **
   8     **-------------------------------------------------------------------**/
   9
  10
  11 /******************************************************************************
  12  *               CLooG : the Chunky Loop Generator (experimental)             *
  13  ******************************************************************************
  14  *                                                                            *
  15  * Copyright (C) 2005 Cedric Bastoul                                          *
  16  *                                                                            *
  17  * This is free software; you can redistribute it and/or modify it under the  *
  18  * terms of the GNU General Public License as published by the Free Software  *
  19  * Foundation; either version 2 of the License, or (at your option) any later *
  20  * version.                                                                   *
  21  *                                                                            *
  22  * This software is distributed in the hope that it will be useful, but       *
  23  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
  24  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
  25  * for more details.                                                          *
  26  *                                                                            *
  27  * You should have received a copy of the GNU General Public License along    *
  28  * with software; if not, write to the Free Software Foundation, Inc.,        *
  29  * 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA                     *
  30  *                                                                            *
  31  * CLooG, the Chunky Loop Generator                                           *
  32  * Written by Cedric Bastoul, Cedric.Bastoul@inria.fr                         *
  33  *                                                                            *
  34  ******************************************************************************/
  35 /* CAUTION: the english used for comments is probably the worst you ever read,
  36  *          please feel free to correct and improve it !
  37  */
  38
  39
  40 # include <stdlib.h>
  41 # include <stdio.h>
  42 # include <ctype.h>
  43 # include "../include/cloog/cloog.h"
  44
  45
  46
  47 /******************************************************************************
  48  *                             Memory leaks hunting                           *
  49  ******************************************************************************/
  50
  51
  52 /**
  53  * These functions and global variables are devoted to memory leaks hunting: we
  54  * want to know at each moment how many CloogMatrix structures are allocated
  55  * (cloog_matrix_allocated) and how many had been freed (cloog_matrix_freed).
  56  * Each time a CloogMatrix structure is allocated, a call to the function
  57  * cloog_matrix_leak_up() must be carried out, and respectively
  58  * cloog_matrix_leak_down() when a CloogMatrix structure is freed. The special
  59  * variable cloog_matrix_max gives the maximal number of CloogMatrix structures
  60  * simultaneously alive (i.e. allocated and non-freed) in memory.
  61  * - April 17th 2005: first version.
  62  */
  63
  64
  65 int cloog_matrix_allocated = 0 ;
  66 int cloog_matrix_freed = 0 ;
  67 int cloog_matrix_max = 0 ;
  68
  69
  70 void cloog_matrix_leak_up()
  71 { cloog_matrix_allocated ++ ;
  72   if ((cloog_matrix_allocated-cloog_matrix_freed) > cloog_matrix_max)
  73   cloog_matrix_max = cloog_matrix_allocated - cloog_matrix_freed ;
  74 }
  75
  76
  77 void cloog_matrix_leak_down()
  78 { cloog_matrix_freed ++ ;
  79 }
  80
  81
  82 /******************************************************************************
  83  *                              PolyLib interface                             *
  84  ******************************************************************************/
  85
  86
  87 /**
  88  * CLooG makes an intensive use of matrix operations and the PolyLib do
  89  * the job. Here are the interfaces to all the PolyLib calls (CLooG uses 18
  90  * PolyLib functions), with or without some adaptations. If another matrix
  91  * library can be used, only these functions have to be changed.
  92  */
  93
  94
  95 /**
  96  * cloog_matrix_print function:
  97  * This function prints the content of a CloogMatrix structure (matrix) into a
  98  * file (foo, possibly stdout).
  99  */
 100 void cloog_matrix_print(FILE * foo, CloogMatrix * matrix)
 101 { Matrix_Print(foo,P_VALUE_FMT,matrix) ;
 102 }
 103
 104
 105 /**
 106  * cloog_matrix_free function:
 107  * This function frees the allocated memory for a CloogMatrix structure
 108  * (matrix).
 109  */
 110 void cloog_matrix_free(CloogMatrix * matrix)
 111 { cloog_matrix_leak_down() ;
 112   Matrix_Free(matrix) ;
 113 }
 114
 115
 116 /**
 117  * cloog_matrix_alloc function:
 118  * This function allocates the memory space for a CloogMatrix structure having
 119  * nb_rows rows and nb_columns columns, it set its elements to 0.
 120  */
 121 CloogMatrix * cloog_matrix_alloc(unsigned nb_rows, unsigned nb_columns)
 122 { cloog_matrix_leak_up() ;
 123   return Matrix_Alloc(nb_rows,nb_columns) ;
 124 }
 125
 126
 127 /******************************************************************************
 128  *                          Structure display function                        *
 129  ******************************************************************************/
 130
 131
 132 /**
 133  * cloog_loop_print_structure function:
 134  * Displays a CloogMatrix structure (matrix) into a file (file, possibly stdout)
 135  * in a way that trends to be understandable without falling in a deep
 136  * depression or, for the lucky ones, getting a headache... It includes an
 137  * indentation level (level) in order to work with others print_structure
 138  * functions. Written by Olivier Chorier, Luc Marchaud, Pierre Martin and
 139  * Romain Tartiere.
 140  * - April 24th 2005: Initial version.
 141  * - June   2nd 2005: (Ced) Extraction from cloog_loop_print_structure and
 142  *                   integration in matrix.c.
 143  * - June  22rd 2005: Adaptation for GMP.
 144  */
 145 void cloog_matrix_print_structure(FILE * file, CloogMatrix * matrix, int level)
 146 { int i, j ;
 147
 148   /* Display the matrix. */
 149   for (i=0; i<matrix->NbRows; i++)
 150   { for(j=0; j<=level; j++)
 151     fprintf(file,"|\t") ;
 152
 153     fprintf(file,"[ ") ;
 154
 155     for (j=0; j<matrix->NbColumns; j++)
 156     { value_print(file,P_VALUE_FMT,matrix->p[i][j]) ;
 157       fprintf(file," ") ;
 158     }
 159
 160     fprintf(file,"]\n") ;
 161   }
 162 }
 163
 164
 165 /******************************************************************************
 166  *                               Reading function                             *
 167  ******************************************************************************/
 168
 169
 170 /**
 171  * cloog_matrix_read function:
 172  * Adaptation from the PolyLib. This function reads a matrix into a file (foo,
 173  * posibly stdin) and returns a pointer this matrix.
 174  * October 18th 2001: first version.
 175  * - April 17th 2005: this function moved from domain.c to here.
 176  * - June  21rd 2005: Adaptation for GMP (based on S. Verdoolaege's version of
 177  *                    CLooG 0.12.1).
 178  */
 179 CloogMatrix * cloog_matrix_read(FILE * foo)
 180 { unsigned NbRows, NbColumns ;
 181   int i, j, n ;
 182   char *c, s[MAX_STRING], str[MAX_STRING] ;
 183   CloogMatrix * matrix ;
 184   Value * p ;
 185
 186   while (fgets(s,MAX_STRING,foo) == 0) ;
 187   while ((*s=='#' || *s=='\n') || (sscanf(s," %u %u",&NbRows,&NbColumns)<2))
 188   fgets(s,MAX_STRING,foo) ;
 189
 190   matrix = cloog_matrix_alloc(NbRows,NbColumns) ;
 191
 192   p = matrix->p_Init ;
 193   for (i=0;i<matrix->NbRows;i++)
 194   { do
 195     { c = fgets(s,MAX_STRING,foo) ;
 196       while ((c != NULL) && isspace(*c) && (*c != '\n'))
 197       c++ ;
 198     }
 199     while (c != NULL && (*c == '#' || *c == '\n'));
 200
 201     if (c == NULL)
 202     { fprintf(stderr, "[CLooG]ERROR: not enough rows.\n") ;
 203       exit(1) ;
 204     }
 205     for (j=0;j<matrix->NbColumns;j++)
 206     { if (c == NULL || *c == '#' || *c == '\n')
 207       { fprintf(stderr, "[CLooG]ERROR: not enough columns.\n") ;
 208         exit(1) ;
 209       }
 210       /* NdCed : Dans le n ca met strlen(str). */
 211       if (sscanf(c,"%s%n",str,&n) == 0)
 212       { fprintf(stderr, "[CLooG]ERROR: not enough rows.\n") ;
 213         exit(1) ;
 214       }
 215       value_read(*(p++),str) ;
 216       c += n ;
 217     }
 218   }
 219
 220   return(matrix) ;
 221 }
 222
 223
 224 /******************************************************************************
 225  *                            Processing functions                            *
 226  ******************************************************************************/
 227
 228
 229 /**
 230  * Function cloog_matrix_normalize:
 231  * This function will modify the constraint system in such a way that when
 232  * there is an equality depending on the element at level 'level', there are
 233  * no more (in)equalities depending on this element. For instance, try
 234  * test/valilache.cloog with options -f 8 -l 9, with and without the call
 235  * to this function. At a given moment, for the level L we will have
 236  * 32*P=L && L>=1 (P is a lower level), this constraint system cannot be
 237  * translated directly into a source code. Thus, we normalize the domain to
 238  * remove L from the inequalities. In our example, this leads to
 239  * 32*P=L && 32*P>=1, that can be transated to the code
 240  * if (P>=1) { L=32*P ; ... }. This function solves the DaeGon Kim bug.
 241  * WARNING: Remember that if there is another call to Polylib after a call to
 242  * this function, we have to recall this function.
 243  *  -June    16th 2005: first version (adaptation from URGent June-7th-2005 by
 244  *                      N. Vasilache).
 245  * - June    21rd 2005: Adaptation for GMP.
 246  * - November 4th 2005: Complete rewriting, simpler and faster. It is no more an
 247  *                      adaptation from URGent.
 248  */
 249 void cloog_matrix_normalize(CloogMatrix * matrix, int level)
 250 { int ref, i, j ;
 251   Value factor_i, factor_ref, temp_i, temp_ref, gcd ;
 252
 253   if (matrix == NULL)
 254   return ;
 255
 256   /* Let us find an equality for the current level that can be propagated. */
 257   for (ref=0;ref<matrix->NbRows;ref++)
 258   if (value_zero_p(matrix->p[ref][0]) && value_notzero_p(matrix->p[ref][level]))
 259   { value_init_c(gcd) ;
 260     value_init_c(temp_i) ;
 261     value_init_c(temp_ref) ;
 262     value_init_c(factor_i) ;
 263     value_init_c(factor_ref) ;
 264
 265     /* Row "ref" is the reference equality, now let us find a row to simplify.*/
 266     for (i=ref+1;i<matrix->NbRows;i++)
 267     if (value_notzero_p(matrix->p[i][level]))
 268     { /* Now let us set to 0 the "level" coefficient of row "j" using "ref".
 269        * First we compute the factors to apply to each row vector element.
 270        */
 271       Gcd(matrix->p[ref][level],matrix->p[i][level],&gcd) ;
 272       value_division(factor_i,matrix->p[ref][level],gcd) ;
 273       value_division(factor_ref,matrix->p[i][level],gcd) ;
 274
 275       /* Maybe we are simplifying an inequality: factor_i must not be <0. */
 276       if (value_neg_p(factor_i))
 277       { value_absolute(factor_i,factor_i) ;
 278         value_oppose(factor_ref,factor_ref) ;
 279       }
 280
 281       /* Now update the vector. */
 282       for (j=1;j<matrix->NbColumns;j++)
 283       { value_multiply(temp_i,factor_i,matrix->p[i][j]) ;
 284         value_multiply(temp_ref,factor_ref,matrix->p[ref][j]) ;
 285         value_substract(matrix->p[i][j],temp_i,temp_ref) ;
 286       }
 287
 288       /* Normalize (divide by GCD of all elements) the updated vector. */
 289       Vector_Normalize(&(matrix->p[i][1]),matrix->NbColumns-1) ;
 290     }
 291
 292     value_clear_c(gcd) ;
 293     value_clear_c(temp_i) ;
 294     value_clear_c(temp_ref) ;
 295     value_clear_c(factor_i) ;
 296     value_clear_c(factor_ref) ;
 297     break ;
 298   }
 299 }
 300
 301
 302 /**
 303  * cloog_matrix_equality_update function:
 304  * this function updates a matrix of equalities where each row corresponds to
 305  * the equality "=0" of an affine expression such that the entry at column
 306  * "row" (="level") is not zero. This matrix is upper-triangular, except the
 307  * row number "level-1" which has to be updated for the matrix to be triangular.
 308  * The first element of each row gives the equality type and is not related to
 309  * the expression which is equal to zero. This function achieves the processing.
 310  * - equal is the matrix to be updated,
 311  * - level gives the row that has to be updated (it is actually row "level-1"),
 312  * - nb_par is the number of parameters of the program.
 313  **
 314  * - September 20th 2005: first version.
 315  */
 316 void cloog_matrix_equality_update(CloogMatrix * equal, int level, int nb_par)
 317 { int i, j ;
 318   Value gcd, factor_level, factor_outer, temp_level, temp_outer ;
 319
 320   value_init_c(gcd) ;
 321   value_init_c(temp_level) ;
 322   value_init_c(temp_outer) ;
 323   value_init_c(factor_level) ;
 324   value_init_c(factor_outer) ;
 325
 326   /* For each previous level, */
 327   for (i=level-2;i>=0;i--)
 328   { /* if the corresponding iterator is inside the current equality and is equal
 329      * to something,
 330      */
 331     if (value_notzero_p(equal->p[level-1][i+1]) &&
 332         value_notzero_p(equal->p[i][0]))
 333     { /* Compute the Greatest Common Divisor. */
 334       Gcd(equal->p[level-1][i+1],equal->p[i][i+1],&gcd) ;
 335
 336       /* Compute the factors to apply to each row vector element. */
 337       value_division(factor_level,equal->p[i][i+1],gcd) ;
 338       value_division(factor_outer,equal->p[level-1][i+1],gcd) ;
 339
 340       /* Now update the row 'level'. */
 341       /* - the iterators, up to level, */
 342       for (j=1;j<=level;j++)
 343       { value_multiply(temp_level,factor_level,equal->p[level-1][j]) ;
 344         value_multiply(temp_outer,factor_outer,equal->p[i][j]) ;
 345         value_substract(equal->p[level-1][j],temp_level,temp_outer) ;
 346       }
 347       /* - between last useful iterator (level) and the first parameter, the
 348        *   matrix is sparse (full of zeroes), we just do nothing there.
 349        * - the parameters and the scalar.
 350        */
 351       for (j=0;j<nb_par+1;j++)
 352       { value_multiply(temp_level,factor_level,
 353                        equal->p[level-1][equal->NbColumns-j-1]) ;
 354         value_multiply(temp_outer,factor_outer,
 355                        equal->p[i][equal->NbColumns-j-1]) ;
 356         value_substract(equal->p[level-1][equal->NbColumns-j-1],
 357                        temp_level,temp_outer) ;
 358       }
 359     }
 360   }
 361
 362   /* Normalize (divide by GCD of all elements) the updated equality. */
 363   Vector_Normalize(&(equal->p[level-1][1]),equal->NbColumns-1) ;
 364
 365   value_clear_c(gcd) ;
 366   value_clear_c(temp_level) ;
 367   value_clear_c(temp_outer) ;
 368   value_clear_c(factor_level) ;
 369   value_clear_c(factor_outer) ;
 370 }
 371
 372
 373 /**
 374  * cloog_matrix_copy function:
 375  * this functions builds and returns a "hard copy" (not a pointer copy) of a
 376  * CloogMatrix data structure.
 377  * - October 26th 2005: first version.
 378  */
 379 CloogMatrix * cloog_matrix_copy(CloogMatrix * matrix)
 380 { int i, j ;
 381   CloogMatrix * copy ;
 382
 383   copy = cloog_matrix_alloc(matrix->NbRows,matrix->NbColumns) ;
 384
 385   for (i=0;i<matrix->NbRows;i++)
 386   for (j=0;j<matrix->NbColumns;j++)
 387   value_assign(copy->p[i][j],matrix->p[i][j]) ;
 388
 389   return copy ;
 390 }
 391
 392
 393 /**
 394  * cloog_matrix_vector_copy function:
 395  * this function rutuns a hard copy of the vector "vector" of length "length"
 396  * provided as input.
 397  * - November 3rd 2005: first version.
 398  */
 399 Value * cloog_matrix_vector_copy(Value * vector, int length)
 400 { int i ;
 401   Value * copy ;
 402
 403   /* We allocate the memory space for the new vector, and we fill it with
 404    * the original coefficients.
 405    */
 406   copy = (Value *)malloc(length * sizeof(Value)) ;
 407   for (i=0;i<length;i++)
 408   { value_init_c(copy[i]) ;
 409     value_assign(copy[i],vector[i]) ;
 410   }
 411
 412   return copy ;
 413 }
 414
 415
 416 /**
 417  * cloog_matrix_vector_simplify function:
 418  * this function simplify an affine expression with its coefficients in
 419  * "vector" of length "length" thanks to an equality matrix "equal" that gives
 420  * for some elements of the affine expression an equality with other elements,
 421  * preferably constants. For instance, if the vector contains i+j+3 and the
 422  * equality matrix gives i=n and j=2, the vector is simplified to n+3 and is
 423  * returned in a new vector.
 424  * - vector is the array of affine expression coefficients
 425  * - equal is the matrix of equalities,
 426  * - length is the vector length,
 427  * - level is a level we don't want to simplify (-1 if none),
 428  * - nb_par is the number of parameters of the program.
 429  **
 430  * - September 20th 2005: first version.
 431  * - November   2nd 2005: (debug) we are simplifying inequalities, thus we are
 432  *                        not allowed to multiply the vector by a negative
 433  *                        constant.Problem found after a report of Michael
 434  *                        Classen.
 435  */
 436 Value * cloog_matrix_vector_simplify(vector, equal, length, level, nb_par)
 437 Value * vector ;
 438 CloogMatrix * equal ;
 439 int length, level, nb_par ;
 440 { int i, j ;
 441   Value gcd, factor_vector, factor_equal, temp_vector, temp_equal, * simplified;
 442
 443   simplified = cloog_matrix_vector_copy(vector,length) ;
 444
 445   value_init_c(gcd) ;
 446   value_init_c(temp_vector) ;
 447   value_init_c(temp_equal) ;
 448   value_init_c(factor_vector) ;
 449   value_init_c(factor_equal) ;
 450
 451   /* For each non-null coefficient in the vector, */
 452   for (i=length-nb_par-2;i>0;i--)
 453   if (i != level)
 454   { /* if the coefficient in not null, and there exists a useful equality */
 455     if ((value_notzero_p(simplified[i])) &&
 456         (value_notzero_p(equal->p[i-1][0])))
 457     { /* Compute the Greatest Common Divisor. */
 458       Gcd(simplified[i],equal->p[i-1][i],&gcd) ;
 459
 460       /* Compute the factors to apply to each row vector element. */
 461       value_division(factor_vector,equal->p[i-1][i],gcd) ;
 462       value_division(factor_equal,simplified[i],gcd) ;
 463
 464       /* We are simplifying an inequality: factor_vector must not be <0. */
 465       if (value_neg_p(factor_vector))
 466       { value_absolute(factor_vector,factor_vector) ;
 467         value_oppose(factor_equal,factor_equal) ;
 468       }
 469
 470       /* Now update the vector. */
 471       /* - the iterators, up to the current level, */
 472       for (j=1;j<=length-nb_par-2;j++)
 473       { value_multiply(temp_vector,factor_vector,simplified[j]) ;
 474         value_multiply(temp_equal,factor_equal,equal->p[i-1][j]) ;
 475         value_substract(simplified[j],temp_vector,temp_equal) ;
 476       }
 477       /* - between last useful iterator (i) and the first parameter, the equal
 478        *   matrix is sparse (full of zeroes), we just do nothing there.
 479        * - the parameters and the scalar.
 480        */
 481       for (j=0;j<nb_par+1;j++)
 482       { value_multiply(temp_vector,factor_vector,simplified[length-1-j]) ;
 483         value_multiply(temp_equal,factor_equal,
 484                        equal->p[i-1][equal->NbColumns-j-1]) ;
 485         value_substract(simplified[length-1-j],temp_vector,temp_equal) ;
 486       }
 487     }
 488   }
 489
 490   /* Normalize (divide by GCD of all elements) the updated vector. */
 491   Vector_Normalize(&(simplified[1]),length-1) ;
 492
 493   value_clear_c(gcd) ;
 494   value_clear_c(temp_vector) ;
 495   value_clear_c(temp_equal) ;
 496   value_clear_c(factor_vector) ;
 497   value_clear_c(factor_equal) ;
 498
 499   return simplified ;
 500 }
 501
 502
 503 /**
 504  * cloog_matrix_simplify function:
 505  * this function simplify all constraints inside the matrix "matrix" thanks to
 506  * an equality matrix "equal" that gives for some elements of the affine
 507  * constraint an equality with other elements, preferably constants.
 508  * For instance, if a row of the matrix contains i+j+3>=0 and the equality
 509  * matrix gives i=n and j=2, the constraint is simplified to n+3>=0. The
 510  * simplified constraints are returned back inside a new simplified matrix.
 511  * - matrix is the set of constraints to simplify,
 512  * - equal is the matrix of equalities,
 513  * - level is a level we don't want to simplify (-1 if none),
 514  * - nb_par is the number of parameters of the program.
 515  **
 516  * - November 4th 2005: first version.
 517  */
 518 CloogMatrix * cloog_matrix_simplify(matrix, equal, level, nb_par)
 519 CloogMatrix * matrix, * equal ;
 520 int level, nb_par ;
 521 { int i, j, k ;
 522   Value * vector ;
 523   CloogMatrix * simplified ;
 524
 525   if (matrix == NULL)
 526   return NULL ;
 527
 528   /* The simplified matrix is such that each row has been simplified thanks
 529    * tho the "equal" matrix. We allocate the memory for the simplified matrix,
 530    * then for each row of the original matrix, we compute the simplified
 531    * vector and we copy its content into the according simplified row.
 532    */
 533   simplified = cloog_matrix_alloc(matrix->NbRows,matrix->NbColumns) ;
 534   for (i=0;i<matrix->NbRows;i++)
 535   { vector = cloog_matrix_vector_simplify(matrix->p[i],equal,matrix->NbColumns,
 536                                           level,nb_par) ;
 537     for (j=0;j<matrix->NbColumns;j++)
 538     value_assign(simplified->p[i][j],vector[j]) ;
 539
 540     cloog_matrix_vector_free(vector,matrix->NbColumns) ;
 541   }
 542
 543   /* After simplification, it may happen that few constraints are the same,
 544    * we remove them here by replacing them with 0=0 constraints.
 545    */
 546   for (i=0;i<simplified->NbRows;i++)
 547   for (j=i+1;j<simplified->NbRows;j++)
 548   { for (k=0;k<simplified->NbColumns;k++)
 549     if (value_ne(simplified->p[i][k],simplified->p[j][k]))
 550     break ;
 551
 552     if (k == matrix->NbColumns)
 553     { for (k=0;k<matrix->NbColumns;k++)
 554       value_set_si(simplified->p[j][k],0) ;
 555     }
 556   }
 557
 558   return simplified ;
 559 }
 560
 561
 562 /**
 563  * cloog_matrix_vector_free function:
 564  * this function clears the elements of a vector "vector" of length "length",
 565  * then frees the vector itself this is useful for the GMP version of CLooG
 566  * which has to clear all elements.
 567  * - October 29th 2005: first version.
 568  */
 569 void cloog_matrix_vector_free(Value * vector, int length)
 570 { int i ;
 571
 572   for (i=0;i<length;i++)
 573   value_clear_c(vector[i]) ;
 574   free(vector) ;
 575 }