clebsch.cpp

   1 #include "clebsch.h"
   2
   3 // implementation of "binomial_t" starts here
   4
   5 clebsch::binomial_t binomial;
   6
   7 int clebsch::binomial_t::operator()(int n, int k) {
   8     if (N <= n) {
   9         for (cache.resize((n + 1) * (n + 2) / 2); N <= n; ++N) {
  10             cache[N * (N + 1) / 2] = cache[N * (N + 1) / 2 + N] = 1;
  11             for (int k = 1; k < N; ++k) {
  12                 cache[N * (N + 1) / 2 + k] = cache[(N - 1) * N / 2 + k]
  13                                            + cache[(N - 1) * N / 2 + k - 1];
  14             }
  15         }
  16     }
  17
  18     return cache[n * (n + 1) / 2 + k];
  19 }
  20
  21 // implementation of "weight" starts here
  22
  23 clebsch::weight::weight(int N) : elem(N), N(N) {}
  24
  25 clebsch::weight::weight(int N, int index) : elem(N, 0), N(N) {
  26     for (int i = 0; index > 0 && i < N; ++i) {
  27         for (int j = 1; binomial(N - i - 1 + j, N - i - 1) <= index; j <<= 1) {
  28             elem[i] = j;
  29         }
  30
  31         for (int j = elem[i] >> 1; j > 0; j >>= 1) {
  32             if (binomial(N - i - 1 + (elem[i] | j), N - i - 1) <= index) {
  33                 elem[i] |= j;
  34             }
  35         }
  36
  37         index -= binomial(N - i - 1 + elem[i]++, N - i - 1);
  38     }
  39 }
  40
  41 clebsch::weight &clebsch::weight::operator=(const clebsch::weight &w) {
  42     int &n = const_cast<int &>(N);
  43     elem = w.elem;
  44     n = w.N;
  45     return *this;
  46 }
  47
  48 int &clebsch::weight::operator()(int k) {
  49     assert(1 <= k && k <= N);
  50     return elem[k - 1];
  51 }
  52
  53 const int &clebsch::weight::operator()(int k) const {
  54     assert(1 <= k && k <= N);
  55     return elem[k - 1];
  56 }
  57
  58 bool clebsch::weight::operator<(const weight &w) const {
  59     assert(w.N == N);
  60     for (int i = 0; i < N; ++i) {
  61         if (elem[i] - elem[N - 1] != w.elem[i] - w.elem[N - 1]) {
  62             return elem[i] - elem[N - 1] < w.elem[i] - w.elem[N - 1];
  63         }
  64     }
  65     return false;
  66 }
  67
  68 bool clebsch::weight::operator==(const weight &w) const {
  69     assert(w.N == N);
  70
  71     for (int i = 1; i < N; ++i) {
  72         if (w.elem[i] - w.elem[i - 1] != elem[i] - elem[i - 1]) {
  73             return false;
  74         }
  75     }
  76
  77     return true;
  78 }
  79
  80 clebsch::weight clebsch::weight::operator+(const weight &w) const {
  81     weight result(N);
  82
  83     transform(elem.begin(), elem.end(), w.elem.begin(), result.elem.begin(), std::plus<int>());
  84
  85     return result;
  86 }
  87
  88 int clebsch::weight::index() const {
  89     int result = 0;
  90
  91     for (int i = 0; elem[i] > elem[N - 1]; ++i) {
  92         result += binomial(N - i - 1 + elem[i] - elem[N - 1] - 1, N - i - 1);
  93     }
  94
  95     return result;
  96 }
  97
  98 long long clebsch::weight::dimension() const {
  99     long long numerator = 1, denominator = 1;
 100
 101     for (int i = 1; i < N; ++i) {
 102         for (int j = 0; i + j < N; ++j) {
 103             numerator *= elem[j] - elem[i + j] + i;
 104             denominator *= i;
 105         }
 106     }
 107
 108     return numerator / denominator;
 109 }
 110
 111 // implementation of "pattern" starts here
 112
 113 clebsch::pattern::pattern(const pattern &p) : elem(p.elem), N(p.N) {}
 114
 115 clebsch::pattern::pattern(const weight &irrep, int index) :
 116         elem((irrep.N * (irrep.N + 1)) / 2), N(irrep.N) {
 117     for (int i = 1; i <= N; ++i) {
 118         (*this)(i, N) = irrep(i);
 119     }
 120
 121     for (int l = N - 1; l >= 1; --l) {
 122         for (int k = 1; k <= l; ++k) {
 123             (*this)(k, l) = (*this)(k + 1, l + 1);
 124         }
 125     }
 126
 127     while (index-- > 0) {
 128         bool b = ++(*this);
 129
 130         assert(b);
 131     }
 132 }
 133
 134 int &clebsch::pattern::operator()(int k, int l) {
 135     return elem[(N * (N + 1) - l * (l + 1)) / 2 + k - 1];
 136 }
 137
 138 const int &clebsch::pattern::operator()(int k, int l) const {
 139     return elem[(N * (N + 1) - l * (l + 1)) / 2 + k - 1];
 140 }
 141
 142 bool clebsch::pattern::operator++() {
 143     int k = 1, l = 1;
 144
 145     while (l < N && (*this)(k, l) == (*this)(k, l + 1)) {
 146         if (--k == 0) {
 147             k = ++l;
 148         }
 149     }
 150
 151     if (l == N) {
 152         return false;
 153     }
 154
 155     ++(*this)(k, l);
 156
 157     while (k != 1 || l != 1) {
 158         if (++k > l) {
 159             k = 1;
 160             --l;
 161         }
 162
 163         (*this)(k, l) = (*this)(k + 1, l + 1);
 164     }
 165
 166     return true;
 167 }
 168
 169 bool clebsch::pattern::operator--() {
 170     int k = 1, l = 1;
 171
 172     while (l < N && (*this)(k, l) == (*this)(k + 1, l + 1)) {
 173         if (--k == 0) {
 174             k = ++l;
 175         }
 176     }
 177
 178     if (l == N) {
 179         return false;
 180     }
 181
 182     --(*this)(k, l);
 183
 184     while (k != 1 || l != 1) {
 185         if (++k > l) {
 186             k = 1;
 187             --l;
 188         }
 189
 190         (*this)(k, l) = (*this)(k, l + 1);
 191     }
 192
 193     return true;
 194 }
 195
 196 int clebsch::pattern::index() const {
 197     int result = 0;
 198
 199     for (pattern p(*this); --p; ++result) {}
 200
 201     return result;
 202 }
 203
 204 clebsch::weight clebsch::pattern::get_weight() const {
 205     clebsch::weight result(N);
 206
 207     for (int prev = 0, l = 1; l <= N; ++l) {
 208         int now = 0;
 209
 210         for (int k = 1; k <= l; ++k) {
 211             now += (*this)(k, l);
 212         }
 213
 214         result(l) = now - prev;
 215         prev = now;
 216     }
 217
 218     return result;
 219 }
 220
 221 double clebsch::pattern::lowering_coeff(int k, int l) const {
 222     double result = 1.0;
 223
 224     for (int i = 1; i <= l + 1; ++i) {
 225         result *= (*this)(i, l + 1) - (*this)(k, l) + k - i + 1;
 226     }
 227
 228     for (int i = 1; i <= l - 1; ++i) {
 229         result *= (*this)(i, l - 1) - (*this)(k, l) + k - i;
 230     }
 231
 232     for (int i = 1; i <= l; ++i) {
 233         if (i == k) continue;
 234         result /= (*this)(i, l) - (*this)(k, l) + k - i + 1;
 235         result /= (*this)(i, l) - (*this)(k, l) + k - i;
 236     }
 237
 238     return std::sqrt(-result);
 239 }
 240
 241 double clebsch::pattern::raising_coeff(int k, int l) const {
 242     double result = 1.0;
 243
 244     for (int i = 1; i <= l + 1; ++i) {
 245         result *= (*this)(i, l + 1) - (*this)(k, l) + k - i;
 246     }
 247
 248     for (int i = 1; i <= l - 1; ++i) {
 249         result *= (*this)(i, l - 1) - (*this)(k, l) + k - i - 1;
 250     }
 251
 252     for (int i = 1; i <= l; ++i) {
 253         if (i == k) continue;
 254         result /= (*this)(i, l) - (*this)(k, l) + k - i;
 255         result /= (*this)(i, l) - (*this)(k, l) + k - i  - 1;
 256     }
 257
 258     return std::sqrt(-result);
 259 }
 260
 261 // implementation of "decomposition" starts here
 262
 263 clebsch::decomposition::decomposition(const weight &factor1, const weight &factor2) :
 264         N(factor1.N), factor1(factor1), factor2(factor2) {
 265     assert(factor1.N == factor2.N);
 266     std::vector<clebsch::weight> result;
 267     pattern low(factor1), high(factor1);
 268     weight trial(factor2);
 269     int k = 1, l = N;
 270
 271     do {
 272         while (k <= N) {
 273             --l;
 274             if (k <= l) {
 275                 low(k, l) = std::max(high(k + N - l, N), high(k, l + 1) + trial(l + 1) - trial(l));
 276                 high(k, l) = high(k, l + 1);
 277                 if (k > 1 && high(k, l) > high(k - 1, l - 1)) {
 278                     high(k, l) = high(k - 1, l - 1);
 279                 }
 280                 if (l > 1 && k == l && high(k, l) > trial(l - 1) - trial(l)) {
 281                     high(k, l) = trial(l - 1) - trial(l);
 282                 }
 283                 if (low(k, l) > high(k, l)) {
 284                     break;
 285                 }
 286                 trial(l + 1) += high(k, l + 1) - high(k, l);
 287             } else {
 288                 trial(l + 1) += high(k, l + 1);
 289                 ++k;
 290                 l = N;
 291             }
 292         }
 293
 294         if (k > N) {
 295             result.push_back(trial);
 296             for (int i = 1; i <= N; ++i) {
 297                 result.back()(i) -= result.back()(N);
 298             }
 299         } else {
 300             ++l;
 301         }
 302
 303         while (k != 1 || l != N) {
 304             if (l == N) {
 305                 l = --k - 1;
 306                 trial(l + 1) -= high(k, l + 1);
 307             } else if (low(k, l) < high(k, l)) {
 308                 --high(k, l);
 309                 ++trial(l + 1);
 310                 break;
 311             } else {
 312                 trial(l + 1) -= high(k, l + 1) - high(k, l);
 313             }
 314             ++l;
 315         }
 316     } while (k != 1 || l != N);
 317
 318     sort(result.begin(), result.end());
 319     for (std::vector<clebsch::weight>::iterator it = result.begin(); it != result.end(); ++it) {
 320         if (it != result.begin() && *it == weights.back()) {
 321             ++multiplicities.back();
 322         } else {
 323             weights.push_back(*it);
 324             multiplicities.push_back(1);
 325         }
 326     }
 327 }
 328
 329 int clebsch::decomposition::size() const {
 330     return weights.size();
 331 }
 332
 333 const clebsch::weight &clebsch::decomposition::operator()(int j) const {
 334     return weights[j];
 335 }
 336
 337 int clebsch::decomposition::multiplicity(const weight &irrep) const {
 338     assert(irrep.N == N);
 339     std::vector<clebsch::weight>::const_iterator it
 340         = std::lower_bound(weights.begin(), weights.end(), irrep);
 341
 342     return it != weights.end() && *it == irrep ? multiplicities[it - weights.begin()] : 0;
 343 }
 344
 345 // implementation of "index_adapter" starts here
 346
 347 clebsch::index_adapter::index_adapter(const clebsch::decomposition &decomp) :
 348         N(decomp.N),
 349         factor1(decomp.factor1.index()),
 350         factor2(decomp.factor2.index()) {
 351     for (int i = 0, s = decomp.size(); i < s; ++i) {
 352         indices.push_back(decomp(i).index());
 353         multiplicities.push_back(decomp.multiplicity(decomp(i)));
 354     }
 355 }
 356
 357 int clebsch::index_adapter::size() const {
 358     return indices.size();
 359 }
 360
 361 int clebsch::index_adapter::operator()(int j) const {
 362     return indices[j];
 363 }
 364
 365 int clebsch::index_adapter::multiplicity(int irrep) const {
 366     std::vector<int>::const_iterator it = std::lower_bound(indices.begin(), indices.end(), irrep);
 367
 368     return it != indices.end() && *it == irrep ? multiplicities[it - indices.begin()] : 0;
 369 }
 370
 371 // implementation of "clebsch" starts here
 372
 373 void clebsch::coefficients::set(int factor1_state,
 374                                 int factor2_state,
 375                                 int multiplicity_index,
 376                                 int irrep_state,
 377                                 double value) {
 378     assert(0 <= factor1_state && factor1_state < factor1_dimension);
 379     assert(0 <= factor2_state && factor2_state < factor2_dimension);
 380     assert(0 <= multiplicity_index && multiplicity_index < multiplicity);
 381     assert(0 <= irrep_state && irrep_state < irrep_dimension);
 382
 383     int coefficient_label[] = { factor1_state,
 384                                 factor2_state,
 385                                 multiplicity_index,
 386                                 irrep_state };
 387     clzx[std::vector<int>(coefficient_label, coefficient_label
 388             + sizeof coefficient_label / sizeof coefficient_label[0])] = value;
 389 }
 390
 391 void clebsch::coefficients::highest_weight_normal_form() {
 392     int hws = irrep_dimension - 1;
 393
 394     // bring CGCs into reduced row echelon form
 395     for (int h = 0, i = 0; h < multiplicity - 1 && i < factor1_dimension; ++i) {
 396         for (int j = 0; h < multiplicity - 1 && j < factor2_dimension; ++j) {
 397             int k0 = h;
 398
 399             for (int k = h + 1; k < multiplicity; ++k) {
 400                 if (fabs((*this)(i, j, k, hws)) > fabs((*this)(i, j, k0, hws))) {
 401                     k0 = k;
 402                 }
 403             }
 404
 405             if ((*this)(i, j, k0, hws) < -EPS) {
 406                 for (int i2 = i; i2 < factor1_dimension; ++i2) {
 407                     for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) {
 408                         set(i2, j2, k0, hws, -(*this)(i2, j2, k0, hws));
 409                     }
 410                 }
 411             } else if ((*this)(i, j, k0, hws) < EPS) {
 412                 continue;
 413             }
 414
 415             if (k0 != h) {
 416                 for (int i2 = i; i2 < factor1_dimension; ++i2) {
 417                     for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) {
 418                         double x = (*this)(i2, j2, k0, hws);
 419                         set(i2, j2, k0, hws, (*this)(i2, j2, h, hws));
 420                         set(i2, j2, h, hws, x);
 421                     }
 422                 }
 423             }
 424
 425             for (int k = h + 1; k < multiplicity; ++k) {
 426                 for (int i2 = i; i2 < factor1_dimension; ++i2) {
 427                     for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) {
 428                         set(i2, j2, k, hws, (*this)(i2, j2, k, hws) - (*this)(i2, j2, h, hws)
 429                                 * (*this)(i, j, k, hws) / (*this)(i, j, h, hws));
 430                     }
 431                 }
 432             }
 433
 434             // next 3 lines not strictly necessary, might improve numerical stability
 435             for (int k = h + 1; k < multiplicity; ++k) {
 436                 set(i, j, k, hws, 0.0);
 437             }
 438
 439             ++h;
 440         }
 441     }
 442
 443     // Gram-Schmidt orthonormalization
 444     for (int h = 0; h < multiplicity; ++h) {
 445         for (int k = 0; k < h; ++k) {
 446             double overlap = 0.0;
 447             for (int i = 0; i < factor1_dimension; ++i) {
 448                 for (int j = 0; j < factor2_dimension; ++j) {
 449                     overlap += (*this)(i, j, h, hws) * (*this)(i, j, k, hws);
 450                 }
 451             }
 452
 453             for (int i = 0; i < factor1_dimension; ++i) {
 454                 for (int j = 0; j < factor2_dimension; ++j) {
 455                     set(i, j, h, hws, (*this)(i, j, h, hws) - overlap * (*this)(i, j, k, hws));
 456                 }
 457             }
 458         }
 459
 460         double norm = 0.0;
 461         for (int i = 0; i < factor1_dimension; ++i) {
 462             for (int j = 0; j < factor2_dimension; ++j) {
 463                 norm += (*this)(i, j, h, hws) * (*this)(i, j, h, hws);
 464             }
 465         }
 466         norm = std::sqrt(norm);
 467
 468         for (int i = 0; i < factor1_dimension; ++i) {
 469             for (int j = 0; j < factor2_dimension; ++j) {
 470                 set(i, j, h, hws, (*this)(i, j, h, hws) / norm);
 471             }
 472         }
 473     }
 474 }
 475
 476 void clebsch::coefficients::compute_highest_weight_coeffs() {
 477     if (multiplicity == 0) {
 478         return;
 479     }
 480
 481     std::vector<std::vector<int> > map_coeff(factor1_dimension,
 482                                              std::vector<int>(factor2_dimension, -1));
 483     std::vector<std::vector<int> > map_states(factor1_dimension,
 484                                               std::vector<int>(factor2_dimension, -1));
 485     int n_coeff = 0, n_states = 0;
 486     pattern p(factor1, 0);
 487
 488     for (int i = 0; i < factor1_dimension; ++i, ++p) {
 489         weight pw(p.get_weight());
 490         pattern q(factor2, 0);
 491         for (int j = 0; j < factor2_dimension; ++j, ++q) {
 492             if (pw + q.get_weight() == irrep) {
 493                 map_coeff[i][j] = n_coeff++;
 494             }
 495         }
 496     }
 497
 498     if (n_coeff == 1) {
 499         for (int i = 0; i < factor1_dimension; ++i) {
 500             for (int j = 0; j < factor2_dimension; ++j) {
 501                 if (map_coeff[i][j] >= 0) {
 502                     set(i, j, 0, irrep_dimension - 1, 1.0);
 503                     return;
 504                 }
 505             }
 506         }
 507     }
 508
 509     double *hw_system = new double[n_coeff * (factor1_dimension * factor2_dimension)];
 510     assert(hw_system != NULL);
 511     memset(hw_system, 0, n_coeff * (factor1_dimension * factor2_dimension) * sizeof (double));
 512
 513     pattern r(factor1, 0);
 514     for (int i = 0; i < factor1_dimension; ++i, ++r) {
 515         pattern q(factor2, 0);
 516
 517         for (int j = 0; j < factor2_dimension; ++j, ++q) {
 518             if (map_coeff[i][j] >= 0) {
 519                 for (int l = 1; l <= N - 1; ++l) {
 520                     for (int k = 1; k <= l; ++k) {
 521                         if ((k == 1 || r(k, l) + 1 <= r(k - 1, l - 1)) && r(k, l) + 1 <= r(k, l + 1)) {
 522                             ++r(k, l);
 523                             int h = r.index();
 524                             --r(k, l);
 525
 526                             if (map_states[h][j] < 0) {
 527                                 map_states[h][j] = n_states++;
 528                             }
 529
 530                             hw_system[n_coeff * map_states[h][j] + map_coeff[i][j]]
 531                                 += r.raising_coeff(k, l);
 532                         }
 533
 534                         if ((k == 1 || q(k, l) + 1 <= q(k - 1, l - 1)) && q(k, l) + 1 <= q(k, l + 1)) {
 535                             ++q(k, l);
 536                             int h = q.index();
 537                             --q(k, l);
 538
 539                             if (map_states[i][h] < 0) {
 540                                 map_states[i][h] = n_states++;
 541                             }
 542
 543
 544                             hw_system[n_coeff * map_states[i][h] + map_coeff[i][j]]
 545                                 += q.raising_coeff(k, l);
 546                         }
 547                     }
 548                 }
 549             }
 550         }
 551     }
 552
 553     int lwork = -1, info;
 554     double worksize;
 555
 556     double *singval = new double[std::min(n_coeff, n_states)];
 557     assert(singval != NULL);
 558     double *singvec = new double[n_coeff * n_coeff];
 559     assert(singvec != NULL);
 560
 561     dgesvd_("A",
 562             "N",
 563             &n_coeff,
 564             &n_states,
 565             hw_system,
 566             &n_coeff,
 567             singval,
 568             singvec,
 569             &n_coeff,
 570             NULL,
 571             &n_states,
 572             &worksize,
 573             &lwork,
 574             &info);
 575     assert(info == 0);
 576
 577     lwork = worksize;
 578     double *work = new double[lwork];
 579     assert(work != NULL);
 580
 581     dgesvd_("A",
 582             "N",
 583             &n_coeff,
 584             &n_states,
 585             hw_system,
 586             &n_coeff,
 587             singval,
 588             singvec,
 589             &n_coeff,
 590             NULL,
 591             &n_states,
 592             work,
 593             &lwork,
 594             &info);
 595     assert(info == 0);
 596
 597     for (int i = 0; i < multiplicity; ++i) {
 598         for (int j = 0; j < factor1_dimension; ++j) {
 599             for (int k = 0; k < factor2_dimension; ++k) {
 600                 if (map_coeff[j][k] >= 0) {
 601                     double x = singvec[n_coeff * (n_coeff - 1 - i) + map_coeff[j][k]];
 602
 603                     if (fabs(x) > EPS) {
 604                         set(j, k, i, irrep_dimension - 1, x);
 605                     }
 606                 }
 607             }
 608         }
 609     }
 610
 611     // uncomment next line to bring highest-weight coefficients into "normal form"
 612     // highest_weight_normal_form();
 613
 614     delete[] work;
 615     delete[] singvec;
 616     delete[] singval;
 617     delete[] hw_system;
 618 }
 619
 620 void clebsch::coefficients::compute_lower_weight_coeffs(int multip_index,
 621                                                         int state,
 622                                                         std::vector<char> &done) {
 623     weight statew(pattern(irrep, state).get_weight());
 624     pattern p(irrep, 0);
 625     std::vector<int> map_parent(irrep_dimension, -1),
 626                      map_multi(irrep_dimension, -1),
 627                      which_l(irrep_dimension, -1);
 628     int n_parent = 0, n_multi = 0;
 629
 630     for (int i = 0; i < irrep_dimension; ++i, ++p) {
 631         weight v(p.get_weight());
 632
 633         if (v == statew) {
 634             map_multi[i] = n_multi++;
 635         } else for (int l = 1; l < N; ++l) {
 636             --v(l);
 637             ++v(l + 1);
 638             if (v == statew) {
 639                 map_parent[i] = n_parent++;
 640                 which_l[i] = l;
 641                 if (!done[i]) {
 642                     compute_lower_weight_coeffs(multip_index, i, done);
 643                 }
 644                 break;
 645             }
 646             --v(l + 1);
 647             ++v(l);
 648         }
 649     }
 650
 651     double *irrep_coeffs = new double[n_parent * n_multi];
 652     assert(irrep_coeffs != NULL);
 653     memset(irrep_coeffs, 0, n_parent * n_multi * sizeof (double));
 654
 655     double *prod_coeffs = new double[n_parent * factor1_dimension * factor2_dimension];
 656     assert(prod_coeffs != NULL);
 657     memset(prod_coeffs, 0, n_parent * factor1_dimension * factor2_dimension * sizeof (double));
 658
 659     std::vector<std::vector<int> > map_prodstat(factor1_dimension,
 660                                                 std::vector<int>(factor2_dimension, -1));
 661     int n_prodstat = 0;
 662
 663     pattern r(irrep, 0);
 664     for (int i = 0; i < irrep_dimension; ++i, ++r) {
 665         if (map_parent[i] >= 0) {
 666             for (int k = 1, l = which_l[i]; k <= l; ++k) {
 667                 if (r(k, l) > r(k + 1, l + 1) && (k == l || r(k, l) > r(k, l - 1))) {
 668                     --r(k, l);
 669                     int h = r.index();
 670                     ++r(k, l);
 671
 672                     irrep_coeffs[n_parent * map_multi[h] + map_parent[i]] += r.lowering_coeff(k, l);
 673                 }
 674             }
 675
 676             pattern q1(factor1, 0);
 677             for (int j1 = 0; j1 < factor1_dimension; ++j1, ++q1) {
 678                 pattern q2(factor2, 0);
 679
 680                 for (int j2 = 0; j2 < factor2_dimension; ++j2, ++q2) {
 681                     if (std::fabs((*this)(j1, j2, multip_index, i)) > EPS) {
 682                         for (int k = 1, l = which_l[i]; k <= l; ++k) {
 683                             if (q1(k, l) > q1(k + 1, l + 1) && (k == l || q1(k, l) > q1(k, l - 1))) {
 684                                 --q1(k, l);
 685                                 int h = q1.index();
 686                                 ++q1(k, l);
 687
 688                                 if (map_prodstat[h][j2] < 0) {
 689                                     map_prodstat[h][j2] = n_prodstat++;
 690                                 }
 691
 692                                 prod_coeffs[n_parent * map_prodstat[h][j2] + map_parent[i]] +=
 693                                         (*this)(j1, j2, multip_index, i) * q1.lowering_coeff(k, l);
 694                             }
 695
 696                             if (q2(k, l) > q2(k + 1, l + 1) && (k == l || q2(k, l) > q2(k, l - 1))) {
 697                                 --q2(k, l);
 698                                 int h = q2.index();
 699                                 ++q2(k, l);
 700
 701                                 if (map_prodstat[j1][h] < 0) {
 702                                     map_prodstat[j1][h] = n_prodstat++;
 703                                 }
 704
 705                                 prod_coeffs[n_parent * map_prodstat[j1][h] + map_parent[i]] +=
 706                                         (*this)(j1, j2, multip_index, i) * q2.lowering_coeff(k, l);
 707                             }
 708                         }
 709                     }
 710                 }
 711             }
 712         }
 713     }
 714
 715     double worksize;
 716     int lwork = -1, info;
 717
 718     dgels_("N",
 719            &n_parent,
 720            &n_multi,
 721            &n_prodstat,
 722            irrep_coeffs,
 723            &n_parent,
 724            prod_coeffs,
 725            &n_parent,
 726            &worksize,
 727            &lwork,
 728            &info);
 729     assert(info == 0);
 730
 731     lwork = worksize;
 732     double *work = new double[lwork];
 733     assert(work != NULL);
 734
 735     dgels_("N",
 736            &n_parent,
 737            &n_multi,
 738            &n_prodstat,
 739            irrep_coeffs,
 740            &n_parent,
 741            prod_coeffs,
 742            &n_parent,
 743            work,
 744            &lwork,
 745            &info);
 746     assert(info == 0);
 747
 748     for (int i = 0; i < irrep_dimension; ++i) {
 749         if (map_multi[i] >= 0) {
 750             for (int j = 0; j < factor1_dimension; ++j) {
 751                 for (int k = 0; k < factor2_dimension; ++k) {
 752                     if (map_prodstat[j][k] >= 0) {
 753                         double x = prod_coeffs[n_parent * map_prodstat[j][k] + map_multi[i]];
 754
 755                         if (fabs(x) > EPS) {
 756                             set(j, k, multip_index, i, x);
 757                         }
 758                     }
 759                 }
 760             }
 761
 762             done[i] = true;
 763         }
 764     }
 765
 766     delete[] work;
 767     delete[] prod_coeffs;
 768     delete[] irrep_coeffs;
 769 }
 770
 771 clebsch::coefficients::coefficients(const weight &irrep, const weight &factor1, const weight &factor2) :
 772         N(irrep.N),
 773         factor1(factor1),
 774         factor2(factor2),
 775         irrep(irrep),
 776         factor1_dimension(factor1.dimension()),
 777         factor2_dimension(factor2.dimension()),
 778         irrep_dimension(irrep.dimension()),
 779         multiplicity(decomposition(factor1, factor2).multiplicity(irrep)) {
 780     assert(factor1.N == irrep.N);
 781     assert(factor2.N == irrep.N);
 782
 783     compute_highest_weight_coeffs();
 784
 785     for (int i = 0; i < multiplicity; ++i) {
 786         std::vector<char> done(irrep_dimension, 0);
 787         done[irrep_dimension - 1] = true;
 788         for (int j = irrep_dimension - 1; j >= 0; --j) {
 789             if (!done[j]) {
 790                 compute_lower_weight_coeffs(i, j, done);
 791             }
 792         }
 793     }
 794 }
 795
 796 double clebsch::coefficients::operator()(int factor1_state,
 797                                          int factor2_state,
 798                                          int multiplicity_index,
 799                                          int irrep_state) const {
 800
 801   assert(0 <= factor1_state && factor1_state < factor1_dimension);
 802   assert(0 <= factor2_state && factor2_state < factor2_dimension);
 803   assert(0 <= multiplicity_index && multiplicity_index < multiplicity);
 804   assert(0 <= irrep_state && irrep_state < irrep_dimension);
 805
 806   int coefficient_label[] = { factor1_state,
 807                               factor2_state,
 808                               multiplicity_index,
 809                               irrep_state };
 810   std::map<std::vector<int>, double>::const_iterator it(
 811             clzx.find(std::vector<int>(coefficient_label, coefficient_label
 812                                        + sizeof coefficient_label / sizeof coefficient_label[0])));
 813
 814   return it != clzx.end() ? it->second : 0.0;
 815 }
 816