src/baseparsuif/suifmath/matrix.cc

   1 /* file "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 "matrix.h"
  13
  14 #define RCS_BASE_FILE matrix_cc
  15
  16 #include "matrix.h"
  17
  18 RCS_BASE("$Id: matrix.cc,v 1.1.1.1 1998/07/07 05:09:27 brm Exp $")
  19
  20 /*
  21  * Private member functions
  22  */
  23
  24 void matrix::init_decomp()
  25 {
  26     lu_decomp = fract_matrix();
  27     factored = FALSE;
  28     interchanges = NULL;
  29     pivot_cols = NULL;
  30 }
  31
  32 void matrix::clear_decomp()
  33 {
  34     lu_decomp.clear();
  35     factored = FALSE;
  36
  37     if (interchanges) {
  38         delete[] interchanges;
  39         interchanges = NULL;
  40     }
  41
  42     if (pivot_cols) {
  43         delete[] pivot_cols;
  44         pivot_cols = NULL;
  45     }
  46 }
  47
  48 void matrix::copy_decomp(const matrix &mat)
  49 {
  50     assert(factored && mat.factored);
  51
  52     lu_decomp = mat.lu_decomp;
  53
  54     int rows = lu_decomp.m();
  55     int cols = lu_decomp.n();
  56
  57     assert(mat.interchanges);
  58     interchanges = new int[rows];
  59     int i;
  60     for (i = 0; i < rows; i++)
  61         interchanges[i] = mat.interchanges[i];
  62
  63     assert(mat.pivot_cols);
  64     pivot_cols = new boolean[cols];
  65     for (i = 0; i < cols; i++)
  66         pivot_cols[i] = mat.pivot_cols[i];
  67 }
  68
  69 void matrix::copy_matrix(const matrix &mat)
  70 {
  71     elems = mat.elems;
  72
  73     if (mat.factored) {
  74         clear_decomp();
  75         factored = TRUE;
  76         copy_decomp(mat);
  77     }
  78     else {
  79         init_decomp();
  80     }
  81 }
  82
  83
  84 //
  85 // Linear Algebra support
  86 //
  87
  88 // Factors column 'j' in elems and places it in lu_decomp
  89 //
  90 void matrix::factor_col_and_insert(int j)
  91 {
  92     assert(j < elems.n());
  93
  94     int rows = elems.m();   // number of rows in resulting decomp
  95     int cols = j;           // column being factored = number of prev columns
  96
  97     fract_vector new_col(elems.col_list[j]);
  98
  99     int cur_pivot_row = 0;
 100     for (int i = 0; i < cols; i++)        // pivot row = number of pivot cols
 101         cur_pivot_row += pivot_cols[i];
 102
 103     assert(cur_pivot_row <= rows);
 104
 105     int new_pivot_row = factor_col(&new_col, cols, cur_pivot_row);
 106     assert((cur_pivot_row <= new_pivot_row && new_pivot_row < rows) ||
 107            (new_pivot_row == rows && cur_pivot_row == new_pivot_row));
 108
 109     if (cur_pivot_row < rows) {
 110         interchanges[cur_pivot_row] = new_pivot_row;
 111
 112         if (new_pivot_row != cur_pivot_row) {
 113
 114             // new_col has swapped rows, now swap prev cols
 115             for(int i = 0; i < cols; i++) {
 116               fract temp = lu_decomp.elt(new_pivot_row, i);
 117               lu_decomp.elt(new_pivot_row, i) = lu_decomp.elt(cur_pivot_row,i);
 118               lu_decomp.elt(cur_pivot_row, i) = temp;
 119             }
 120         }
 121
 122         pivot_cols[cols] = (new_col[cur_pivot_row] != fract(0));
 123     }
 124     else pivot_cols[cols] = FALSE;
 125
 126     lu_decomp.col_list[j] = new_col;
 127 }
 128
 129
 130 // Factor the given fract_vector column. Returns the row that the pivot row
 131 // is exchanged with.
 132 //
 133 int matrix::factor_col(fract_vector *vec, int cols, int pivot_row)
 134 {
 135     int rows = elems.m();
 136
 137     apply_L(vec, cols);   // Apply the L we have so far to the column vec
 138
 139     if (pivot_row == rows) return rows;
 140     assert(pivot_row < rows);
 141
 142     int new_pivot_row = pivot_row;
 143
 144     fract max = vec->elt(pivot_row);
 145
 146     if (max == fract(0)) {  // Find new pivot row to interchange
 147         for (int i = pivot_row+1; i < rows; i++) {
 148             max = vec->elt(i).abs();
 149
 150             if (max > fract(0)) {
 151                 new_pivot_row = i;
 152                 break;
 153             }
 154         }
 155     }
 156
 157     assert(pivot_row <= new_pivot_row && new_pivot_row < rows);
 158
 159     if (new_pivot_row != pivot_row) {  // swap
 160         fract temp = vec->elt(pivot_row);
 161         vec->elt(pivot_row) = vec->elt(new_pivot_row);
 162         vec->elt(new_pivot_row) = temp;
 163     }
 164
 165     if (vec->elt(pivot_row) != fract(0)) {
 166
 167         if (pivot_row == cols) {
 168             for (int i = pivot_row+1; i < rows; i++) {
 169                 vec->elt(i) /= vec->elt(pivot_row);
 170             }
 171         }
 172         else {
 173             assert(pivot_row < cols);
 174             for (int i = pivot_row+1; i < rows; i++) {
 175                 lu_decomp.elt(i,pivot_row) = vec->elt(i) / vec->elt(pivot_row);
 176                 vec->elt(i) = 0;
 177             }
 178         }
 179     }
 180
 181     return new_pivot_row;
 182 }
 183
 184 // Apply the L matrix to the column.
 185 //
 186 void matrix::apply_L(fract_vector *vec, int cols)
 187 {
 188     int rows = elems.m();
 189
 190     // Perform interchanges we've done so far
 191     //
 192     for (int i = 0; i < rows; i++) {
 193         if (interchanges[i] != i) {
 194             assert(interchanges[i] > i);
 195             fract temp = vec->elt(i);
 196             vec->elt(i) = vec->elt(interchanges[i]);
 197             vec->elt(interchanges[i]) = temp;
 198         }
 199     }
 200
 201     // Operations on column as dictated by L.
 202     //
 203     for (int j = 0; j < cols; j++) {
 204         for (int i = j+1; i < rows; i++) {
 205             vec->elt(i) -= vec->elt(j) * lu_decomp.elt(i,j);
 206         }
 207     }
 208 }
 209
 210
 211 // Backsolve the column 'vec' using U.
 212 // All free variables are assumed to be zero except the var specified by
 213 // 'free_var'.
 214 //
 215 fract_vector matrix::backsolve_U(fract_vector &vec, int free_var,
 216                                  boolean *valid)
 217 {
 218     assert(factored);
 219
 220     int rows = lu_decomp.m();
 221     int cols = lu_decomp.n();
 222
 223     fract_vector result(cols);
 224
 225     // First examine all zero rows:
 226     int last_row = last_nonzero_Urow();
 227
 228     for (int ii = last_row + 1; ii < rows; ii++) {
 229        if (vec[ii] != fract(0)) {
 230            *valid = FALSE;
 231            return result;  // All zeros
 232        }
 233     }
 234
 235     int i = last_row;
 236     for (int j = cols-1; j >= 0; j--) {
 237         if (!pivot_cols[j]) {
 238             if(j == free_var) result[j] = fract(1);
 239             else result[j] = fract(0);
 240         }
 241
 242         else {
 243             fract temp = vec[i];
 244
 245             for (int jj = j+1; jj < cols; jj++) {
 246                 temp -= lu_decomp.elt(i,jj) * result[jj];
 247             }
 248             result[j] = temp/lu_decomp.elt(i,j);
 249             i--;
 250         }
 251     }
 252
 253     *valid = TRUE;
 254     return result;
 255 }
 256
 257 int matrix::last_nonzero_Urow()
 258 {
 259     assert(factored);
 260
 261     int result = -1;
 262     for (int i = 0; i < lu_decomp.n(); i++)
 263         result += pivot_cols[i];
 264
 265     return result;
 266 }
 267
 268
 269 /*
 270  * Public member functions
 271  */
 272
 273 //
 274 // Constructors
 275 //
 276
 277 matrix::matrix()
 278 {
 279     elems = fract_matrix();
 280     init_decomp();
 281 }
 282
 283 matrix::matrix(int r, int c)
 284 {
 285     elems = fract_matrix(r, c);
 286     init_decomp();
 287 }
 288
 289 matrix::matrix(const integer_matrix &mat)
 290 {
 291     elems = fract_matrix(mat);
 292     init_decomp();
 293 }
 294
 295 matrix::matrix(const integer_matrix &mat, int div)
 296 {
 297     elems = fract_matrix(mat, div);
 298     init_decomp();
 299 }
 300
 301 matrix::matrix(const fract_matrix &mat)
 302 {
 303     elems = mat;
 304     init_decomp();
 305 }
 306
 307 matrix::matrix(const fract_matrix &mat, int c)
 308 {
 309     elems = fract_matrix(mat, c);
 310     init_decomp();
 311 }
 312
 313
 314 matrix::matrix(const matrix &mat)
 315 {
 316     init_decomp();
 317     copy_matrix(mat);
 318 }
 319
 320 matrix::matrix(const matrix &mat, int c)
 321 {
 322     elems = fract_matrix(mat.elems, c);
 323     init_decomp();
 324 }
 325
 326
 327 //
 328 // Equality and assignment operators:
 329 //
 330
 331 boolean matrix::operator==(const matrix &mat) const
 332 {
 333     return (elems == mat.elems);
 334 }
 335
 336 matrix &matrix::operator=(const matrix &mat)
 337 {
 338     copy_matrix(mat);
 339     return *this;
 340 }
 341
 342
 343 //
 344 // Matrix-Matrix operations:
 345 //
 346
 347 matrix matrix::operator*(matrix &mat)
 348 {
 349     return matrix(elems * mat.elems);
 350 }
 351
 352 matrix matrix::operator+(matrix &mat)
 353 {
 354     return matrix(elems + mat.elems);
 355 }
 356
 357 matrix matrix::operator-(matrix &mat)
 358 {
 359     return matrix(elems - mat.elems);
 360 }
 361
 362 matrix &matrix::operator+=(matrix &mat)
 363 {
 364     elems += mat.elems;
 365     factored = FALSE;
 366
 367     return *this;
 368 }
 369
 370 matrix &matrix::operator-=(matrix &mat)
 371 {
 372     elems -= mat.elems;
 373     factored = FALSE;
 374
 375     return *this;
 376 }
 377
 378
 379 //
 380 // Element-wise operations:
 381 //
 382 matrix matrix::operator*(fract val)
 383 {
 384     return matrix(elems * val);
 385 }
 386
 387 matrix matrix::operator/(fract val)
 388 {
 389     return matrix(elems / val);
 390 }
 391
 392 matrix matrix::operator+(fract val)
 393 {
 394     return matrix(elems + val);
 395 }
 396
 397 matrix matrix::operator-(fract val)
 398 {
 399     return matrix(elems - val);
 400 }
 401
 402 matrix &matrix::operator+=(fract val)
 403 {
 404     elems += val;
 405     factored = FALSE;
 406
 407     return *this;
 408 }
 409
 410 matrix &matrix::operator-=(fract val)
 411 {
 412     elems -= val;
 413     factored = FALSE;
 414
 415     return *this;
 416 }
 417
 418 matrix &matrix::operator*=(fract val)
 419 {
 420     elems *= val;
 421     factored = FALSE;
 422
 423     return *this;
 424 }
 425
 426 matrix &matrix::operator/=(fract val)
 427 {
 428     elems /= val;
 429     factored = FALSE;
 430
 431     return *this;
 432 }
 433
 434
 435 //
 436 // Matrix-vector
 437 //
 438 fract_vector matrix::operator*(fract_vector &vec)
 439 {
 440     return (elems * vec);
 441 }
 442
 443
 444 //
 445 // Other useful functions:
 446 //
 447 void matrix::ident(int n)
 448 {
 449     elems.ident(n);
 450     factored = FALSE;
 451 }
 452
 453 matrix matrix::transpose()
 454 {
 455     return matrix(elems.transpose());
 456 }
 457
 458
 459 //
 460 // Linear Algebra Operations
 461 //
 462
 463 boolean matrix::is_pivot(int i)
 464 {
 465     factor();
 466     return pivot_cols[i];
 467 }
 468
 469
 470 int matrix::rank()
 471 {
 472     int val = 0;
 473     int cols = elems.n();
 474
 475     factor();
 476     for (int i = 0; i < cols; i++)
 477         if (pivot_cols[i]) val++;
 478
 479     return (val);
 480 }
 481
 482
 483 void matrix::factor()
 484 {
 485     if (factored) return;
 486
 487     int rows = m();
 488     int cols = n();
 489
 490     clear_decomp();
 491     lu_decomp = fract_matrix(0, cols);  // Say rows are size 0 so no storage
 492     lu_decomp.rows = rows;              // is allocated.
 493
 494     interchanges = new int[rows];
 495     pivot_cols = new boolean[cols];
 496
 497     for (int r = 0; r < rows; r++)
 498         interchanges[r] = r;
 499
 500     for (int j = 0; j < cols; j++) {
 501         factor_col_and_insert(j);
 502     }
 503
 504     factored = TRUE;
 505 }
 506
 507 matrix matrix::inverse()
 508 {
 509     int rows = elems.m();
 510     int cols = elems.n();
 511
 512     assert_msg(rows == cols,
 513         ("Matrix of size %d x %d is not square", rows, cols));
 514
 515     factor();
 516
 517     int i;
 518     for (i = 0; i < m(); i++) {
 519         assert_msg(pivot_cols[i],
 520             ("Matrix is singular, col %d is not a pivot col", i));
 521     }
 522
 523     fract_matrix result(0, cols);         // Say rows are size 0 so no storage
 524     result.rows = rows;                   // is allocated.
 525
 526     for (i = 0; i < rows; i++) {
 527         fract_vector vec(rows);
 528         vec[i] = fract(1);
 529
 530         apply_L(&vec, cols);
 531         boolean valid;
 532         fract_vector solve_col = backsolve_U(vec, -1, &valid);
 533         assert(valid);
 534
 535         result.col_list[i] = solve_col;
 536     }
 537
 538     return matrix(result);
 539 }
 540
 541 vector_space matrix::kernel()
 542 {
 543     int rows = m();
 544     int cols = n();
 545
 546     factor();
 547
 548     vector_space result(cols);
 549     fract_vector zeros(rows);
 550
 551     for (int free_var = 0; free_var < cols; free_var++) {
 552         // the kernel has a dimension for each non-pivot column
 553         //
 554         if (!pivot_cols[free_var]) {
 555             boolean valid;
 556             fract_vector solve_col = backsolve_U(zeros, free_var, &valid);
 557             assert(valid);
 558             valid = result.insert(solve_col);
 559             assert(valid);
 560         }
 561     }
 562
 563     return result;
 564 }
 565
 566 fract_vector matrix::particular_solution(const fract_vector &vec,
 567                                          boolean *valid)
 568 {
 569     int cols = n();
 570     factor();
 571
 572     fract_vector new_vec(vec);
 573
 574     apply_L(&new_vec, cols);
 575     fract_vector result = backsolve_U(new_vec, -1, valid);
 576
 577     return result;
 578 }
 579
 580
 581 vector_space matrix::range()
 582 {
 583     return (vector_space(*this));
 584 }
 585
 586
 587 vector_space matrix::domain()
 588 {
 589     matrix trans = transpose();
 590     return (vector_space(trans));
 591 }
 592
 593
 594 matrix matrix::operator%(integer_row &rw)
 595 {
 596     return matrix(elems % rw);
 597 }
 598
 599 matrix matrix::del_row(int i, int j)
 600 {
 601     return matrix(elems.del_row(i,j));
 602 }
 603
 604
 605 matrix matrix::del_col(int i, int j)
 606 {
 607     return matrix(elems.del_col(i,j));
 608 }
 609
 610
 611 matrix matrix::del_columns(integer_row &rw)
 612 {
 613     return matrix(elems.del_columns(rw));
 614 }
 615
 616
 617 matrix matrix::insert_col(int i)
 618 {
 619     return matrix(elems.insert_col(i));
 620 }
 621
 622 matrix matrix::swap_col(int i, int j)
 623 {
 624     return matrix(elems.swap_col(i,j));
 625 }
 626
 627 matrix matrix::swap_row(int i, int j)
 628 {
 629     return matrix(elems.swap_row(i,j));
 630 }
 631
 632 fract_vector matrix::get_row(int i)
 633 {
 634     return elems.get_row(i);
 635 }
 636
 637 fract_vector matrix::get_col(int i)
 638 {
 639     return elems.get_col(i);
 640 }
 641
 642 void matrix::set_col(int i, fract_vector &vec)
 643 {
 644     elems.set_col(i, vec);
 645     factored = FALSE;
 646 }
 647
 648 void matrix::set_row(int i, fract_vector &vec)
 649 {
 650     elems.set_row(i, vec);
 651     factored = FALSE;
 652 }
 653
 654 matrix matrix::resize_offset(int r1, int r2, int c1, int c2, int fill)
 655 {
 656     return matrix(elems.resize_offset(r1, r2, c1, c2, fill));
 657 }
 658
 659 matrix matrix::r_merge(matrix &mat)
 660 {
 661     return matrix(elems.r_merge(mat.elems));
 662 }
 663
 664
 665 matrix matrix::c_merge(matrix &mat)
 666 {
 667     return matrix(elems.c_merge(mat.elems));
 668 }
 669
 670 immed_list *matrix::cvt_immed_list()
 671 {
 672     return elems.cvt_immed_list();
 673 }
 674
 675 matrix cvt_matrix(immed_list *il)
 676 {
 677     return matrix(cvt_fract_matrix(il));
 678 }
 679
 680 void matrix::print(FILE *fp)
 681 {
 682     elems.print(fp);
 683 }
 684
 685 void matrix::print_full(FILE *fp)
 686 {
 687     fprintf(fp, "Elements:\n");
 688     elems.print(fp);
 689
 690     if (factored) {
 691         fprintf(fp, "LU decomp:\n");
 692         lu_decomp.print(fp);
 693
 694         assert(interchanges);
 695         fprintf(fp, "Interchanges [");
 696         int i;
 697         for (i = 0; i < lu_decomp.m(); i++)
 698             fprintf(fp, "%s%d", (i == 0 ? "" : " "), interchanges[i]);
 699         fprintf(fp, "]\n");
 700
 701         assert(pivot_cols);
 702         fprintf(fp, "Pivot Columns [");
 703         for (i = 0; i < lu_decomp.n(); i++)
 704             fprintf(fp, "%s%d", (i == 0 ? "" : " "), pivot_cols[i]);
 705         fprintf(fp, "]\n");
 706     }
 707
 708     else {
 709         fprintf(fp, "Not factored\n");
 710     }
 711 }
 712
 713 void matrix::clear()
 714 {
 715     elems.clear();
 716     clear_decomp();
 717 }