一： 关于二进制的读取
1)--------------------------------------------------------------------------------------
FOR的输出分 有格式form = 'formatted'、无格式form = 'unformatted'两种，前者是默认输出格式，即如果open语句里不声明form的话，那就是formatted。
无格式又分 直接存取access = 'direct'、间接存取access = 'seuqential' 两种，后者是无格式文件的默认格式。几个open示例：有格式文件：
open( 1, file = '...', status = 'new', form = 'formatted' )无格式/顺序文件：
open( 1, file = '...', status = 'new', form = 'unformatted', access = 'sequential' )无格式/直接文件： open( 1, file = '...', status = 'new', form = 'unformatted', access = 'direct', recl = ... )
2)--------------------------------------------------------------------------------------
如果是向 有格式文件 输出，则write语句写成：write( 1,* ) var； 或如，write( 1,'(5i2)' ) var
如果是向 无格式/顺序文件 输出，则write语句应写成： write( 1 ) var
如果是向 无格式/直接文件 输出，则write语句应写成： write( 1, rec = k ) var
3)--------------------------------------------------------------------------------------
对于 无格式/直接文件 ， open中的recl指明了一个输出记录的长度。它具体等于多少要看输出write语句如何写，同时write语句中的记录号rec也要写正确。
比如有一个m(20,30)的二维数组，可以有如下几种写法（注意recl, rec, m的变化）：
a) 把整个数组看作一个记录（1个记录，1个记录长度20*30）：
open ( 1, file = '...', status = 'new', form = 'unformatted', access = 'direct', recl = 20*30 )
write( 1, rec = 1 ) m
b) 把一行看作一个记录（20个记录，1个记录长度30）：
open ( 1, file = '...', status = 'new', form = 'unformatted', access = 'direct', recl = 30 )
do k = 1, 20
write( 1, rec = k ) ( m( k, p ), p = 1, 30 )
end do
c) 把一个数组元素看作一个记录（20*30个记录，1个记录长度1）：
open ( 1, file = '...', status = 'new', form = 'unformatted', access = 'direct', recl = 1 )
k = 0
do i = 1, 20
do j = 1, 30
k = k + 1
write( 1, rec = k ) m( i, j ) ！rec表示输出数据的位置
或者
write( 1, rec = (i-1)*30+j ) m( i, j )
end do
end do
注：二维数组是先存列再存行的，所以上述a)打印m元素的顺序与c)不同，而与下述d)相同
d) 把一个数组元素看作一个记录：
open ( 1, file = '...', status = 'new', form = 'unformatted', access = 'direct', recl = 1 )
k = 0
do j = 1, 30
do i = 1, 20
k = k + 1
write( 1, rec = k ) m( i, j )
或者
write( 1, rec = (j-1)*20+i ) m( i, j )
end do
end do
在一些Fortran编译器中，对recl解释成字长而非字节，因此还需要将recl乘以4才能正确读取。后面将会详细讲解gfortran和ifort在这点上的区别。
4)--------------------------------------------------------------------------------------
二进制存放的数据还有一个反位的问题。一般在UNIX机上是二进制数是高位结束，在WIN和LINUX中是低位结束（没记错的话）。如果读数时出现很大的怪数，比如1.e+35之类，或者执行时显示读取文件到末尾，或者是记录不够长，并且你又确定程序没错，那么多半是反位问题的原因。这时如果你用的编译器有反位选项，可以加上，比如ifort编译器可以写：
ifort -convert big_endian xxx.f90 或
ifort -convert little_endian xxx.f90
也可以直接在FOR程序的OPEN语句中说明，比如：
open( 1, file='...', form='unformatted', access='direct', recl=..., convert='big_endian' ) 或
open( 1, file='...', form='unformatted', access='direct', recl=..., convert='little_endian') big_endian 适用于在本地PC读取UNIX服务器上生成的二进制数据。little_endian适用于在UNIX服务器读取本地PC生成的二进制数据。
二： 不同编译器的差别
1. 程序举例
program test
implicit none
real(kind=4) :: sst(144,73),q(144,73)
integer :: i,j,recl1
recl1=144*73*4
sst(:,:)=10.0
q(:,:)=20
open(11,file='out.grd',form='unformatted',access='direct',status='replace',recl=recl1)
write(11,rec=1) ((sst(i,j),i=1,144),j=1,73)
write(11,rec=2) ((q(i,j),i=1,144),j=1,73)
close(11)
end program
2. 分析
kind=4表示单精度数据占用4个字节，一般默认也是4个字节； recl的设置，不同的编译器是不同的：
1> 在pgi和gfortran编译器下， recl1=144*73*4，生成的二进制文件大小是84096字节（84096=144*73*4*2）。
2> 在ifort编译器下，recl1=144*73*1, 生成的二进制文件大小同上。
因此说明不同编译器在读取（写入）二进制文件时有*4或不*4的问题，pgi、gfortran需要，而ifort不需要。
如果是双精度数据时，在pgi和gfortran编译器下，recl=144*73*8，而ifort编译器下recl=144*73*2。
如果用顺序（sequential）的方式读取二进制文件，文件开始和结束会自动加上文件标志符，而不同的编译器处理这点是有差别的，所以，当涉及到多个编译器时，建议用direct的方式读取二进制数据，因为可控性强。
另外，如果一次输出或读取的数据很大时，可能会超过recl（默认是kind=4）的数据区间，这时需要将数据分段输出或读取。
三：参考资源
http://bbs.lasg.ac.cn/bbs/thread-13889-1-9.html
http://bbs.06climate.com/forum.php?mod=viewthread&tid=14763
2014年5月5日星期一
2013年5月9日星期四
线性方程组求解与矩阵求逆
1. 线性方程组
在大地测量数据处理方面，线性方程组的求解无疑是最重要的问题之一。线性方程组的通用形式可写为\( AX=B \)，系数矩阵A的维数为 \( (M，N) \)，\( M \)为观测方程的个数（等于观测值的个数，1个观测值一个方程），\(N\)为需要求解的系数个数，因此向量\( X \) 的长度为\( N \)，观测值向量\( B \)的长度为\( M \)。根据\(M\)和\( N \)的大小观测，可以将线性方程组分为三类：
在大地测量数据处理方面，线性方程组的求解无疑是最重要的问题之一。线性方程组的通用形式可写为\( AX=B \)，系数矩阵A的维数为 \( (M，N) \)，\( M \)为观测方程的个数（等于观测值的个数，1个观测值一个方程），\(N\)为需要求解的系数个数，因此向量\( X \) 的长度为\( N \)，观测值向量\( B \)的长度为\( M \)。根据\(M\)和\( N \)的大小观测，可以将线性方程组分为三类：
- \( M > N \), 超定方程组（也称矛盾方程组），此时通常求方程的LS解；
- \( M = N \)，适定方程组，有唯一解；
- \( M < N \)，欠定方程组，有无数组解，但通常求其最小范数解。
接下来主要讨论超定方程组和适定方程组的情况。
对于超定方程组，根据LS准则，可得到法方程：\( A^TAX=A^TB \)，令\( C=A^TA，D=A^TB \)，则\( C \)的维数为\((N，N)\)，D的维数为\((N，1)\)，此时法方程的求解与适定方程组一致。下面的讨论以适定方程组为例。
2. 线性方程组的求解与求逆的关系
线性方程组\( AX=B \)的求解在数学上表示为\( X=A^{-1} \times B \)，但是在实际数值计算过程中，并不是一定要去计算\( A^{-1} \)，如常用的消去法，三角分解前推回代法等。
矩阵求逆等价于求解方程组\( AX=I \)，理论上与求解线性方程组一致。
在数值计算过程中，我们必须考虑计算的速度以及对内存的耗用情况，此时求逆和线性方程组的求解是两个不同的概念。
在求解线性方程组之前，必须明确你所需要的结果。如果只要最后结果\( X \)，则尽量不要去求逆；
如果系数阵\( A \)的求逆无法避免，尽量选择快速的求逆算法。当然，对于小矩阵，就不用考虑两者的差别了。
用matlab函数表示以上情况：
\( X=A\B \) (基于LU分解，前推回代直接求解参数向量\( X \))；
\( X=A^{-1} \cdot B \) （基于求逆求解参数向量\( X \) ）；
通常第一种方法速度快，占用内存少，但无法给出\( A^{-1} \)的信息。
3. 线性方程组的常用数值解法
通常把线性方程组的数值解法分为直接法和间接法。所谓直接法，就是经有限次运算即可求得方程组的精确解的方法；迭代法将求解的方程组转化为一个无限序列，其极限就是方程组的解。直接法一般来说用于求解系数矩阵阶数不是太高的问题，工作量较小，精度较高，但程序较复杂；迭代法主要用于一些高阶问题，一般程序较为简单，但工作量有时较大，耗用内存较多。
1> 直接法
a. 克莱姆法则
\( x_i=\frac{det(A_i)}{detA} \)
算法复杂度：\( (n+1) \times [(n-1) \times n! + n!] \)
运算量非常大，理论上可行，实际计算中很少用；
三角形方程组的求解简单快速，所以通常将系数阵A转化为三角阵再进行运算。不同的分解形式对应不同的分解方法。求解下三角阵L的过程，称为前推\( (X_1-X_2---X_n) \);求解上三角阵U的过程，称为回代\( (X_n---X_2-X_1)\)。前推和回代的计算量一样，为\( \frac{n^2}{2} \)次乘法和加减法以及n次除法。
b. 高斯消去法
利用消去法，将方程组变换为上三角形方程组；整个消去过程需完成\( O(\frac{1}{3}n^3) \)次乘法和加法运算； 回代过程需\( \frac{n^2}{2} \)次乘法和加减法以及n次除法，故用消去法求解n阶方程组的计算量约为\( O(\frac{2}{3}n^3+2n^2) \)；
从高斯消去法进一步完善，产生了列主元和全主元消去法，这两种方法在消去过程中的计算量进一步增加，但方法上更严谨。
每一步消去过程对应系数阵A左乘了一个初等变换矩阵，这一系列初等变换矩阵的乘积对应一个下三角矩阵L，故高斯消去法理论上等价与LU分解法。
c. LU分解法
将系数阵A直接分解成两个三角阵的乘积，根据矩阵运算的规则计算每个三角阵的元素。分解后\( AX=B \)的计算转化为两个三角形方程组的求解：
\( LY=B \) （前推，计算出Y）
\( UX=Y \) （回代，计算出X）
算法复杂度为 \( O(\frac{2}{3}n^3+2n^2) \)，所以直接分解法与消去法的计算工作量基本上是相同的。
求解逆阵的算法复杂度为\( O(2n^3) \)，是求解线性方程组的3倍。
d. 正交分解法
将矩阵A分解为一个正交阵和一个上三角阵的乘积，\( A=QR \)。根据分解方法，通常有Householder方法（镜像映射法）和Schmidt正交化方法。
运算复杂度为：\( O(\frac{4}{3}n^3+2n^2) \)，与上面两种方法相比，运算量多一倍。
该方法在求解超定方程组时，较为优越，原因如下：
\( A^TAX=A^TB \) with \( A=QR \), then \( R^TQ^TQRX=R^TQ^TB \) and finally, \( RX=Q^TB \)，因此，QR分解大大减少了计算量。
e. 对称正定矩阵的平方根法和\( LDL^T\)分解法
对称正定矩阵是特殊矩阵，可以用平方根法和\( LDL^T\)分解法（通常也叫做Cholesky decomposition），其运算量和存储量较普通消去法和LU分解法均节省一半，为\( O(\frac{1}{3}n^3+2n^2) \)。
平方根分解：\( A=LL^T \)
\( LDL^T\)分解：\( L \) 为单位下三角阵
2> 迭代法
对于阶数不高的线性方程组，直接法有效，但对于阶数很高的矩阵，有时需要迭代法来解决。迭代的规则不同，迭代法也就不同。通常分为：线性迭代和非线性迭代；一阶迭代和P阶迭代；定常迭代和非定常迭代。具体，可参考教科书，如“数值计算方法（冯康）”。
共轭梯度法：非线性迭代法，适用于对称正定阵，不需选取任何迭代参数，使用方便。
此外，对于高阶方程组的处理通常还采用分块矩阵处理方法。节省内存需要量，但增加了计算量，需要更多的计算时间（由于内外存交换数据需要时间）。
4. 矩阵求逆
除了必须要求矩阵的逆，否则尽量不采用求逆阵的方法来求解线性方程组。
求逆阵的算法复杂度为\( O(2n^3) \)，大约是同等方法求线性方程组\( O(\frac{2}{3}n^3) \)的3倍。
矩阵乘积的算法复杂度和矩阵求逆是等价的。
矩阵求逆当前最快方法的复杂度为\( O(n^{2.37}) \)
高斯消去法，LU分解法，SVD等，复杂度都为\( O(n^3) \)，差别在于前面的系数。
特殊矩阵的处理通常会有特殊的方法，能大大减少运算量和存储量。
\( A^TAX=A^TB \) with \( A=QR \), then \( R^TQ^TQRX=R^TQ^TB \) and finally, \( RX=Q^TB \)，因此，QR分解大大减少了计算量。
e. 对称正定矩阵的平方根法和\( LDL^T\)分解法
对称正定矩阵是特殊矩阵，可以用平方根法和\( LDL^T\)分解法（通常也叫做Cholesky decomposition），其运算量和存储量较普通消去法和LU分解法均节省一半，为\( O(\frac{1}{3}n^3+2n^2) \)。
平方根分解：\( A=LL^T \)
\( LDL^T\)分解：\( L \) 为单位下三角阵
2> 迭代法
对于阶数不高的线性方程组，直接法有效，但对于阶数很高的矩阵，有时需要迭代法来解决。迭代的规则不同，迭代法也就不同。通常分为：线性迭代和非线性迭代；一阶迭代和P阶迭代；定常迭代和非定常迭代。具体，可参考教科书，如“数值计算方法（冯康）”。
共轭梯度法：非线性迭代法，适用于对称正定阵，不需选取任何迭代参数，使用方便。
此外，对于高阶方程组的处理通常还采用分块矩阵处理方法。节省内存需要量，但增加了计算量，需要更多的计算时间（由于内外存交换数据需要时间）。
4. 矩阵求逆
除了必须要求矩阵的逆，否则尽量不采用求逆阵的方法来求解线性方程组。
求逆阵的算法复杂度为\( O(2n^3) \)，大约是同等方法求线性方程组\( O(\frac{2}{3}n^3) \)的3倍。
矩阵乘积的算法复杂度和矩阵求逆是等价的。
矩阵求逆当前最快方法的复杂度为\( O(n^{2.37}) \)
高斯消去法，LU分解法，SVD等，复杂度都为\( O(n^3) \)，差别在于前面的系数。
特殊矩阵的处理通常会有特殊的方法，能大大减少运算量和存储量。
2013年5月7日星期二
lapack 线性方程组最小二乘求解函数dgels
1. DGELS函数简介
DGELS函数是利用最小二乘求解线性方程组（超定或欠定）AX=B或A'X=B的函数，处于driver level。利用LU分解或者QR分解进行求解。
2. DGELS函数调用格式
DGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
输入参数：
TRANS（CHAR）： 值为'N'，表示No transpose, 对应AX=B； 值为'T',表示Transpose，对应A'X=B;
M,N（INT）: 系数矩阵A的维数，M>N,超定（下面只以超定方程组为例）,M<N，欠定；
NRHS（INT）： 观测值向量B的列数（通常为1）；
A（Double）：系数阵；
LDA（INT）：同DGESV函数中的LDA，LDA>=max(1,M)，通常等于M；
B（Double）：观测值向量，当TRANS=‘N’时，B（M，NRHS），当TRANS=‘T’时，B（N，NRHS）；
LDB（INT）：同DGESV函数中的LDB，LDB>=max(1,M)，通常等于M；
LWORK（INT）：Length of Work, 数组WORK的长度，见http://www.netlib.org/lapack/lug/node118.html；
输出参数：
B（Double）：B(1:N)为计算的结果；B（N+1：M）误差信息，具体见函数说明；
A（Double）：M>N, QR分解矩阵信息；M<N, LU分解矩阵信息；
WORK（Double）：维数为LWORK，存储空间(WORKspace)信息，具体见说明http://www.netlib.org/lapack/lug/node120.html
INFO：=0，成功执行；=-i，第i个参数illegal；=i，A的第i个对角元素为0，A非满秩阵，无LS解。
3. DGELS函数举例：
4. 同类型函数
- DGELSY：利用完全正交分解计算AX=B的LS解
- DGELSS：利用奇异值分解计算AX=B的LS解
- DGELSD：利用奇异值分解（divide and conquer）计算AX=B的LS解
5. 参考资料
2013年5月6日星期一
lapack 线性方程组求解函数dgesv介绍
1. DGESV函数简介
DGESV函数是lapack函数库求解线性方程组AX=B的函数，处于drive level。利用LU分解进行求解，A=PLU，P为置换矩阵，L为单位下三角矩阵，U为上三角矩阵。
2. dgesv函数参数介绍
a) 形式
DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
输入参数：
N（INT）： 线性方程组的数目，矩阵A的维数；
NRHS（INT）： 矩阵B的列数（一般为1，单列）；
A（Double）： 系数矩阵A，维数（LDA，N），通常LDA=N；
LDA（INT）：Leading Dimension of A, LDA>=max(1,N);
B（Double）： 矩阵B，维数（LDB，NRHS），通常LDB=N，NRHS=1；
LDB（INT）： Leading Dimension of B，LDB>=max(1,N);
输出结果：
A（Double）：存储LU分解的L和U矩阵，L为单位下三角阵且对角线的单位值1没有存储；
IPIV（INT）：维数为N，利用该矩阵的值可以推求置换矩阵P，starting from a unit matrix E, E矩阵的第i行与E矩阵的第IPIV（i）行进行交换。
如IPIV=[2 2 3 4]，则P矩阵的结果为：
.gif)
B（Double）：存储最终的求解结果X值；
3. dgesv函数举例
http://www.nag.co.uk/lapack-ex/node5.html
最终返回结果：
B=X= [ 1.0000 -1.0000 3.0000 -5.0000];
A=PLU=[
5.2500 -2.9500 -0.9500 -3.8000
0.3429 3.8914 2.3757 0.4129
0.3010 -0.4631 -1.5139 0.2948
-0.2114 -0.3299 0.0047 0.1314]
4. dgesvx函数
dgesvx是expert dgesv函数，该函数可以返回更多的中间结果，如计算矩阵A的条件数，估计解的误差范围等等，相应的，函数也增加了很多参数，所以占用的内存也就更多。具体使用可以参见user book http://www.netlib.org/lapack/lug/node26.html和应用举例 http://www.nag.co.uk/lapack-ex/node6.html#DGESVX
DGESV函数是lapack函数库求解线性方程组AX=B的函数，处于drive level。利用LU分解进行求解，A=PLU，P为置换矩阵，L为单位下三角矩阵，U为上三角矩阵。
2. dgesv函数参数介绍
a) 形式
DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
输入参数：
N（INT）： 线性方程组的数目，矩阵A的维数；
NRHS（INT）： 矩阵B的列数（一般为1，单列）；
A（Double）： 系数矩阵A，维数（LDA，N），通常LDA=N；
LDA（INT）：Leading Dimension of A, LDA>=max(1,N);
B（Double）： 矩阵B，维数（LDB，NRHS），通常LDB=N，NRHS=1；
LDB（INT）： Leading Dimension of B，LDB>=max(1,N);
输出结果：
A（Double）：存储LU分解的L和U矩阵，L为单位下三角阵且对角线的单位值1没有存储；
IPIV（INT）：维数为N，利用该矩阵的值可以推求置换矩阵P，starting from a unit matrix E, E矩阵的第i行与E矩阵的第IPIV（i）行进行交换。
如IPIV=[2 2 3 4]，则P矩阵的结果为：
.gif)
INFO（INT）：求解状态参数。INFO=0，成功执行；INFO<0，则INFO=-i，第i个参数是illegal value；INFO>0, U(i,i)=0, LU分解成功执行，但是U是奇异的，最终结果无法计算；
DGESV函数调用了Lapack Computation level的两个函数DGETRF和DGETRS，前者是进行LU分解的函数，后者是利用分解的结果进行前向和后向推求最终结果X。简单介绍下这两个函数的调用格式：
DGETRF（M，N，A，LDA，IPIV，INFO），矩阵A的维度为（M，N），其它参数同上，返回结果为A（PLU），IPIV以及INFO。
DGETRS（TRANS，N，NRHS，A，LDA，IPIV，B，LDB，INFO），求解AX=B或A'X=B，该函数只有在INFO=0时才可以进行。输入变量'TRANS'为字符型，其意义如下：'N'对应AX=B； ‘T’对应A'X=B; 'C' 对应A'X=B (conjugate transpose)，其它变量的意义同DGESV函数。
3. dgesv函数举例
http://www.nag.co.uk/lapack-ex/node5.html
最终返回结果：
B=X= [ 1.0000 -1.0000 3.0000 -5.0000];
A=PLU=[
5.2500 -2.9500 -0.9500 -3.8000
0.3429 3.8914 2.3757 0.4129
0.3010 -0.4631 -1.5139 0.2948
-0.2114 -0.3299 0.0047 0.1314]
IPIV=[2 2 3 4]; 从一个单位阵开始，单位阵的第一行和第2行交换，第二行和第2行交换，第三行和第3行交换，第四行和第4行交换，
所以：
P=[0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];
L=[1 0 0 0; 0.3429 1 0 0; 0.3010 -0.4631 1 0; -0.2114 -.3299 0.0047 1];
U=[5.2500 -2.9500 -0.9500 -3.8000; 0 3.8914 2.3757 0.4129; 0 0 -1.5139 0.2948; 0 0 0 0.1314]经过验证PLU=A；
4. dgesvx函数
dgesvx是expert dgesv函数，该函数可以返回更多的中间结果，如计算矩阵A的条件数，估计解的误差范围等等，相应的，函数也增加了很多参数，所以占用的内存也就更多。具体使用可以参见user book http://www.netlib.org/lapack/lug/node26.html和应用举例 http://www.nag.co.uk/lapack-ex/node6.html#DGESVX
2013年5月4日星期六
fortran编译器
Windows和Linux下Fortran编译器众多，这里简要介绍一些常用的编译器。
一. Windows平台下
在Win平台下，用的较多的是基于微软Visual Studio而集成的一整套IDE、商业编辑器：CVF和IVF。
1. CVF（Compaq Visual Fortran）
最常见的是CVF6.5，基于VS6，是win98时代的东西，更新到V6.6c被Intel收购了，即后来的IVF。CVF只支持到Fortran 95标准。
2. IVF（Intel Visual Fortran）
CVF的后继者，支持扩展的2003标准。IVF对Fortran77，90，95的兼容性不错。
二. Linux平台下
Linux下的fortran编译器很丰富，而且几乎都是免费的，如g77, gfortran, ifort 以及g95。
g77以前是GNU里的成员，后来被gfortran替代了。g95是在gfortan的发展过程中独立出去的，现在就只有作者一个人在维护。ifort是Intel的产品，对于单独的用户是免费的。这四个编译器用的较多的是gfortran和ifort，关于这两个编译器的讨论也集中在对老版fortran代码的兼容性（容错性）以及编译后的程序运行效率问题。总体来说，ifort的兼容性更好些，效率问题不太清晰。以下下面这两篇帖子对这个问题做过一些讨论，希望读完后自己能有一个总体的评判。
一. Windows平台下
在Win平台下，用的较多的是基于微软Visual Studio而集成的一整套IDE、商业编辑器：CVF和IVF。
1. CVF（Compaq Visual Fortran）
最常见的是CVF6.5，基于VS6，是win98时代的东西，更新到V6.6c被Intel收购了，即后来的IVF。CVF只支持到Fortran 95标准。
2. IVF（Intel Visual Fortran）
CVF的后继者，支持扩展的2003标准。IVF对Fortran77，90，95的兼容性不错。
二. Linux平台下
Linux下的fortran编译器很丰富，而且几乎都是免费的，如g77, gfortran, ifort 以及g95。
g77以前是GNU里的成员，后来被gfortran替代了。g95是在gfortan的发展过程中独立出去的，现在就只有作者一个人在维护。ifort是Intel的产品，对于单独的用户是免费的。这四个编译器用的较多的是gfortran和ifort，关于这两个编译器的讨论也集中在对老版fortran代码的兼容性（容错性）以及编译后的程序运行效率问题。总体来说，ifort的兼容性更好些，效率问题不太清晰。以下下面这两篇帖子对这个问题做过一些讨论，希望读完后自己能有一个总体的评判。
注：关于ifort在ubuntu下的安装可以参加另一篇博文：
2013年3月25日星期一
ubuntu下安装Lapack
1. 关于Lapack和Blas
Lapack (Linear Algebra Package)是一个以Fortran90编写，用于数值计算的函数集。例如，用于解线性方程组、计算特征向量、计算矩阵QR分解、奇异值分解等问题。
Blas(Basic Linear Algebra Subprograms)，基础线性代数程序集，用于基本的矢量和矩阵运算。最初发布于1979年，并用于建立更大的数值程序包（Lapack）。一级Blas负责矢量-矢量运算，二级Blas负责矢量-矩阵运算，而三级Blas负责矩阵-矩阵运算。由于Blas的高效，在高质量线性代数软件中常常使用它。
2. Lapack安装
1> 下载Lapack
http://www.netlib.org/lapack/ 最新版本3.4.2
2> 解压
tar -zxf lapack-3.4.2.tgz
3> 编辑make.inc文件
进入解压后的文件夹，会有一个make.inc.example文件，将该文件重命名为make.inc文件；并查看该文件，默认编译器为gfortran，根据需要修改对应的参数。在ubuntu下，采用gfortran作为编译器，一般不许要修改该文件。
4> 编辑Makefile文件
将下列语句（设定生成的库，默认不包含blaslib）
lib: lapacklib tmglib
# lib: blaslib variants lapacklib tmglib
改成：
# lib: lapacklib tmglib
lib: blaslib variants lapacklib tmglib
生成的库中包含了blaslib。
5> make
安装并测试Lapack。
6> 将生成的静态库(.a)移到/usr/lib文件夹下
生成的静态库包含libblas.a liblapack.a libtmglib等三个库。
sudo mv *.a /usr/lib
3. Lapack使用
在编译程序时，语句如下（其中test.f90文件中包含了lapack里的routines）：
gfortran -o test test.f90 -llapack -lblas
4. Lapack概述
1> Levels of Routines
Lapack的程序分三个层次，分别是driver， computational和auxiliary。driver routines解决一个具体的问题，如求解线性方程组或计算一个实对称阵的特征值。如果用户的问题刚好可以用该层次的程序解决，推荐用该层次的程序。这类程序通常调用几个computational routines或者一些其他的子程序。computational routines解决单独的一个计算任务，如LU分解。用户可以根据自身的需求调用该类程序。auxiliary routines这里不做阐述。一般用户一般调用driver和computational routines较多。
2> 数据类型和函数命名规则
Lapack函数命名格式为XYYZZZ，其中：
X表示数据类型：
补充：
Lapack是用fortran写的，Lapacke是它的C语言接口，先安装Lapack,然后安装Lapacke,那么在程序中可以直接调用C函数来实现需要的功能。
Lapack (Linear Algebra Package)是一个以Fortran90编写，用于数值计算的函数集。例如，用于解线性方程组、计算特征向量、计算矩阵QR分解、奇异值分解等问题。
Blas(Basic Linear Algebra Subprograms)，基础线性代数程序集，用于基本的矢量和矩阵运算。最初发布于1979年，并用于建立更大的数值程序包（Lapack）。一级Blas负责矢量-矢量运算，二级Blas负责矢量-矩阵运算，而三级Blas负责矩阵-矩阵运算。由于Blas的高效，在高质量线性代数软件中常常使用它。
2. Lapack安装
1> 下载Lapack
http://www.netlib.org/lapack/ 最新版本3.4.2
2> 解压
tar -zxf lapack-3.4.2.tgz
3> 编辑make.inc文件
进入解压后的文件夹，会有一个make.inc.example文件，将该文件重命名为make.inc文件；并查看该文件，默认编译器为gfortran，根据需要修改对应的参数。在ubuntu下，采用gfortran作为编译器，一般不许要修改该文件。
4> 编辑Makefile文件
将下列语句（设定生成的库，默认不包含blaslib）
lib: lapacklib tmglib
# lib: blaslib variants lapacklib tmglib
改成：
# lib: lapacklib tmglib
lib: blaslib variants lapacklib tmglib
生成的库中包含了blaslib。
5> make
安装并测试Lapack。
6> 将生成的静态库(.a)移到/usr/lib文件夹下
生成的静态库包含libblas.a liblapack.a libtmglib等三个库。
sudo mv *.a /usr/lib
3. Lapack使用
在编译程序时，语句如下（其中test.f90文件中包含了lapack里的routines）：
gfortran -o test test.f90 -llapack -lblas
4. Lapack概述
1> Levels of Routines
Lapack的程序分三个层次，分别是driver， computational和auxiliary。driver routines解决一个具体的问题，如求解线性方程组或计算一个实对称阵的特征值。如果用户的问题刚好可以用该层次的程序解决，推荐用该层次的程序。这类程序通常调用几个computational routines或者一些其他的子程序。computational routines解决单独的一个计算任务，如LU分解。用户可以根据自身的需求调用该类程序。auxiliary routines这里不做阐述。一般用户一般调用driver和computational routines较多。
2> 数据类型和函数命名规则
Lapack函数命名格式为XYYZZZ，其中：
X表示数据类型：
- S 表示单精度；
- D 表示双精度；
- C 表示复数；
- Z 表示双精度复数；
YY表示矩阵类型：
- DI 表示对角阵；
- GE 表示一般矩阵；
- GT 表示一般三角阵；
- OR 表示正交阵；
- PO 表示正定阵；
- TR表示三角阵；
ZZZ表示具体的计算问题：
- SV或SVX（求解线性方程组），前者是simple版，后者是expert版（输出一些中间结果，耗费内存更多）
- QR（QR分解）
- LU （LU分解）
- ...
3> 几个常用的函数名
更多内容参见Lapack User's Guide, http://www.netlib.org/lapack/lug/
- dgesv 求解线性方程组AX=B
- dgels 线性方程组AX=B的LS解
- xyyTRF 将矩阵yy进行三角分解（根据矩阵yy的类型，选择不同的分解方法，如GE矩阵用LU分解，对称正定矩阵采用Cholesky分解）
- xyyTRS 基于xyyTRF分解的结果求解AX=B
- xyyTRI 基于xyyTRF分解的结果求解$A^{-1}$
因此，dgesv等价与dgetrf和dgetrs的组合，dposv等价与dpotrf和dpotrs的组合，dgels等价与dgeqrf和dtrtrs的组合。
更多内容参见Lapack User's Guide, http://www.netlib.org/lapack/lug/
5. 不错的资源
补充：
Lapack是用fortran写的，Lapacke是它的C语言接口，先安装Lapack,然后安装Lapacke,那么在程序中可以直接调用C函数来实现需要的功能。
2013年1月11日星期五
fortran中字符和数字之间的相互转换
fortran中自带的ichar和char函数只能处理单个字符和数字之间的转化，在实际应用中经常会遇到多个字符和多位数字之间的转换，如：
- 字符串'500' 转换为整型500；
- 整型500转换为字符型'500'；
- 浮点型500.00转换为字符型；
1. 字符型转换为整型
例：'500'----->500
character(len=30) :: a
integer b
a='500'
read(a,*)b
write(*,'(I3)')b
2. 字符型转化为浮点型
例：'500.00'---->500.00
character(len=30) :: a
real b
a='500.00'
read(a,*)b
write(*,'(F6.2)')b
3. 整型转换为字符型
例：500----->'500'
character(len=3) :: a
integer b
b=500
write(a,'(I3)')b
write(*,*)a
4. 浮点型转换为字符型
例：500.13----->'500.13'
character(len=6) :: a
real b
b=500.13
write(a,'(F6.2)')b
write(*,*)a
2012年10月15日星期一
gfortran 编写动态库时遇到 ‘relocation R_X86_64_32 against’错误
如题，用gfortran编写动态链接库时遇到 ‘relocation R_X86_64_32 against ...’错误，做如下处理即可：
在生成.o文件时，用：gfortran -c -fPIC *.f90
在生成.o文件时，用：gfortran -c -fPIC *.f90
在生成.so文件时，用：gfortran -shared -fPIC -o libxxxx.so *.o
2012年7月11日星期三
ubuntu下安装ifort xe
Intel Fortan Compiler简称ifort, windows下的ifort是收费的，但是linux系统下提供免费的ifort,可以在下面的链接中下载需要的版本（必须先注册，随后会收到官网发来的邮件，里面提供了接下来安装需要的series-number）http://software.intel.com/en-us/articles/non-commercial-software-download/
在安装ifort之前，需要先安装一些软件包
sudo apt-get install build-essential
sudo apt-get install g++
sudo apt-get install gcc-multilib
sudo apt-get install rpm
sudo apt-get install openjdk-6-jre-headless
在安装ifort之前，需要先安装一些软件包
sudo apt-get install build-essential
sudo apt-get install g++
sudo apt-get install gcc-multilib
sudo apt-get install rpm
sudo apt-get install openjdk-6-jre-headless
sudo apt-get install libstdc++6
将下载下来的文件解压 tar -zxvf l...tgz
安装
进入上面解压过后的文件夹，./install.sh
安装一共分六步，根据自己的需要选择设置信息，一般一路enter下去即可。
修改运行环境信息
将source /opt/intel/bin/ifortvars.sh ia32 （或者根据安装最后一步的提示将红色部分替换掉）添加到 ～/.bashrc文件里
检测安装是否成功
2012年1月3日星期二
fortran 编写动态库 linux系统
基于linux系统，编写fortran库函数过程如下：
以aa.f90, bb.f90, main.f90三个文件为例
1. 生成.o文件
gfortran -c aa.f90 bb.f90 将aa.f90和bb.f90文件生成.o文件
2. 生成动态库文件
gfortran -shared -fPIC -o libtest.so *.o 基于.o文件生成动态库文件(.so文件)
3. 将动态库文件移动到/usr/lib目录下
sudo mv libtest.so /usr/lib
4. 调用动态库文件
gfortran -c main.f90 将main.f90文件生成.o文件
gfortran -o main main.o -ltest 生成可执行文件main，其中-ltest为调用库文件libtest.so
以aa.f90, bb.f90, main.f90三个文件为例
1. 生成.o文件
gfortran -c aa.f90 bb.f90 将aa.f90和bb.f90文件生成.o文件
2. 生成动态库文件
gfortran -shared -fPIC -o libtest.so *.o 基于.o文件生成动态库文件(.so文件)
3. 将动态库文件移动到/usr/lib目录下
sudo mv libtest.so /usr/lib
4. 调用动态库文件
gfortran -c main.f90 将main.f90文件生成.o文件
gfortran -o main main.o -ltest 生成可执行文件main，其中-ltest为调用库文件libtest.so
Fortran 库函数（持续更新中）
1. FORTRAN调用shell函数
I=SYSTEM(‘command’)
例：I=SYSTEM(‘ls *.txt > filelist’)将当前目录下.txt文件名保存到filelist
2. FORTRAN取整函数
AINT(x[,kind]) 对x取整，并转换为实数
ANINT(x[,kind]) 对x四舍五入取整，并转换为实数
CEILING(x) 求大于等于x的最小整数
FLOOR(x) 求小于等于x的最大整数
IFIX(x) 将x转换为整数
INT(x) 将x转换为整型
例： a=1.56 b=-2.5
AINT对应值 1.000 -2.000
ANINT对应值 2.000 -3.000
CEILING对应值 2 -2
FLOOR对应值 1 -3
IFIX对应值 1 -2
INT对应值 1 -2
3. 计算时间统计
call CPU_TIME(t1) 放在开始计算前
call CPU_TIME(t2) 放在计算结束后
t=t2-t1
I=SYSTEM(‘command’)
例：I=SYSTEM(‘ls *.txt > filelist’)将当前目录下.txt文件名保存到filelist
2. FORTRAN取整函数
AINT(x[,kind]) 对x取整，并转换为实数
ANINT(x[,kind]) 对x四舍五入取整，并转换为实数
CEILING(x) 求大于等于x的最小整数
FLOOR(x) 求小于等于x的最大整数
IFIX(x) 将x转换为整数
INT(x) 将x转换为整型
例： a=1.56 b=-2.5
AINT对应值 1.000 -2.000
ANINT对应值 2.000 -3.000
CEILING对应值 2 -2
FLOOR对应值 1 -3
IFIX对应值 1 -2
INT对应值 1 -2
3. 计算时间统计
call CPU_TIME(t1) 放在开始计算前
call CPU_TIME(t2) 放在计算结束后
t=t2-t1
4. 统计函数
(1) 计算最大值，最小值
MAX(v1,v2,...)找出v1,v2...中的最大值，此函数没法判断一个数组的最大值
MAXVAL(A,dim) 根据dim值找出一个数组中每一维的最大值
MAXVAL(A) 返回整个数组的最大值
MAXLOC(A,dim) 根据dim值找出一个数组中每一维最大值的位置
同理：MINVAL,MINLOC相应的范围数组的最小值和最小值的位置
例： 数组A（5,2），值为/1.2,3.5,4.6,8.7,9.6
2.3,63,3.9,6.8,1.6/
MAXVAL(A) 结果为63
MAXVAL(A,1) 结果为9.6 63
MAXLOC(A,1) 结果为5 2
(2) 求和函数
SUM(A,dim)根据dim求出每一dim数组元素的和
例： 数组A（5,2），值为/2.3,3.6,4.7,5.8,6.9
1.0,1.3,3.6,7.3,4.5/
SUM(A) 结果为41.00
SUM(A,2) 结果为3.3 4.9 8.3 13.1 11.4
SUM(A,1) 结果为23.3 17.7
订阅：
博文 (Atom)