LAPACK QR因子分解

我試圖從LAPACK zgeqrf和zungqr例程中獲得Q-Matrix。LAPACK QR因子分解

我有一個Nw-by-Nw複數矩陣，其列上有非正交向量。這是我的C++代碼。（矩陣名爲vr_tr）

//QR Fact. 
complex<double> TAU[Nw*Nw]; 
zgeqrf_(&Nw, &Nw, vr_tr, &Nw, TAU, &wkopt, &lwork, &info); 
lwork = (int)wkopt.real(); 
work = new complex<double> [lwork]; 
zgeqrf_(&Nw, &Nw, vr_tr, &Nw, TAU, work, &lwork, &info); 
zungqr_(&Nw, &Nw, &Nw, vr_tr, &Nw, TAU, &wkopt, &lwork, &info); 
lwork = (int)wkopt.real(); 
work = new complex<double> [lwork]; 
zungqr_(&Nw, &Nw, &Nw, vr_tr, &Nw, TAU, work, &lwork, &info); 

//Checking if vr_tr * vr_tr' = I 
complex<double> one = 1,zer=0,res[Nw*Nw]; char Tchar = 'T', Nchar = 'N'; 
zgemm_(&Nchar, &Tchar, &Nw, &Nw, &Nw, &one, vr_tr, &Nw, vr_tr, &Nw, &zer, res, &Nw); 
for(i=0;i<Nw;i++){ 
    for(j=0;j<Nw;j++){ 
     cout<<res[i*Nw+j]<<"\t"; 
    } 
    cout<<"\n\n"; 
}

什麼運行此代碼後，我得到的是不是單位矩陣，因爲它應該是，因爲它應該從zgeqrf計算QR分解，並從zungqr和Q獲得Q-矩陣-matrix必須是正交的，所以Q*Q'=I。

這段代碼有什麼問題？

來源

2017-08-07 Alireza

出於好奇，當高級C++包裝器存在時，爲什麼你會用lapack打擾，像[Armadillo]（http://arma.sourceforge.net/docs.html#qr）？ –

@TheQuantumPhysicist習慣，因爲我工作的人使用它。儘管沒嘗試過犰狳。它可以執行所有LAPACK例程嗎？ – Alireza

檢查我提供的鏈接中的文檔。當我還在學術界進行研究時，它有我需要的所有例程，甚至更多。例如，LAPACK沒有一個通用的矩陣求冪函數（除非你想通過計算特徵值並假設你有一個對稱矩陣來完成它），但是Armadillo確實有這個功能。您可以將Armadillo鏈接到您希望的任何LAPACK實現，並且可以爲您節省很多時間。 Armadillo唯一的問題是它不支持集羣（AFAIK）。否則，你應該去做。 –

我使用zunmqr而不是znugqr它得到了修復。

雖然我真的不知道爲什麼zungqr沒有工作。

來源

2017-08-07 18:43:11 Alireza

LAPACK QR因子分解

回答

相關問題