2013-02-11 135 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