2013-04-05 47 views
1

我在FORTRAN 77中編寫了一個模擬程序,其中我需要從某些實驗數據中獲取參數值。數據來自互聯網數據庫,所以我已經預先下載了它,但是沒有簡單的數學模型可以用來連續提供數值 - 我只有離散的數據點。但是,我需要知道這個參數是否在x軸上有任何值,而不僅僅是我從數據庫中得到的離散值。實現高性能FORTRAN代碼查找的智能方法

爲了簡化,你可以說,我知道的f(x)x所有整數值的數值,需要一種方法來找到f(x)任何實際x值(從不以外的最小或最大x我的知識)。

我的想法是獲取數據並進行線性插值,以便能夠獲取參數值;在僞代碼:

double xd = largest_data_x_lower_than(x) 

double slope = (f(xd+dx)-f(xd))/dx // dx is the distance between two x values 
double xtra = x-xd 

double fofx = f(xd)+slope*xtra 

要實現這一點,我需要對數據點進行某種查找。我可以通過從數據庫中獲取所有整數x的值來簡化xd的查找,以便xd = int(x)dx = 1,但我仍然不知道如何實現f(xd)的查找。

什麼是一個很好的實現方法?

在一次模擬運行期間,數值將被提取爲10^7到10^9次,因此性能至關重要。換句話說,每當我需要f(xd)的值時,從IO中讀取數據不是一種選擇。

我目前在每行中有一對(製表符分隔)x,f(x)的文本文件中的數據點,因此解決方案的獎勵點也提供了從那裏獲取數據的平滑方式,需要是。

回答

2

你說你有所有整數的值。你是否有i, f(i)所有整數iMN?然後將值f(i)讀取到數組y尺寸爲M:N。除非值的數量是巨大的。對於MN之間的實際值,很容易將其索引到數組中並在最接近的一對值之間進行插值。

爲什麼使用FORTRAN 77? Fortran 90/95/2003已經和我們在一起幾年了...

編輯:在評論中回答問題,重新如何在FORTRAN 77中只讀取一次數據值,而不必將它們作爲在一個長長的呼叫鏈中的爭論。技術1:在程序啓動時,將它們讀入數組中,該數組位於已命名的公共塊中。技巧2:第一次調用返回f(x)的函數時,將值讀入本地變量中,該變量也在SAVE語句中。使用SAVEd的邏輯來指定函數是否在第一次調用時。一般來說,我更喜歡技術2更「本地化」,但它不是線程安全的。如果你是在做並行模擬,第一種技術可以在程序進入多線程之前的啓動階段完成。

以下是使用SAVE的一個示例:。 (在Fortran 95符號...轉換爲FORTRAN 77)。將讀取的數據放入IF塊中的數組中。

+0

FORTRAN 77由我正在工作的團隊強加於我。這是一個臨時職位,所以我決定迫使團隊繼續前行的好處不僅僅是學習:P如何將值讀入內存*只有一次*,而不必一直傳遞數組?這種查找是相當順利的 - 如果可能的話,我想避免在四個或五個地方添加一些額外的函數參數... – 2013-04-06 01:21:32

+0

它是並行的,使用MPI - 換句話說,我運行該程序在一些單獨的進程中(但我不在進程內執行任何線程分支)。我認爲這個通用模塊將是最好的解決方案,以確保我能夠以安全的方式讀取文件 - 明天我會試試看,看看我能否實現它。好的; – 2013-04-07 21:19:31

+0

好的;我現在已經決定採用一種可以在我的並行環境中安全線程工作的方法。我將讀入文件,並在主線程中創建一個足夠精細的x分辨率插值,然後將該線程廣播到所有其他線程。該向量將駐留在一個公共塊中,以避免將其作爲參數傳遞。如果遇到問題,這可能與方法的選擇無關,而是對我來說是一個笨拙的程序員,所以我現在會認爲這已經解決了。謝謝您的幫助! – 2013-04-08 16:42:48

0

您可能想要一種插值或適合數據的方法,但您需要更具體地說明數據的維度,數據的行爲方式,訪問數據的方式(例如,也許您的下一個請求總是接近最後一個請求),網格的製作方式(均勻間隔,隨機或其他方式),以及您需要的數據能夠知道哪種方法最適合您。

但是,如果現有數據集非常密集且接近線性,那麼您當然可以進行線性插值。

+0

您是否將「線性插值」與「線性擬合」混淆?數據遠不是線性的,所以線性*擬合*不存在問題,但是如果我能找到一個很好的方法來查找最接近給定x值的兩個數據點,它們之間的線性插值就好了。 (我可以確保x-space中的數據足夠密集。)'x'和'f(x)'都是標量,所以具有所有x值和所有函數值的矩陣將是2-by -N,其中N是數據點的數量。 – 2013-04-05 23:05:42

0

使用您的數據庫(文件),您可以創建一個數組fvals,其中fvals(ii)是函數f(xmin + (ii-1) * dx)。 x值xx與您的陣列索引之間的映射爲ii = floor((xx - xmin)/dx) + 1。一旦知道ii,您可以使用它周圍的點進行插值:可以使用iiii+1或某些更高階多項式插值進行線性插值。對於後者,您可以使用來自Numerical Recipes的相應polint程序。請參閱第103頁。