在超节L^T分解中,我如何乘以cholmod_factor L?我不想转换成单纯的,因为超节点表示会导致更快的反解,而且我也不想复制因子,因为两个副本可能不适合RAM。
发布于 2014-08-04 23:48:07
最后,我从t_cholmod_change_factor.c中超节点到单纯形辅助函数中的一个很好的注释中理解了超节点表示。我解释了这一评论,并在下面添加了一些细节:
超节点Cholesky分解被表示为超节点块的集合。超节点块的条目按列的主要顺序排列,如这个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;
}https://stackoverflow.com/questions/25122076
复制相似问题