Fortran筆記：
這裡會放一些以前上課時沒教到，但又蠻有趣的用法，未來累積的量如果太多，會再重新編排。
請愛用Ctrl+F搜尋關鍵字
- 存取資料長度不一的二進位檔
- 指標（Pointer）
- 輸出處理
- 直接於執行程式時傳入變數（./a.out 變數1 變數2）
- 對BIG/LITTLE Endian的處理
- 文字處理
- 透過Fortran執行系統指令
- 亂數
- 呼叫以C(或其他語言)編譯的Library
- Namelist的製作與使用
- 製作DLL檔
以前在學校學到用Fortran處理的二進位檔都是資料長度一致的陣列（比如說4-byte的實數），後來才慢慢接觸到包含不同資料類型的二進位檔（例如加了Header的資料，像是RIPv4的輸出，ESRI的Shapefile等等），以下是這陣子想到的處理方式：
a. 以循序檔的方式建立的檔案（如RIPv4的輸出）
當以sequential的方式寫入檔案時，每個write指令所輸出的變數前後各會被加上一個4-byte的header，以及一個4-byte的terminator，它們一對數值相同的unsigned integer，標出這對flag裡面所包含的資料總長度，若你寫入資料的程式碼寫成：
real(kind=8) :: double = 3.14159d0
integer(kind=4) :: int32 = 999
open(10,file="temp.dat",form="unformatted")
|
write(10) double
write(10) int32,"abc",double
|
! temp.dat之大小為39個byte
! 這裡的write因為是要寫入unformatted的二進位檔，所以不需要指定格式
! 若你直接用文字編輯器打開它，你會發現除了abc以外其他都是亂碼
! 因為一個字本身就是一個byte（不考慮特例），所以寫入時不用經過轉換
|
如果你嘗試將頭四個byte讀成整數，你會發現他的值是8，因為他後面放的雙精度時數大小就是8個byte；同理，將第17~20個位元組取出轉成整數你會得到15。
要讀取這類資料，因為header/terminator的關係，所以一個write要對應一個read（要用access="direct"的隨機存取來跳過也是可以，但是就是要知道它們是在第幾個byte的位置上）
| open(10,file="temp.dat",form="unformatted") |
read(10) double
read(10) int32,char3,double |
b. 資料長度已知，不是循序檔
假設我知道這筆資料的長度，也知道它是不是用sequential的方式寫出，那麼：
! 資料結構：int32,char3,double,int32,int32,int32，總長度4+3+8+4+4+4
open(10,file="temp.dat",form="unformatted",access="direct",recl=27)
read(10,rec=1) int32,char3,double,(int32(i),i=1,3)
|
這種寫法的壞處是比較沒有彈性，在recl指定了一組record的資料長度，故所有要讀取的變數都要塞在read的那一行，但好處則是因為是用隨機存取，所以可以快速調出想要的資料。
c. 資料長度未知，不是循序檔（如Shapefile的.shp檔）
其實也不一定是未知長度才能用這個方法，.shp檔用方法b其實也可以讀，但如果不想把要讀取的變數都塞在一行的話，用access="stream"會是比較好的選擇：
! 資料結構：int32,char3,double,int32, int32,int32，總長度4+3+8+4+4+4
open(10,file="temp.dat",form="unformatted",access="stream")
read(10) int32,char3
read(10) double,(int32(i),i=1,3) |
這裡你的read要拆成幾個放都可以，反正只要後面讀取的變數有對應到資料結構的話就ok。
如果配合POS的話，甚至可以像access="direct"一樣跳到要讀的位置（該位置以byte表示，不過POS這東西似乎是compiler-dependent，我用PGF90不能過），詳細用法可以參考底下fortranwiki這個連結。
參考資料：http://fortranwiki.org/fortran/show/Stream+Input+Output
http://mitgcm.org/download/daily_snapshot/MITgcm/doc/OutputFiles
在C語言中，指標是個很重要的功能，若想修改傳入function的變數就得利用指標來傳遞記憶體位置（傳址，by
reference），否則除了陣列之外都是以值來傳遞（傳值，by value），Fortran則是全部都用傳址的方式來傳遞變數，這功能
直到Fortran90以後才出現，在這之前想要使用此功能的話需要另外安裝CRAY pointer這個延伸套件，使用上有些要注意的事項：
a. point與target的屬性：
指標變數要加上pointer的屬性，而被指向的變數則要加上target的屬性。
b. 型別要一致：
real要對real，integer 要對integer，雙精度的實數其指標也要宣告為雙精度。
c. 維度要注意：
指標可以不用像陣列變數一樣宣告大小沒關係，但是若要指向多筆資料，可能要宣告其維度。
d. 指向符號為"=>"：
當指向一陣列變數時，可以只指向該陣列之某一範圍。
e. 取消指標與變數的連結
/ 查詢指標是否有指向任何變數：
使用nullify來取消指標與變數之間的連結；使用associated()來查詢指標是否有指向任何變數。
real(kind=8),target :: a(3,3) = reshape( [11,21,31,12,22,32,13,23,33] , [3,3])
real(kind=8),pointer
:: ptrA,ptrB(:),ptrC(:,:)
|
ptrA => a(2 , 2 )
ptrB => a( : , 2 )
ptrC => a(2:3 , : )
nullify(ptrA)
if(associated(ptrB)) print*,"ptrB is => sth."
|
註：若只是想讓一個變數擁有多個名稱的話，也可以用equivlence()這個指令，它們將會指向同一個記憶體區塊，詳細用法及限制請參考該語法之說明。
指標在C語言裡的另一個重要的功能就是用來做記憶體的動態配置，在Fortran90以前由於沒有allocate()/deallocate()可
以用，想要動態配置記憶體的話得安裝CRAY
Pointer延伸套件，之後就可以利用語法類似於C的malloc/free做記憶體的配置與釋放，在F90也是有支援這用法，我第一次看到這種寫法是
在Vis5D的GrADS格式轉換工具裡（1994年的程式碼）。具體的用法可以參考http://gcc.gnu.org/onlinedocs/gfortran/MALLOC.html
參考資料：http://www.personal.psu.edu/jhm/f90/lectures/42.html
a. 純文字之格式化輸出
在透過print或是write輸出文字時，除了方便的自由格式之外，另一個選擇就是指定輸出格式，在指定輸出格式時，需要指定字元數量，不過老實說
這樣有點麻煩，每次都要去計算字元數，比較簡便的方法就是不指定字元數量，直接用"a"就可以輸出所有字元了：
print*,"Paragraph 1-","paragraph 2"
print "(a15,a11)","Paragraph 1-","paragraph 2"
print "(a,a)","Paragraph 1-","Paragraph 2"
! 輸出結果
Paragraph 1-Paragraph 2
Paragraph 1-Paragraph 2
Paragraph 1-Paragraph 2 |
b. 輸出整數前補零
有時在輸出內容時，會為了排版好看或其他原因而需要在整數結果前面自動補零，讓不同位數的資料看起來一樣長，那麼可以在格式後加上欲「輸出之位數」.
「欲補零之數量」，而在指定格式時，跟例a一樣，標明位數與否都是可以的：
c. 10進位轉2進位/8進位/16進位之輸出（DEC to
BIN/OCT/HEX）
透過指定格式"b"可將10進位之整數轉為2進位之值，用"o"可以轉為8進位，"z"可轉換為16進位之值，若想讓結果前面自動補零，同樣也可以依
照上面對整數的作法來處理：
print
"(z)",5566
print "(z6.6)",5566
print "(o)",5566 |
d. 不自動換行
在C語言中，printf要配合"\n"才會有換行的作用，而在Fortran中write或print指令都有自動換行的功能，如果不希望輸出自動換行的話，只要在write裡加入advance="no"即可（此為90標準，同樣也可以應用在read上），另外也可以在輸出格式裡的最末端加上"$"，不過這似乎就不是標準了：
print
"(a$)", "Please enter your name:"
! 以下為advance=no的用法
write(*,*,advance="no") "Please enter your name:"
read(*,*) YourName
|
如此一來讓使用者輸入的地方就會跟在提示語句後面，在版面上會好看一點。
透過呼叫getarg或是這個函數，就可以將指令列後所接續的變數讀入，在這之前可以用iargc()來檢查使用者是不是有輸入正確數量的變數：
（新版的程式建議使用get_command_argument這個函數，該函數為Fortran
2003之標準）
character :: Var1,Var2*2
if ( iargc()
/= 2 ) then
print*,"ERROR : Two arguments required"
stop
endif
call getarg(
1, Var1)
call getarg(
2, Var2)
print*,"變數1為：",Var1,"變數2為：",Var2
|
參考資料：http://gcc.gnu.org/onlinedocs/gfortran/GETARG.html
http://gcc.gnu.org/onlinedocs/gfortran/GET_005fCOMMAND_005fARGUMENT.html
這東西老實說有點麻煩…如果你電腦記憶體設計是LITTLE_ENDIAN的配置，而想要讀取BIG_ENDIAN的資料（或反之），那麼大概有三種
方法可以來處理這個問題：
a. 利用CONVERT specifier
指定convert="BIG_ENDIAN"即可，此外還有SWAP,LITTLE_ENDIAN等標籤可使用。不過在我的vim裡convert
這個字並不會變色，據gilocustom表示，這可能是因為它還不在標準裡（所查閱的文件為90/95/2003/2008標準的草稿）
| open(10,file="Test.dat",form="unformatted",access="stream",convert="BIG_ENDIAN")
|
b. 在編譯時加上flag
利用這個方法，會自動對調Endianness
gfortran -fconvert
pgf90 -byteswapio
c. 透過副程式來轉換
在正常情況下，用到這方法的機會很少，但是如果一個資料裡面有BIG_ENDIAN的資料，也有LITTLE_ENDIAN的資料…那就只能用這個方
法了，Shapefile（.shp）就是這麼奇怪的檔案。
可以下載這隻由NCAR/CGD/CAS科學家David Stepaniak所寫的副程式來做轉換：
http://www.cgd.ucar.edu/cas/software/endian.html這程式
會對實數做轉換，原理其實就是先將資料由TRANSFER轉換至整數型別，再透過呼叫mvbits來做位元組的搬移，後面的搬到前面，前面的搬到後面，搞
定後再把型別轉換原本資料的型別，以下是我簡化過的Code：
! AUTHOR:David Stepaniak, NCAR/CGD/CAS
! DATE INITIATED: 29 April 2003
! LAST MODIFIED: 3 Augest 2012 by Cypresslin
function
rSwap(rIN)
implicit none
integer :: i_element, i_element_br
real,intent(in) :: rIN
real,intent(out):: rSwap
integer :: i
i_element = TRANSFER( rIN, 0 )
do i=1,4
call
mvbits(i_element,8*(4-i),8,i_element_br,8*(i-1))
enddo
rSwap = TRANSFER( i_element_br, 0.0 )
return
|
註：字串資料由於一個字元為一個byte，所以它不會受到byte-ordering的影響
參考資料：http://gcc.gnu.org/onlinedocs/gfortran/CONVERT-specifier.html
http://www.cgd.ucar.edu/cas/software/endian.html
a. 字串之連接
只要用","即可，但是逗號的用法在以逗號分隔變數的函數（如open），或是指定變數時中會出問題，此時只要用兩個跳脫符號"//"即可連接兩字
串：
print*,"str1"//"str2"
open(10,file="Header"//cFileName,status="old") |
a.數字/字串轉換
數字與字串間的轉換，沒有類似C那種atoi()等函數可以使用，必須靠read跟write來完成：
integer :: num
character :: txt*4 = "1000"
read(txt,"(i)") num
write(txt,"(i4)") num
|
b.英文大小寫轉換
此轉換可以靠iachar()函數得出該英文字母的ASCII代碼（http://www.asciitable.com/）再經過代碼的加（大寫轉小寫）或減（小寫轉大寫），
找到對應的字母代碼後，用achar()函數把它轉換成字母即可，以下示範如何將字串由小寫轉大寫：
character
:: txt*5="a1Bc$"
integer :: code
do i=1,len(txt)
code = iachar( txt(i:i)
)
if( code >= 97 .and. code <= 122 ) then
txt(i:i) =achar( code -32 )
endif
enddo |
註：iachar/achar與ichar/char其實做的事情一樣，但是前者是根據ASCII
Code去做轉換（故稱iAchar與Achar，A for
ASCII），而後者是用系統本身的字碼表去轉換，故考慮可移植性(portability)請使用iachar/achar。
其實在上例中的if判斷式裡可以直接透過比較運算子讓字元跟字元比較，程式會自動用內碼來比，而不需要事先轉換：
if(
txt(i:i) >= "a" .and. txt(i:i) <= "z" ) then
txt(i:i) = achar( iachar( txt(i:i) - 32 )
endif |
註：進行這種比較時，電腦會用系統本身的字碼表去處理，所以如果你的電腦使用的不是ASCII代碼而是EBCDIC，那結果將會不同。為了可移植性請用後
述之方法。
LLT/LLE/LGT/LGE這幾個字元比大小的函式（其大或小的定義不論在任何電腦上都是依照ASCII碼）它們的名稱其實就是在一般常用的比較
運算子前面加上L（lexically）而已，蠻好記的：
if(
LGE(txt(i:i),"a") .and. LLE( txt(i:i),"z") ) then
txt(i:i) = achar( ichar( txt(i:i) - 32 )
endif
|
除了用ASCII碼來做大小寫的轉換之外，也可以用index這個字串搜尋功能，配合預先寫好26個字母的字串來做代換，下例中假設txt的宣告跟值
都與先前相同：
character
:: a2zL*26="abcdefghijklmnopqrstuvwxyz"
character :: a2zU*26="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
integer :: num
do i=1, len(txt)
num = index( a2zL,txt(i:i) )
if ( num > 0 ) txt(i:i) = a2zU( num : num )
enddo
|
參考資料：http://gcc.gnu.org/onlinedocs/gfortran/LLT.html
http://gcc.gnu.org/onlinedocs/gfortran/IACHAR.html
在程式執行時，想要執行例如關機、新增資料夾等等之類的一些系統指令（Shell Commands）可以用call system()執行，而在Fortran 2008之後的標準有EXECUTE_COMMAND_LINE可以用。
call system("mkdir test")
call system("cd ..; ./a.out")
|
註：其實如果可以的話，盡量避免使用這個指令比較好，因為不同平台指令不同，所以會影響到程式的可攜性（當然如果沒打算移植程式碼的話就沒差）；還有就是可能會有權限上的問題。
參考資料：
http://gcc.gnu.org/onlinedocs/gfortran/SYSTEM.html
其實這裡所謂的亂數應該稱為「假亂數（pseudorandom number）」，因為電腦在產生亂數時其實是依據固定的演算法配合亂數種子去算出來的，故使用一樣的亂數種子就可以得出一樣的答案。
亂數的產生其實很簡單，只要call random_number()即可，若希望每次產生的亂數都不同，就得要在產生亂數前先產生亂數種子call random_seed()，如下：
real :: num(3),x
call random_seed()
! call init_random_seed()
call random_number(num)
call random_number(x)
print*,num,x
! 產生固定範圍(1-6)的亂數
print*,int(x*6)+1 |
random_number()裡面可以放陣列或是單一變數，以上面的例子來說，丟入陣列num那它就會將該陣列裡放滿亂數。
就如同前面所說的，在產生亂數前若不製作新的亂數種子的話，會因為在每次執行時都使用相同的亂數種子而得到一樣的答案，這樣其實有個好處就是可以確保每次都可以拿到一樣的值，利於對採用亂數做運算的程式的除錯。
這裡還有一個問題就是亂數種子的產生通常是用主機上的時間，所以在極短的時間內（例如一秒內，也就是主機系統時鐘變動前）執行多次亂數產生程式，會因為用來當種子的時間都一樣，所以得到一樣的結果。這種問題可以靠空迴圈或是sleep之類的函數強制讓程式暫停一下下，如此一來程式再次執行時就能夠抓到不同的系統時間，缺點就是會較為耗時。另外也許可以利用GUID的HASH Code當作亂數子（VB .NET與C#上有這種做法，但可惜在Fortran上我沒有實做過）。
若還有興趣研究，這裡有篇來自NUMERICAL RECIPES IN FORTRAN 77的文章可以參考。
參考資料：
http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fNUMBER.html
http://nuclear.fis.ucm.es/COMP-PHYS/RANDOM/RandomNumbers.pdf
有時候會拿到以其他語言寫成的副程式，然而想把它用Fortran改寫可能又不大合乎效益（而且是還要看得懂！），遇上這種情況直接呼叫它可能會省事點。這就是所謂的混合式程式設計（Mixed-Language Programming），關於這個主題大概可以寫一整面，就慢慢來吧！
假設我們現在有一個以C寫成的副程式：
#include<stdio.h>
void sub_(int *input){
printf("%d\n",*input);
}
|
這支簡單的副程式就只是讀入一個整數，然後把它印出來而已，用來呼叫它的主程式可以寫成這樣：
integer :: i=100
call sub(i)
end |
編譯時只要先把C寫成的副程式編譯成object file，再把主程式跟他一起編譯即可：
[root@HPC ~] cc -c sub.c
[root@HPC ~] pgf90 main.f90 sub.o
[root@HPC ~] ./a.out
100
|
要是在副程式的名字後面沒有加底線，編譯時應該會看到以下的錯誤訊息：
/tmp/ccwMUsZi.o(.text+0x19): In function `MAIN__':
: undefined reference to `sub_'
collect2: ld returned 1 exit status |
參考資料：https://engineering.purdue.edu/ECN/Support/KB/Docs/CallingCFromFortran
若想要讓程式裡的參數變動，大概有下列幾個方式：
- 參數寫在程式裡，壞處是每次要用就重新編譯一次。
- 利用read讀入使用者想設定的值，壞處是萬一可變動的參數很多，或者不是每次都會改到這麼多參數就有點麻煩（其實這可以用輸出輸入重新導向來解決）
- 使用namelist格式，讀入設定檔。
Namelist格式在F90後被列為標準，如果有在跑模式的話對這個東西應該不陌生，模式的大大小小設定幾乎都可以在 namelist 裡面調整。沒看過也沒關係，只要把它想成是一個程式的設定檔就好。
一個namelist檔的格式大概是這樣的：
&DOMAIN
SIZEX = 401,
SIZEZ = 121,
ds = 100.0000
/
&POISSON
ISSUCCESSIVE = F,
BONDTYP = 3,
ERROR = 4.00000003E-005
/
|
在namelist裡可以定義多個區塊，每個區塊用「&區塊名稱」開始，用「/」結束，在這區塊內就可以放置變數。這個檔案可以依照格式手動寫一份，或者是用程式幫你產生，如下：
implicit none
integer :: sizeX=401,sizeZ=121
real :: ds=100.0
logical :: isSuccessive=.false.
character :: BondTyp = "3"
real(kind=8) :: error=0.00004d0
! 建立namelist group，有點像是struct，一個名稱裡面跟著一些變數
namelist /DOMAIN/ sizeX,sizeZ,ds
namelist /POISSON/ isSuccessive,BondTyp,error
open(10,file="namelist.input")
! namelist group的寫入，只要在write裡面用nml指定該group的名稱即可
write(10,nml=DOMAIN)
write(10,nml=POISSON)
end |
要注意欲寫入之變數是不能使用parameter的。
而namelist的讀入只要依樣畫葫蘆把write改成read就好，不過要小心資料型別轉換問題，例如在上面那個程式中用字元寫出BondTyp，但是在寫出的namelist裡（如範例資料），BondTyp =3，這會導致程式認為3不是個字元而是整數，故讀入時要改用整數來讀，或修改資料他變成BondTyp= "3"。
integer :: BondTyp
! 其他宣告與寫入時差不多，只是沒有指定值而已，故略過
open(10,file="namelist.input")
read(10,nml=DOMAIN)
read(10,nml=POISSON)
end |
如果遇上entity name is not member of group.的錯誤，請檢察變數的數量與名稱是否皆與namelist中所設定的一致。
參考資料：http://www.lahey.com/docs/lfprohelp/F95ARNAMELISTStmt.htm
<-Memo
