Code/Sandbox/Plugins/LodGeneratorPlugin/Util/SparseMatrix.cpp

   1 #include "StdAfx.h"
   2 #include "SparseMatrix.h"
   3
   4 namespace LODGenerator {
   5
   6
   7     void CSparseLine::add(int index, double val) {
   8         CCoeff* coeff = NULL;
   9
  10         for(int ii=0; ii < nb_coeffs_; ii++) {
  11             if(coeff_[ii].index == index) {
  12                 coeff = &(coeff_[ii]);
  13                 break;
  14             }
  15         }
  16         if(coeff != NULL) {
  17             coeff->a += val;
  18         } else {
  19             nb_coeffs_++;
  20             if(nb_coeffs_ > capacity_) {
  21                 grow();
  22             }
  23             coeff = &(coeff_[nb_coeffs_ - 1]);
  24             coeff->a     = val;
  25             coeff->index = index;
  26         }
  27     }
  28
  29     void CSparseLine::set(int index, double val) {
  30         CCoeff* coeff = NULL;
  31
  32         for(int ii=0; ii < nb_coeffs_; ii++) {
  33             if(coeff_[ii].index == index) {
  34                 coeff = &(coeff_[ii]);
  35                 break;
  36             }
  37         }
  38         if(coeff != NULL) {
  39             coeff->a = val;
  40         } else {
  41             nb_coeffs_++;
  42             if(nb_coeffs_ > capacity_) {
  43                 grow();
  44             }
  45             coeff = &(coeff_[nb_coeffs_ - 1]);
  46             coeff->a     = val;
  47             coeff->index = index;
  48         }
  49     }
  50
  51         int CSparseLine::find(int i)
  52         {
  53
  54         for(int ii=0; ii < nb_coeffs_; ii++) {
  55             if(coeff_[ii].index == i) {
  56                 return ii;
  57             }
  58         }
  59                 return -1;
  60         }
  61
  62         int CSparseLine::findSorted(int i)
  63         {
  64                 if(i<coeff_[0].index)
  65                         return -1;
  66                 if(i>coeff_[nb_coeffs_-1].index)
  67                         return -1;
  68
  69                 int left=0, right=nb_coeffs_-1;
  70                 if(coeff_[left].index == i)
  71                         return left;
  72                 if(coeff_[right].index == i)
  73                         return right;
  74                 int middle = (left+right)/2;
  75
  76                 while (left!=middle)
  77                 {
  78                         if(coeff_[middle].index == i)
  79                                 return middle;
  80                         if(coeff_[middle].index < i)
  81                                 left = middle;
  82                         else
  83                                 right = middle;
  84                         middle = (left+right)/2;
  85                 }
  86                 return -1;
  87         }
  88
  89         void CSparseLine::push_element(int index, double val)
  90         {
  91                 nb_coeffs_++;
  92                 if(nb_coeffs_ > capacity_) {
  93                         grow();
  94                 }
  95                 CCoeff* coeff = &(coeff_[nb_coeffs_ - 1]);
  96                 coeff->a     = val;
  97                 coeff->index = index;
  98         }
  99
 100         void CSparseLine::reserve(int nb_elements)
 101         {
 102                 if (nb_elements<capacity_) return;
 103
 104                 int old_capacity = capacity_;
 105                 capacity_ = nb_elements;
 106         CAllocator<CCoeff>::reallocate(
 107             coeff_, old_capacity, capacity_
 108         );
 109         }
 110
 111         void CSparseLine::erase_element(int index)
 112         {
 113                 if (index!=nb_coeffs_-1)
 114                         memmove(&coeff_[index], &coeff_[index+1], (nb_coeffs_-index-1)*sizeof(coeff_[0]));
 115                 nb_coeffs_--;
 116         }
 117
 118     double CSparseLine::get(int index) const {
 119         CCoeff* coeff = NULL;
 120         // Search for a_{index}
 121         for(int ii=0; ii < nb_coeffs_; ii++) {
 122             if(coeff_[ii].index == index) {
 123                 coeff = &(coeff_[ii]);
 124                 break;
 125             }
 126         }
 127         if(coeff != NULL) {
 128             return coeff->a ;
 129         } else {
 130             return 0;
 131         }
 132     }
 133
 134     void CSparseLine::grow() {
 135         int old_capacity = capacity_;
 136                 if (capacity_>1000) {
 137                         capacity_ = (int)(capacity_ * 1.2);
 138                 } else {
 139                         capacity_ = capacity_ * 2;
 140                 }
 141         CAllocator<CCoeff>::reallocate(
 142             coeff_, old_capacity, capacity_
 143         );
 144     }
 145
 146     class CCoeffIndexCompare {
 147     public:
 148         bool operator()(const CCoeff& c1, const CCoeff& c2) {
 149             return c1.index < c2.index;
 150         }
 151     };
 152
 153     void CSparseLine::sort() {
 154         CCoeff* begin = coeff_;
 155         CCoeff* end   = coeff_ + nb_coeffs_;
 156         std::sort(begin, end, CCoeffIndexCompare());
 157     }
 158
 159     void CSparseLine::removeDuplicates() {
 160
 161                 if (nb_coeffs_==0) return;
 162
 163                 int nbDifferentCCoeffs = 1;
 164                 for (int i=1; i<nb_coeffs_; i++)
 165                         if (coeff_[i-1].index!=coeff_[i].index)
 166                                 nbDifferentCCoeffs++;
 167
 168                 CCoeff* newCCoeffs = CAllocator<CCoeff>::allocate(nbDifferentCCoeffs);
 169                 newCCoeffs[0] = coeff_[0];
 170                 int k=1;
 171                 for (int i=1; i<nb_coeffs_; i++)
 172                 {
 173                         if (coeff_[i-1].index!=coeff_[i].index)
 174                         {
 175                                 newCCoeffs[k] = coeff_[i];
 176                                 k++;
 177                         }
 178                         else
 179                                 newCCoeffs[k-1].a = min(coeff_[i].a, newCCoeffs[k-1].a);
 180                 }
 181                 CAllocator<CCoeff>::deallocate(coeff_);
 182                 coeff_ = newCCoeffs;
 183                 nb_coeffs_ = nbDifferentCCoeffs;
 184                 capacity_ = nbDifferentCCoeffs;
 185     }
 186
 187     CSparseMatrix::CSparseMatrix(int m, int n, EStorage storage) {
 188         storage_ = eS_NONE;
 189         allocate(m,n,storage,false);
 190     }
 191
 192     CSparseMatrix::CSparseMatrix(
 193         int n, EStorage storage, bool symmetric_storage
 194     ) {
 195         storage_ = eS_NONE;
 196         allocate(n,n,storage,symmetric_storage);
 197     }
 198
 199         CSparseMatrix::CSparseMatrix(const CSparseMatrix& rhs)
 200         {
 201                 storage_ = eS_NONE;
 202                 allocate(rhs.m(), rhs.n(), rhs.storage(), rhs.is_symmetric());
 203                 memcpy(diag_, rhs.diag_, diag_size_*sizeof(diag_[0]));
 204
 205                 if (rows_are_stored_)
 206                 {
 207                         for (int i=0; i<m_; i++)
 208                         {
 209                                 row_[i].coeff_ = new CCoeff[std::max(2,rhs.row_[i].nb_coeffs())];
 210                                 memcpy(row_[i].coeff_, rhs.row_[i].coeff_, rhs.row_[i].nb_coeffs()*sizeof(row_[i].coeff_[0]));
 211                                 row_[i].nb_coeffs_ = rhs.row_[i].nb_coeffs();
 212                                 row_[i].capacity_ = std::max(2, rhs.row_[i].nb_coeffs());
 213                         }
 214                 }
 215                 if (columns_are_stored_)
 216                 {
 217                         for (int i=0; i<n_; i++)
 218                         {
 219                                 column_[i].coeff_ = new CCoeff[std::max(2,rhs.column_[i].nb_coeffs())];
 220                                 memcpy(column_[i].coeff_, rhs.column_[i].coeff_, rhs.column_[i].nb_coeffs()*sizeof(column_[i].coeff_[0]));
 221                                 column_[i].nb_coeffs_ = rhs.column_[i].nb_coeffs();
 222                                 column_[i].capacity_ = std::max(2, rhs.column_[i].nb_coeffs());
 223                         }
 224                 }
 225         }
 226
 227         CSparseMatrix& CSparseMatrix::operator=(const CSparseMatrix& rhs)
 228         {
 229                 storage_ = eS_NONE;
 230                 allocate(rhs.m(), rhs.n(), rhs.storage(), rhs.is_symmetric());
 231                 memcpy(diag_, rhs.diag_, diag_size_*sizeof(diag_[0]));
 232
 233                 if (rows_are_stored_)
 234                 {
 235                         for (int i=0; i<m_; i++)
 236                         {
 237                                 row_[i].coeff_ = new CCoeff[std::max(2,rhs.row_[i].nb_coeffs())];
 238                                 memcpy(row_[i].coeff_, rhs.row_[i].coeff_, rhs.row_[i].nb_coeffs()*sizeof(row_[i].coeff_[0]));
 239                                 row_[i].nb_coeffs_ = rhs.row_[i].nb_coeffs();
 240                                 row_[i].capacity_ = std::max(2, rhs.row_[i].nb_coeffs());
 241                         }
 242                 }
 243                 if (columns_are_stored_)
 244                 {
 245                         for (int i=0; i<n_; i++)
 246                         {
 247                                 column_[i].coeff_ = new CCoeff[std::max(2,rhs.column_[i].nb_coeffs())];
 248                                 memcpy(column_[i].coeff_, rhs.column_[i].coeff_, rhs.column_[i].nb_coeffs()*sizeof(column_[i].coeff_[0]));
 249                                 column_[i].nb_coeffs_ = rhs.column_[i].nb_coeffs();
 250                                 column_[i].capacity_ = std::max(2, rhs.column_[i].nb_coeffs());
 251                         }
 252                 }
 253                 return *this;
 254         }
 255
 256     CSparseMatrix::~CSparseMatrix() {
 257         deallocate();
 258     }
 259
 260     CSparseMatrix::CSparseMatrix() {
 261         m_ = 0;
 262         n_ = 0;
 263         diag_size_ = 0;
 264
 265         row_ = NULL;
 266         column_ = NULL;
 267         diag_ = NULL;
 268
 269         storage_ = eS_NONE;
 270         rows_are_stored_ = false;
 271         columns_are_stored_ = false;
 272         symmetric_storage_ = false;
 273         symmetric_tag_ = false;
 274     }
 275
 276     int CSparseMatrix::nnz() const {
 277         int result = 0;
 278         if(rows_are_stored()) {
 279             for(int i=0; i<m(); i++) {
 280                 result += row(i).nb_coeffs();
 281             }
 282         } else if(columns_are_stored()) {
 283             for(int j=0; j<n(); j++) {
 284                 result += column(j).nb_coeffs();
 285             }
 286         } else {
 287             assert(false);
 288         }
 289         return result;
 290     }
 291
 292     void CSparseMatrix::deallocate() {
 293         m_ = 0;
 294         n_ = 0;
 295         diag_size_ = 0;
 296
 297         delete[] row_;
 298         delete[] column_;
 299         delete[] diag_;
 300         row_ = NULL;
 301         column_ = NULL;
 302         diag_ = NULL;
 303
 304         storage_ = eS_NONE;
 305         rows_are_stored_ = false;
 306         columns_are_stored_ = false;
 307         symmetric_storage_ = false;
 308     }
 309
 310     void CSparseMatrix::allocate(
 311         int m, int n, EStorage storage, bool symmetric_storage
 312     ) {
 313         if(storage_ != eS_NONE) {
 314             deallocate();
 315         }
 316
 317         m_ = m;
 318         n_ = n;
 319         diag_size_ = min(m,n);
 320         symmetric_storage_ = symmetric_storage;
 321         symmetric_tag_ = false;
 322         storage_ = storage;
 323         switch(storage) {
 324         case eS_NONE:
 325             assert(false);
 326             break;
 327         case eS_ROWS:
 328             rows_are_stored_    = true;
 329             columns_are_stored_ = false;
 330             break;
 331         case eS_COLUMNS:
 332             rows_are_stored_    = false;
 333             columns_are_stored_ = true;
 334             break;
 335         case eS_ROWS_AND_COLUMNS:
 336             rows_are_stored_    = true;
 337             columns_are_stored_ = true;
 338             break;
 339         }
 340         diag_ = new double[diag_size_];
 341         for(int i=0; i<diag_size_; i++) {
 342             diag_[i] = 0.0;
 343         }
 344
 345         if(rows_are_stored_) {
 346             row_ = new Row[m];
 347         } else {
 348             row_ = NULL;
 349         }
 350
 351         if(columns_are_stored_) {
 352             column_ = new Column[n];
 353         } else {
 354             column_ = NULL;
 355         }
 356     }
 357
 358     void CSparseMatrix::sort() {
 359         if(rows_are_stored_) {
 360             for(int i=0; i<m_; i++) {
 361                 row(i).sort();
 362             }
 363         }
 364         if(columns_are_stored_) {
 365             for(int j=0; j<n_; j++) {
 366                 column(j).sort();
 367             }
 368         }
 369     }
 370
 371     void CSparseMatrix::zero() {
 372         if(rows_are_stored_) {
 373             for(int i=0; i<m_; i++) {
 374                 row(i).zero();
 375             }
 376         }
 377         if(columns_are_stored_) {
 378             for(int j=0; j<n_; j++) {
 379                 column(j).zero();
 380             }
 381         }
 382         for(int i=0; i<diag_size_; i++) {
 383             diag_[i] = 0.0;
 384         }
 385     }
 386
 387     void CSparseMatrix::clear() {
 388         if(rows_are_stored_) {
 389             for(int i=0; i<m_; i++) {
 390                 row(i).clear();
 391             }
 392         }
 393         if(columns_are_stored_) {
 394             for(int j=0; j<n_; j++) {
 395                 column(j).clear();
 396             }
 397         }
 398         for(int i=0; i<diag_size_; i++) {
 399             diag_[i] = 0.0;
 400         }
 401     }
 402
 403         static void mult_rows_symmetric(
 404         const CSparseMatrix& A, const double* x, double* y
 405     ) {
 406         int m = A.m();
 407         for(int i=0; i<m; i++) {
 408             y[i] = 0;
 409             const CSparseMatrix::Row& Ri = A.row(i);
 410             for(int ij=0; ij<Ri.nb_coeffs(); ij++) {
 411                 const CCoeff& c = Ri.coeff(ij);
 412                 y[i]       += c.a * x[c.index];
 413                 if(i != c.index) {
 414                     y[c.index] += c.a * x[i];
 415                 }
 416             }
 417         }
 418     }
 419
 420     static void mult_rows(
 421         const CSparseMatrix& A, const double* x, double* y
 422     ) {
 423         int m = A.m();
 424         for(int i=0; i<m; i++) {
 425             y[i] = 0;
 426             const CSparseMatrix::Row& Ri = A.row(i);
 427             for(int ij=0; ij<Ri.nb_coeffs(); ij++) {
 428                 const CCoeff& c = Ri.coeff(ij);
 429                 y[i] += c.a * x[c.index];
 430             }
 431         }
 432     }
 433
 434
 435     static void mult_cols_symmetric(
 436         const CSparseMatrix& A, const double* x, double* y
 437     ) {
 438         ::memset(y, 0, A.m() * sizeof(double));
 439         int n = A.n();
 440         for(int j=0; j<n; j++) {
 441             const CSparseMatrix::Column& Cj = A.column(j);
 442             for(int ii=0; ii<Cj.nb_coeffs(); ii++) {
 443                 const CCoeff& c = Cj.coeff(ii);
 444                 y[c.index] += c.a * x[j]      ;
 445                 if(j != c.index) {
 446                     y[j]       += c.a * x[c.index];
 447                 }
 448             }
 449         }
 450     }
 451
 452     static void mult_cols(
 453         const CSparseMatrix& A, const double* x, double* y
 454     ) {
 455         ::memset(y, 0, A.m() * sizeof(double));
 456         int n = A.n();
 457         for(int j=0; j<n; j++) {
 458             const CSparseMatrix::Column& Cj = A.column(j);
 459             for(int ii=0; ii<Cj.nb_coeffs(); ii++) {
 460                 const CCoeff& c = Cj.coeff(ii);
 461                 y[c.index] += c.a * x[j];
 462             }
 463         }
 464     }
 465
 466     void mult(const CSparseMatrix& A, const double* x, double* y) {
 467         if(A.rows_are_stored()) {
 468             if(A.has_symmetric_storage()) {
 469                 mult_rows_symmetric(A, x, y);
 470             } else {
 471                 mult_rows(A, x, y);
 472             }
 473         } else {
 474             if(A.has_symmetric_storage()) {
 475                 mult_cols_symmetric(A, x, y);
 476             } else {
 477                 mult_cols(A, x, y);
 478             }
 479         }
 480     }
 481
 482     static void mult_xpose_rows(
 483         const CSparseMatrix& A, const double* x, double* y
 484     ) {
 485         ::memset(y, 0, A.n() * sizeof(double));
 486         int m = A.m();
 487         for(int i=0; i<m; i++) {
 488             const CSparseMatrix::Row& Ri = A.row(i);
 489             for(int ij=0; ij<Ri.nb_coeffs(); ij++) {
 490                 const CCoeff& c = Ri.coeff(ij);
 491                 y[c.index] += c.a * x[i];
 492             }
 493         }
 494     }
 495
 496     static void mult_xpose_cols(
 497         const CSparseMatrix& A, const double* x, double* y
 498     ) {
 499         int n = A.n();
 500         for(int j=0; j<n; j++) {
 501             y[j] = 0.0;
 502             const CSparseMatrix::Column& Cj = A.column(j);
 503             for(int ii=0; ii<Cj.nb_coeffs(); ii++) {
 504                 const CCoeff& c = Cj.coeff(ii);
 505                 y[j] += c.a * x[c.index];
 506             }
 507         }
 508     }
 509
 510     void mult_transpose(
 511         const CSparseMatrix& A, const double* x, double* y
 512     ) {
 513         if(A.rows_are_stored()) {
 514             if(A.has_symmetric_storage()) {
 515                 mult_rows_symmetric(A, x, y);
 516             } else {
 517                 mult_xpose_rows(A, x, y);
 518             }
 519         } else {
 520             if(A.has_symmetric_storage()) {
 521                 mult_cols_symmetric(A, x, y);
 522             } else {
 523                 mult_xpose_cols(A, x, y);
 524             }
 525         }
 526     }
 527
 528 }
 529
 530
 531