求解方程的帶狀矩陣系統

我需要求解二維泊松方程，也就是說，對於AX = B的方程組，其中A是一個n×n矩陣，B是一個n×n的矩陣， 1矢量。作爲二維泊松問題的離散化矩陣，我知道只有5個對角線不會爲空。 Lapack不提供解決這個特定問題的函數，但它具有求解帶狀矩陣方程組的功能，即DGBTRF（用於LU分解）和DGBTRS。現在，5個對角線是：主對角線，主對角線上方和下方的第一個對角線，以及對角線上下兩對角線與主對角線對角線。在閱讀關於波段存儲的lapack文檔後，我瞭解到我必須創建一個（3 * m + 1）n矩陣來存儲帶內存儲格式，我們稱之爲矩陣AB。現在的問題：求解方程的帶狀矩陣系統

1）dgbtrs和dgbtrs_有什麼區別？英特爾MKL提供了這兩種功能，但我不明白爲什麼

2）dgbtrf要求帶寬存儲矩陣是一個數組。我應該按行還是按列對AB進行線性化？

3）這是調用這兩個函數的正確方法嗎？

int n, m; 
double *AB; 
/*... fill n, m, AB, with appropriate numbers */ 

int *pivots; 
int nrows = 3 * m + 1, info, rhs = 1; 
dgbtrf_(&n, &n, &m, &m, AB, &nrows, pivots, &info); 
char trans = 'N'; 
dgbtrs_(&trans, &n, &m, &m, &rhs, AB, &nrows, pivots, B, &n, &info);

來源

2012-06-06 Patrik

它還提供DGBTRS和DGBTRS_。那些是你不應該關心的fortran administrativa。請致電dgbtrs（原因在於，在某些體系結構中，Fortran例程名稱附加了下劃線，而其他名稱不能，名稱可能是大寫或小寫 - 英特爾MKL #define是dgbtrs的右列）。
LAPACK例程需要列主矩陣（即Fortran樣式）：將列存儲在其他列之後。您必須使用的帶狀存儲不難：http://www.netlib.org/lapack/lug/node124.html。
對我來說這似乎很好，但請事先嚐試一下小問題（順便提一下，總是一個好主意）。還要確保你處理非零值info（這是錯誤報告的方式）。

更好的風格是使用MKL_INT而不是純int，這是一個typedef向右類型（可以是在某些體系結構不同）。

在致電dgbtrf之前，請確保爲pivots分配內存。

來源

2012-06-06 19:45:52

關於＃2，是不是會搞亂LU分解？我的意思是，沒有選擇告訴dgbtrf AB是行主要的。另外，我如何設定兩種情況下的主要維度？ – Patrik

是的，你會解決A'x = B。讓我糾正我的答案。對於主要維度，它們通常等於行數，除非您正在做一些有趣的事情（如將子矩陣作爲參數傳遞）。再一次，解決LAPACK問題的最好方法是嘗試手動解決小問題。 –

這可能是主題。但對於泊松方程，基於FFT的解決方案要快得多。只需做你的勢場的二維FFT，除以 - （k^2 + lambda^2），然後做IFFT。 λ是避免k = 0發散的小數。 5對角線方程是泊松方程的帶限制近似，它通過有限差分近似微分算子。

http://en.wikipedia.org/wiki/Screened_Poisson_equation

來源

2013-08-28 19:56:21 fchen

求解方程的帶狀矩陣系統

回答

相關問題