src/baseparsuif/suifmath/fract_matrix.cc

   1 /* file "fract_matrix.cc" */
   2
   3 /*  Copyright (c) 1994,95 Stanford University
   4
   5     All rights reserved.
   6
   7     This software is provided under the terms described in
   8     the "suif_copyright.h" include file. */
   9
  10 #include <suif_copyright.h>
  11
  12 #pragma implementation "fract_matrix.h"
  13
  14 #define RCS_BASE_FILE fract_matrix_cc
  15
  16 #include "fract_matrix.h"
  17 #include <string.h>
  18
  19 RCS_BASE("$Id: fract_matrix.cc,v 1.1.1.1 1998/07/07 05:09:27 brm Exp $")
  20
  21 /*
  22  * Private member functions
  23  */
  24
  25 void fract_matrix::create_columns()
  26 {
  27     col_list = new fract_vector[cols];
  28
  29     for (int i = 0; i < cols; i++) {
  30         col_list[i] = fract_vector(rows);
  31     }
  32 }
  33
  34
  35 void fract_matrix::copy_fract_matrix(const fract_matrix &mat)
  36 {
  37     assert(rows == mat.m() && cols == mat.n());
  38
  39     for (int j = 0; j < cols; j++) {
  40         col_list[j] = mat.col_list[j];
  41     }
  42 }
  43
  44 void fract_matrix::copy_int_matrix(const integer_matrix &mat)
  45 {
  46     assert(rows == mat.m() && cols == mat.n());
  47
  48     for (int j = 0; j < cols; j++)
  49         for (int i = 0; i < rows; i++) {
  50             elt(i,j) = mat.r(i).c(j);
  51         }
  52 }
  53
  54
  55 /*
  56  * Public member functions
  57  */
  58
  59 //
  60 // Constructors
  61 //
  62
  63 fract_matrix::fract_matrix()
  64 {
  65     rows = cols = 0;
  66     col_list = NULL;
  67 }
  68
  69
  70 fract_matrix::fract_matrix(int r, int c)
  71 {
  72     assert(r >= 0 && c >= 0);
  73
  74     rows = r;
  75     cols = c;
  76     create_columns();
  77 }
  78
  79 fract_matrix::fract_matrix(const integer_matrix &mat)
  80 {
  81     assert(mat.m() >= 0 && mat.n() >= 0);
  82
  83     rows = mat.m();
  84     cols = mat.n();
  85     create_columns();
  86     copy_int_matrix(mat);
  87 }
  88
  89 fract_matrix::fract_matrix(const integer_matrix &mat, int div)
  90 {
  91     assert(mat.m() >= 0 && mat.n() >= 0);
  92
  93     rows = mat.m();
  94     cols = mat.n();
  95     create_columns();
  96     copy_int_matrix(mat);
  97
  98     *this /= fract(div);
  99 }
 100
 101 fract_matrix::fract_matrix(const fract_matrix &mat)
 102 {
 103     rows = mat.rows;
 104     cols = mat.cols;
 105     create_columns();
 106     copy_fract_matrix(mat);
 107 }
 108
 109 fract_matrix::fract_matrix(const fract_matrix &mat, int c)
 110 {
 111     assert(c >= 0);
 112
 113     rows = mat.rows;
 114     cols = c;
 115
 116     col_list = new fract_vector[cols];
 117
 118     if (cols <= mat.cols) {
 119         for (int i = 0; i < cols; i++) col_list[i] = mat.col_list[i];
 120     }
 121     else {
 122         int i;
 123         for (i = 0; i< mat.cols; i++) col_list[i] = mat.col_list[i];
 124         for (i = mat.cols; i< cols; i++) col_list[i] = fract_vector(rows);
 125     }
 126 }
 127
 128 //
 129 // Equality and assignment operators:
 130 //
 131
 132 boolean fract_matrix::operator==(const fract_matrix &mat) const
 133 {
 134     if ((rows != mat.rows) || (cols != mat.cols)) return FALSE;
 135
 136     for (int i = 0; i < cols; i++) {
 137         if (mat.col_list[i] != col_list[i]) return FALSE;
 138     }
 139
 140     return TRUE;
 141 }
 142
 143 fract_matrix & fract_matrix::operator=(const fract_matrix &mat)
 144 {
 145     if (&mat == this) return *this;
 146
 147     if ((rows != mat.rows) || (cols != mat.cols)) {
 148         clear();
 149         rows = mat.rows;
 150         cols = mat.cols;
 151         create_columns();
 152     }
 153
 154     copy_fract_matrix(mat);
 155     return *this;
 156 }
 157
 158
 159 //
 160 // Matrix-Matrix operations:
 161 //
 162
 163 fract_matrix fract_matrix::operator*(const fract_matrix &mat) const
 164 {
 165     assert(cols == mat.rows);
 166
 167     fract_matrix result(rows, mat.cols);
 168
 169     for (int i = 0; i < rows; i++) {
 170         for (int j = 0; j < mat.cols; j++) {
 171             fract temp = fract(0);
 172             for (int k = 0; k < cols; k++) {
 173                 temp += rc(i,k) * mat.rc(k,j);
 174             }
 175             result.elt(i,j) = temp;
 176         }
 177     }
 178
 179     return result;
 180 }
 181
 182 fract_matrix fract_matrix::operator+(const fract_matrix &mat) const
 183 {
 184     fract_matrix result = *this;
 185     result += mat;
 186
 187     return result;
 188 }
 189
 190 fract_matrix fract_matrix::operator-(const fract_matrix &mat) const
 191 {
 192     fract_matrix result = *this;
 193     result -= mat;
 194
 195     return result;
 196 }
 197
 198 fract_matrix & fract_matrix::operator+=(const fract_matrix &mat)
 199 {
 200     assert(rows == mat.rows);
 201     assert(cols == mat.cols);
 202
 203     for (int j = 0; j < cols; j++)
 204         for (int i = 0; i < rows; i++)
 205             elt(i,j) += mat.rc(i,j);
 206
 207     return *this;
 208 }
 209
 210 fract_matrix & fract_matrix::operator-=(const fract_matrix &mat)
 211 {
 212     assert(rows == mat.rows);
 213     assert(cols == mat.cols);
 214
 215     for (int j = 0; j < cols; j++)
 216         for (int i = 0; i < rows; i++)
 217             elt(i,j) -= mat.rc(i,j);
 218
 219     return *this;
 220 }
 221
 222 //
 223 // Element-wise operations:
 224 //
 225
 226 fract_matrix fract_matrix::operator*(fract val) const
 227 {
 228     fract_matrix result = *this;
 229
 230     result *= val;
 231     return result;
 232 }
 233
 234 fract_matrix fract_matrix::operator/(fract val) const
 235 {
 236     fract_matrix result = *this;
 237
 238     result /= val;
 239     return result;
 240 }
 241
 242 fract_matrix fract_matrix::operator+(fract val) const
 243 {
 244     fract_matrix result = *this;
 245
 246     result += val;
 247     return result;
 248 }
 249
 250 fract_matrix fract_matrix::operator-(fract val) const
 251 {
 252     fract_matrix result = *this;
 253
 254     result -= val;
 255     return result;
 256 }
 257
 258 fract_matrix &fract_matrix::operator+=(fract val)
 259 {
 260     for (int j = 0; j < cols; j++)
 261         for (int i = 0; i < rows; i++)
 262             elt(i,j) += val;
 263
 264     return *this;
 265 }
 266
 267 fract_matrix &fract_matrix::operator-=(fract val)
 268 {
 269     for (int j = 0; j < cols; j++)
 270         for (int i = 0; i < rows; i++)
 271             elt(i,j) -= val;
 272
 273     return *this;
 274 }
 275
 276 fract_matrix &fract_matrix::operator*=(fract val)
 277 {
 278     for (int j = 0; j < cols; j++)
 279         for (int i = 0; i < rows; i++)
 280             elt(i,j) *= val;
 281
 282     return *this;
 283 }
 284
 285 fract_matrix &fract_matrix::operator/=(fract val)
 286 {
 287     for (int j = 0; j < cols; j++)
 288         for (int i = 0; i < rows; i++)
 289             elt(i,j) /= val;
 290
 291     return *this;
 292 }
 293
 294 //
 295 // Matrix-vector
 296 //
 297
 298 fract_vector fract_matrix::operator*(const fract_vector &vec) const
 299 {
 300     fract_vector result(rows);
 301     assert(cols == vec.n());
 302
 303     for (int i = 0; i < rows; i++) {
 304         fract f(0);
 305         for (int j = 0; j < cols; j++) {
 306             f += rc(i,j) * vec.c(j);
 307         }
 308         result[i] = f;
 309     }
 310
 311     return result;
 312 }
 313
 314 //
 315 // Other useful functions:
 316 //
 317 void fract_matrix::ident(int n)
 318 {
 319     if ((rows != n) || (cols != n)) {
 320         clear();
 321         rows = n;
 322         cols = n;
 323         create_columns();
 324     }
 325
 326     for (int j = 0; j < cols; j++) {
 327         elt(j,j) = fract(1);
 328     }
 329 }
 330
 331 fract_matrix fract_matrix::transpose() const
 332 {
 333     fract_matrix result(cols, rows);
 334
 335     for (int j = 0; j < cols; j++)
 336         for (int i = 0; i < rows; i++)
 337             result.elt(j,i) = rc(i,j);
 338
 339     return result;
 340 }
 341
 342 fract_matrix fract_matrix::operator%(const integer_row &rw) const
 343 {
 344     fract_matrix result(rows, cols);
 345
 346     for (int j = 0; j < cols; j++) {
 347         result.col_list[j] = col_list[j] % rw;
 348     }
 349
 350     return result;
 351 }
 352
 353 fract_matrix fract_matrix::del_row(int i, int j) const
 354 {
 355     assert (i <= j);
 356     assert((i >= 0) && (j < rows));
 357
 358     fract_matrix result(rows - (j - i + 1), cols);
 359
 360     for (int k = 0; k < cols; k++)
 361         result.col_list[k] = col_list[k].del_pos(i, j);
 362
 363     return result;
 364 }
 365
 366
 367 fract_matrix fract_matrix::del_col(int i, int j) const
 368 {
 369     assert (i <= j);
 370     assert((i >= 0) && (j < cols));
 371
 372     fract_matrix result(rows, cols - (j - i + 1));
 373
 374     int count = 0;
 375     int k;
 376     for (k = 0; k < i; k++)
 377         result.col_list[count++] = col_list[k];
 378
 379     for (k = j+1; k < cols; k++)
 380         result.col_list[count++] = col_list[k];
 381
 382     assert(count == result.cols);
 383     return result;
 384 }
 385
 386
 387 fract_matrix fract_matrix::del_columns(const integer_row &rw) const
 388 {
 389     assert(rw.n() == cols);
 390
 391     int col_count = 0;
 392
 393     for (int i = 0; i < cols; i++) {
 394         if (rw.c(i)) col_count++;
 395     }
 396
 397     fract_matrix result(rows, col_count);
 398     int count = 0;
 399
 400     for (int j = 0; j < cols; j++)
 401         if (rw.c(j)) result.col_list[count++] = col_list[j];
 402
 403     assert(count == result.cols);
 404     return result;
 405 }
 406
 407 fract_matrix fract_matrix::insert_col(int i) const
 408 {
 409     assert(i >= 0 && i <= cols);
 410     fract_matrix result(rows, cols+1);
 411
 412     int count = 0;
 413     int j;
 414     for (j = 0; j < i; j++)
 415         result.col_list[count++] = col_list[j];
 416
 417     result.col_list[count++] = fract_vector(rows);
 418
 419     for (j = i+1; j < cols; j++)
 420         result.col_list[count++] = col_list[j];
 421
 422     assert(count == result.cols);
 423     return result;
 424 }
 425
 426 fract_matrix fract_matrix::swap_col(int i, int j) const
 427 {
 428     assert(i >= 0 && i < cols);
 429     assert(j >= 0 && j < cols);
 430     fract_matrix result = *this;
 431
 432     result.col_list[i] = col_list[j];
 433     result.col_list[j] = col_list[i];
 434
 435     return result;
 436 }
 437
 438 fract_matrix fract_matrix::swap_row(int i, int j) const
 439 {
 440     assert(i >= 0 && i < rows);
 441     assert(j >= 0 && j < rows);
 442     fract_matrix result(rows, cols);
 443
 444     for(int k = 0; k < cols; k++)
 445         result.col_list[j] = col_list[j].swap_pos(i, j);
 446
 447     return result;
 448 }
 449
 450
 451 fract_vector fract_matrix::get_row(int i) const
 452 {
 453     assert(i >= 0 && i < rows);
 454     fract_vector the_row(cols);
 455
 456     for (int j = 0; j < cols; j++) {
 457         the_row[j] = rc(i,j);
 458     }
 459
 460     return the_row;
 461 }
 462
 463
 464 fract_vector fract_matrix::get_col(int i) const
 465 {
 466     assert(i >= 0 && i < cols);
 467     return col_list[i];
 468 }
 469
 470
 471 void fract_matrix::set_col(int i, const fract_vector &vec)
 472 {
 473     assert(i >= 0 && i < cols);
 474     assert(vec.n() == rows);
 475
 476     col_list[i] = vec;
 477 }
 478
 479
 480 void fract_matrix::set_row(int i, const fract_vector &vec)
 481 {
 482     assert(i >= 0 && i < rows);
 483     assert(vec.n() == cols);
 484
 485     for (int j = 0; j < cols; j++) {
 486         elt(i,j) = vec.c(j);
 487     }
 488 }
 489
 490
 491 fract_matrix fract_matrix::resize_offset(int r1, int r2, int c1, int c2,
 492                                          int fill) const
 493 {
 494     assert(rows-r1+r2 >=0);
 495     assert(cols-c1+c2 >=0);
 496
 497     fract_matrix result(rows-r1+r2, cols-c1+c2);
 498
 499     int j;
 500     for (j = 0; j < result.cols; j++) result.col_list[j] = fill;
 501
 502     for (j = MAX(0,c1); j < cols+MIN(0,c2); j++)
 503         for (int i = MAX(0,r1); i < rows+MIN(0,r2); i++)
 504             result.elt(i-r1,j-c1) = rc(i,j);
 505
 506     return result;
 507 }
 508
 509 fract_matrix fract_matrix::r_merge(const fract_matrix &mat) const
 510 {
 511     assert(cols == mat.cols);
 512
 513     return (resize_offset(0, mat.rows, 0, 0) +
 514             mat.resize_offset(-rows, 0, 0, 0));
 515 }
 516
 517
 518 fract_matrix fract_matrix::c_merge(const fract_matrix &mat) const
 519 {
 520     assert(rows == mat.rows);
 521
 522     return (resize_offset(0, 0, 0, mat.cols) +
 523             mat.resize_offset(0, 0, -cols, 0));
 524 }
 525
 526 immed_list *fract_matrix::cvt_immed_list() const
 527 {
 528     immed_list *il = new immed_list;
 529
 530     il->append(immed("rows"));
 531     il->append(immed(rows));
 532     il->append(immed("cols"));
 533     il->append(immed(cols));
 534
 535     // Find lcm of denoms
 536     int curr_lcm = 1;
 537     for (int j = 0; j < cols; j++) {
 538         for (int i = 0; i < rows; i++) {
 539             curr_lcm = lcm(curr_lcm, rc(i,j).denom());
 540         }
 541     }
 542
 543     if (curr_lcm != 1) {
 544         il->append(immed("divisor"));
 545         il->append(immed(curr_lcm));
 546     }
 547
 548     il->append(immed("matrix"));
 549     for (int i = 0; i < rows; i++) {
 550         for (int j = 0; j < cols; j++) {
 551             fract temp = rc(i,j);
 552             temp *= curr_lcm;         // fract "*=" operator reduces fraction
 553             assert(temp.denom() == 1);
 554             il->append(immed(temp.num()));
 555         }
 556     }
 557
 558     return il;
 559 }
 560
 561 fract_matrix cvt_fract_matrix(immed_list *il)
 562 {
 563     assert(il->count() >= 5);  // at least: "rows" #rows "cols" #cols "matrix"
 564
 565     immed temp = il->pop();
 566     assert(temp.is_string() && !strcmp(temp.string(), "rows"));
 567     immed r = il->pop();
 568
 569     temp = il->pop();
 570     assert(temp.is_string() && !strcmp(temp.string(), "cols"));
 571     immed c = il->pop();
 572
 573     assert(r.is_integer() && c.is_integer());
 574
 575     fract_matrix result(r.integer(), c.integer());
 576     temp = il->pop();
 577     assert(temp.is_string());
 578
 579     // Grab divisor, if any
 580     //
 581     int divisor = 1;
 582     if (!strcmp(temp.string(), "divisor")) {
 583         assert(il->count() >= 2);
 584         immed div = il->pop();
 585         assert(div.is_integer());
 586         divisor = div.integer();
 587         temp = il->pop();
 588         assert(temp.is_string());
 589     }
 590
 591     assert(!strcmp(temp.string(), "matrix"));
 592     assert(il->count() >= result.rows * result.cols);
 593
 594     for (int i = 0; i < result.rows; i++) {
 595         for (int j = 0; j < result.cols; j++) {
 596             immed val = il->pop();
 597             assert(val.is_integer());
 598             result.elt(i,j) = fract(val.integer(), divisor);
 599         }
 600     }
 601
 602     return result;
 603 }
 604
 605 void fract_matrix::print(FILE *fp) const
 606 {
 607     for (int i = 0; i < rows; i++) {
 608         for (int j = 0; j < cols; j++) {
 609             rc(i,j).print(fp);
 610             fprintf(fp, " ");
 611         }
 612         fprintf(fp, "\n");
 613     }
 614 }
 615
 616 void fract_matrix::clear()
 617 {
 618     if (col_list) {
 619         delete[] col_list;
 620         col_list = NULL;
 621     }
 622
 623     rows = cols = 0;
 624 }