2013-02-11 120 views
1

我最近安装了Cholmod以便在一些C++代码中执行稀疏胆固醇分解。我想然后使用分解来计算矩阵逆(I有以下问题:使用Cholmod和Cholmod-Extra计算稀疏矩阵的逆

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

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

其中I想要计算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++中稀疏矩阵的逆。

干杯 林德利

回答

1

我以前用过JAMA。它有乔列斯基分解。

+0

JAMA支持稀疏矩阵吗?该代码似乎在密集矩阵中执行普通的cholesky分解,在这种情况下,我只是使用lapack或其他东西。 – LindleyLentati 2013-02-11 19:03:56

+0

美国医学会杂志使用TNT(模板数值工具包),它支持稀疏矩阵http://math.nist.gov/tnt/overview.html#structures,所以我假设它 – Smash 2013-02-11 19:47:54