Умножить на супернодальный L в CHOLMOD?

Question

Как я могу умножить на cholmod_factor L в супернодальной L L^T факторизации? Я бы предпочел не преобразовывать в симплициальные, так как супернодальное представление приводит к более быстрым обратным вычислениям, и я бы предпочел не делать копию фактора, поскольку две копии могут не помещаться в ОЗУ.

Answer 1

Я получил представление о супернодальном представлении из хорошего комментария в хелперной функции от супернодального к симплициальному в t_cholmod_change_factor.c, Я перефразирую комментарий и добавляю некоторые детали ниже:

Супер-узловая факторизация Холецкого представлена ​​как совокупность супер-узловых блоков. Записи супернодального блока расположены в главном порядке столбцов, как этот суперузел 6x4:

t - - -    (row s[pi[snode+0]])
t t - -    (row s[pi[snode+1]])
t t t -    (row s[pi[snode+2]])
t t t t    (row s[pi[snode+3]])
r r r r    (row s[pi[snode+4]])
r r r r    (row s[pi[snode+5]])

• Есть неиспользуемые записи (обозначенные дефисами), чтобы сделать матрицу прямоугольной.
• Индексы столбцов являются последовательными.
• Первый ncols индексы строк - это те же самые последовательные индексы столбцов. Более поздние индексы строк могут ссылаться на любой ряд ниже t треугольник.
• super член имеет одну запись для каждого суперузла; это относится к первому столбцу, представленному суперузлом.
• pi член имеет одну запись для каждого суперузла; это относится к первому индексу в s член, где вы можете посмотреть номера строк.
• px член имеет одну запись для каждого суперузла; это относится к первому индексу в x член, где хранятся записи. Опять же, это не упакованное хранилище.

Следующий код для умножения на cholmod_factor *L кажется, работает (я забочусь только о int индексы и реальные записи двойной точности):

cholmod_dense *mul_L(cholmod_factor *L, cholmod_dense *d) {
  int rows = d->nrow, cols = d->ncol;
  cholmod_dense *ans = cholmod_allocate_dense(rows, cols, rows,
      CHOLMOD_REAL, &comm);
  memset(ans->x, 0, 8 * rows * cols);

  FOR(i, L->nsuper) {
    int *sup = (int *)L->super;
    int *pi = (int *)L->pi;
    int *px = (int *)L->px;
    double *x = (double *)L->x;
    int *ss = (int *)L->s;

    int r0 =  pi[i], r1 =  pi[i+1], nrow = r1 - r0;
    int c0 = sup[i], c1 = sup[i+1], ncol = c1 - c0;
    int px0 = px[i];

    /* TODO: Use BLAS instead. */
    for (int j = 0; j < ncol; j++) {
      for (int k = j; k < nrow; k++) {
        for (int l = 0; l < cols; l++) {
          ((double *)ans->x)[l * rows + ss[r0 + k]] +=
              x[px0 + k + j * nrow] * ((double *)d->x)[l*rows+c0 + j];
        }
      }
    }
  }
  return ans;
}