gpt4 book ai didi

c++ - 使用 Cholmod 和 Cholmod-Extra 计算稀疏矩阵的逆

转载 作者:太空宇宙 更新时间:2023-11-04 11:59:15 25 4
gpt4 key购买 nike

我最近安装了 Cholmod 以便在某些 C++ 代码中执行稀疏 cholesky 分解。然后我想使用分解来计算矩阵逆(我有以下问题:

d^T . (A^-1 + B^-1)^-1 . d

其中d是一个 vector ^T表示转置,AB是稀疏的)

我想计算 B 的实际倒数,然后可以线性求解涉及总和的解。调用它的代码如下:

cholmod_common Common, *cm ;
cm = &Common ;
cholmod_start (cm) ;
cm->print=5;
Common.supernodal = CHOLMOD_SUPERNODAL ;
double *Tx, x;
int *Ti, *Tj, *Rdeg, *Cdeg,j;
cholmod_triplet *T ;
size_t nrow;
size_t ncol;
size_t nnz ;

nrow=Csize;
ncol=Csize;
nnz=numpulsars*numpulsars*numcoeff;

T = cholmod_allocate_triplet(nrow,ncol,nnz,0,CHOLMOD_REAL, cm) ;
Ti=(int*)T->i;
Tj=(int*)T->j;
Tx=(double*)T->x;

for(int i=0;i<numpulsars;i++){
for(int j=0;j<numpulsars;j++){
if(i==j){
pcorr=1.0;
}
else{
angle=pulsarCartesian[i][0]*pulsarCartesian[j][0] +pulsarCartesian[i][1]*pulsarCartesian[j][1]+pulsarCartesian[i][2]*pulsarCartesian[j][2];
pcorr=(1.5*(0.5*(1-angle))*log(0.5*(1-angle)) - 0.25*0.5*(1-angle) + 0.5);
}

for(int c1=0; c1<numcoeff; c1++){
val= pcorr*powercoeff[c1];
Ti[T->nnz]=i*numcoeff+c1;
Tj[T->nnz]=j*numcoeff+c1;
Tx[T->nnz]=val;
(T->nnz)++;


}
}
}


cholmod_sparse *PPFMSparse;
cholmod_factor *L ;
cholmod_dense *spinvK;
PPFMSparse=cholmod_triplet_to_sparse(T,T->nnz,cm);
// cholmod_print_sparse(PPFMSparse, "PPFMSparse", cm);
L = cholmod_analyze (PPFMSparse, cm) ;
cholmod_factorize (PPFMSparse, L, cm) ;

cholmod_sparse *PPFMinv;

PPFMinv=cholmod_spinv(L,cm);
// cholmod_print_sparse(PPFMinv, "PPFMinv", cm);
spinvK = cholmod_sparse_to_dense(PPFMinv, &Common) ;
cholmod_print_dense(spinvK, "spinvK", cm);
cholmod_free_sparse(&PPFMinv,cm);
cholmod_free_factor (&L, cm) ;
cholmod_free_sparse(&PPFMSparse,cm);
cholmod_free_triplet (&T, cm) ;
cholmod_free_dense (&spinvK, cm) ;
cholmod_finish(cm);

我找到了 https://cholmod-extra.readthedocs.org/en/latest/functions.html这个函数本来是用来计算倒数的,但它给了我一些与倒数平方相关的东西,而不是倒数。我只是想知道是否有人有使用它的经验,或者用 C++ 计算稀疏矩阵的逆矩阵的等效方法。

干杯林德利

最佳答案

我以前用过JAMA .它具有 Cholesky 分解。

关于c++ - 使用 Cholmod 和 Cholmod-Extra 计算稀疏矩阵的逆,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/14818465/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com