linux如何编译fortran子程序,linux下Fortran编译Lapack、clapack库及使用的方法

apt安装：

$ sudo apt install libblas-dev liblapacke-dev liblapack-dev libtmglib-dev libatlas-dev libatlas-base-dev libblacs-mpi-dev libopenblas-base libopenblas-dev libf2c2-dev

一、获取lapack源代码

wget http://www.netlib.org/lapack/lapack-3.6.1.tgz

二、解压后编译

cd /mybk/lapack-3.6.1

cp make.inc.example make.inc

make后得到liblapack.a  librefblas.a  libtmglib.a

三、拷贝到系统库路径如/opt/local/lib

四、测试：

tetst.f

! gfortran -o test test.f -L/opt/local/lib -llapack -lrefblas

program main

implicit none

real :: a(3,3),b(3)

integer :: v(3),iflag

!external sgesv

a=reshape([2.0,0.0,0.0,0.0,3.0,0.0,0.0,0.0,4.0],[3,3])

b=[998.0,999.0,1000.0]

print *,'a=',a

print *,'b=',b

call sgesv(3,1,a,3,v,b,3,iflag)

print *,'solve=',b

end program main

gfortran -o test test.f -L/opt/local/lib -llapack -lrefblas

$ ./test

a=   2.00000000       0.00000000       0.00000000       0.00000000       3.00000000       0.00000000       0.00000000       0.00000000       4.00000000

b=   998.000000       999.000000       1000.00000

solve=   499.000000       333.000000       250.000000

test2.s

! gfortran -o test2 test2.f -L/opt/local/lib -llapack -lrefblas

program Console1

external dgesv

integer n, lda, nrhs, ldb

parameter (n=2,lda=2,nrhs=1,ldb=2)

double precision A(lda,n)

double precision b(ldb,nrhs)

character byebye

integer ipiv(n), info, i, j

A(1,1)=1

A(1,2)=2

A(2,1)=3

A(2,2)=4

B(1,1)=5

B(2,1)=6

write(*,*) 'Hello World'

call dgesv(n,nrhs,A,lda,ipiv,b,ldb,info)

write(*,*) 'INFO =', info

write(*,*) ((A(i,j),i=1,lda),j=1,n)

write(*,*) ((B(i,j),i=1,ldb),j=1,nrhs)

write(*,*) "END OF PROGRAM..."

end program Console1

$ gfortran -o test2 test2.f -L/opt/local/lib -llapack -lrefblas

$ ./test2

Hello World

INFO =           0

3.0000000000000000       0.33333333333333331        4.0000000000000000       0.66666666666666674

-3.9999999999999987        4.4999999999999991

或直接安装

sudo apt install libopenblas-base libopenblas-dev liblapacke liblapacke-dev liblapack3 liblapack-dev libtmglib3 libtmglib-dev

编译： gfortran -o test2 test2.f -llapack

五、CLAPACK的安装与简单使用

CLAPACK是没有源的，只能从其主页上下载源码自己编译http://www.netlib.org/clapack/clapack.tgz，解压。在其目录下的INSTALL目录中有make.inc的范本文件，比如在linux下就是make.inc.LINUX，将这个文件拷到CLAPACK目录下并改名为make.inc。

首先编译F2CLIBS，用于将fortran转换为c语言，在CLAPACK目录下，make f2clib，在

CLAPACK/F2CLIBS下就会生成libf2c.a。

编译tmglib，在CLAPACK目录下，make tmglib，会生成CLAPACK/tmglib_LINUX.a。

编译blas，在CLAPACK目录下，make blaslib，会生成CLAPACK/blas_LINUX.a。

最后是编译CLAPACK，make，会生成CLAPACK/lapack_LINUX.a

CLAPACK下的clapack.h就是所需要的头文件，除此之外还需要的一个头文件是F2CLIBS/f2c.h。

为方便使用，将必要的include文件和lib文件集中到一个目录下(如/mybk/CLAPACK)：

~/prg/Fortran95/lapack$ tree /mybk/CLAPACK

/mybk/CLAPACK

├── include

│   ├── blaswrap.h

│   ├── clapack.h

│   └── F2CLIBS

│   └── f2c.h

└── lib

├── blas_LINUX.a

├── lapack_LINUX.a

├── libcblaswr.a

├── libf2c.a

└── tmglib_LINUX.a

现在就通过使用CLAPACK中的一个函数sgesv_解线性方程组来学习一下使用的方法。

包括此函数在内的所有函数可以在CLAPACK/SRC下找到源代码，并在代码中有函数参数的说明信息。sgesv_的代码文件就是sgesv.c。

int sgesv_(integer *n, integer *nrhs, real *a, integer *lda, integer *ipiv, real *b, integer *ldb, integer *info)

sgesv_的功能是使用LU分解法解线性方程组AX=B，其中A是一个n*n的方阵。

integer *n, 方程的个数，也就是A的行数和列数

integer *nrhs, B的列数

real *a, 存储矩阵A数据的一维数组，在fortran中，数组是列主序存储，在此a中的二维数据也必须是列主序

integer *lda, 等于n

integer *ipiv, 一个输出数据数组，数组大小是n，具体什么功能不太明白，但是似乎不影响最后结果，谁明白请告诉我

real *b,存储矩阵B数据的一维数组，在fortran中，数组是列主序存储，在此b中的二维数据也必须是列主序

integer *ldb, 等于n

integer *info,输出参数，如果返回此参数为0，表示函数正常退出，否则表示出错。

代码：

// g++ -o clapackdemo clapackdemo.cpp -I/mybk/CLAPACK/include /mybk/CLAPACK/lib/lapack_LINUX.a /mybk/CLAPACK/lib/blas_LINUX.a /mybk/CLAPACK/lib/tmglib_LINUX.a /mybk/CLAPACK/lib/libf2c.a

#include

#include

using namespace std;

#include

//因为程序是C++，而CLAPACK是C语言写的，所以在此处用extern关键字

extern"C"

{

#include

}

int main(void){

integer M=3;

integer N=1;

real a[9]={4,3,11,2,-1,0,-1,2,3};

real b[3]={2,10,8};

integer lda;

integer ldb;

integer INFO;

lda=M;

ldb=M;

integer ipiv[M];

sgesv_(&M, &N, a, &lda,ipiv, b, &ldb, &INFO);

if(INFO==0)

for(int i=0;i

cout<

else

cout<

return EXIT_SUCCESS;

}

编译运行：

$ g++ -o clapackdemo clapackdemo.cpp -I/mybk/CLAPACK/include /mybk/CLAPACK/lib/lapack_LINUX.a /mybk/CLAPACK/lib/blas_LINUX.a /mybk/CLAPACK/lib/tmglib_LINUX.a /mybk/CLAPACK/lib/libf2c.a

$ ./clapackdemo

-14

56

54

c语言openblas：

//g++ -o openblastest1 openblastest1.cpp `pkg-config --cflags --libs blas-netlib`

#include

#include

#include

#include

// 随机生成20以内的给定尺寸数组

static void RandomFill(std::vector& numbers,size_t size);

// 打印数组元素的函数

static void Print(const std::vector& numbers);

// vector是一维的，输出是个矩阵，那输出的时候就要指定有几行几列

static void Print(const std::vector& numbers, int rows, int cols);

// 寻找数组中最大的那个元素的索引和值

static void TestLevel1();

// 测试Level2里面最常用的函数：向量和矩阵的乘积

static void TestLevel2();

static void TestLevel3();

int main(int argc, const char * argv[]) {

TestLevel1();

std::cout <

TestLevel2();

std::cout <

TestLevel3();

return 0;

}

void RandomFill(std::vector& numbers, size_t size) {

// 预分配size的缓冲区，这样性能相对更好一点

numbers.resize(size);

for (size_t i = 0; i != size; ++ i) {

numbers[i] = static_cast(rand() % 20);

}

}

static void TestLevel1() {

const int VECTOR_SIZE = 4;

std::vector fv1;

RandomFill(fv1, VECTOR_SIZE);

Print(fv1);

/**

从数组里面找出最大的那个数的索引：cblas_isamax()属于Level1的函数

VECTOR_SIZE 数组长度

fv1.data 数组缓冲区指针的首地址，怎么获得fv1内部缓冲区呢，用fv1.data

third params: 跳跃数量，两个元素之间间隔几个元素，也就是每处理一个元素之后，+1得到下一个元素，如果想跳过一个元素 +2

*/

size_t maxIndex = cblas_isamax(VECTOR_SIZE, fv1.data(), 1);

std::cout << maxIndex << std::endl;

std::cout << fv1[maxIndex] << std::endl;

}

/*

Prints:

7 9 13 18

3

18

*/

/*

测试Level2里面最常用的函数：向量和矩阵的乘积

*/

static void TestLevel2()

{

// 假设我们有一个三行二列的矩阵

const int M = 3;

const int N = 2;

/*

A(M*N),x(N*1), y(M*1)

我们定义三个矩阵a、x、y

a：M*N的矩阵：维度为M*N矩阵

x：N*1的矩阵，实际上是一个向量，维度为N矩阵

y：M*1的矩阵，实际上是长度为M的向量，维度为M的矩阵

*/

std::vector a;

std::vector x;

std::vector y;

RandomFill(a, M * N);

RandomFill(x, N);

RandomFill(y, M);

std::cout << "A" << std::endl;

Print(a, M, N);

std::cout << "x" << std::endl;

Print(x);

std::cout << "y" << std::endl;

Print(y);

/*

我们的目标是想计算这么一个公式：

y := alpha * A * x + beta * y

A:是一个矩阵，x是一个向量，所以我希望说去计算一个矩阵和向量的乘积。alpha是一个乘积的缩放，

beta是对y的缩放，

相当于把y里面的数字乘以beta，再加上A矩阵和向量的乘积。

那这边有一个特例，假如我y里面都是0，或这beta是0的情况下，我就可以把公式看成：

// y := alpha * A * x

这个函数名称为：cblas_sgemv()

// s:single 单精度浮点数

// ge: 是一个乘法

// m: matrix

// v: vector

*/

/**

参数解释：

param CblasRowMajor 行主序还是列主序，默认行主序，何为主序：即数组存储元素的方式--按行存储还是按列存储，行主序：00，01，列主序00，10

param CblasNoTrans 矩阵是否需要转置，不需要转置，如果需要转置的话，运算的时候它会自动做转置

param M 矩阵的行数

param N 矩阵的列数

param 1.0f alpha ，我们设为1

param a.data a矩阵的缓冲区首地址

param lda a矩阵的列数

param x.data x矩阵的缓冲区首地址

param 1 x里面每次跳跃累加的个数，默认为1

param 2.0f beta对y的缩放值

param y.data y矩阵的缓冲区首地址

param 1 y里面每次跳跃累加的个数，默认为1

*/

int lda = N;

cblas_sgemv(CblasRowMajor, CblasNoTrans, M, N, 1.0f, a.data(), lda, x.data(), 1, 2.0f, y.data(), 1);

std::cout << "result y" << std::endl;

Print(y);

}

/*

Prints:

A

7 9

13 18

10 12

x

4 18

y

3 9 0

result y

196 394 256

*/

void Print(const std::vector& numbers) {

for (float number : numbers) {

std::cout << number << ' ';

}

std::cout << std::endl;

}

void Print (const std:: vector& numbers ,int rows, int cols) {

for (int row =0; row != rows; ++ row) {

for (int col = 0; col != cols; ++ col) {

// 取出每一列的数字

std::cout << numbers[row * cols + col] << ' ';

}

// 没输出一行之后换一行

std::cout << std::endl;

}

}

/**

计算两个矩阵的乘积

*/

static void TestLevel3() {

// 我们希望计算两个矩阵的乘积，我们就需要定义三个参数M、N、K。

const int M = 3;

const int N = 2;

const int K = 4;

std::vector a;

std::vector b;

std::vector c;

RandomFill(a, M * K);

RandomFill(b, K * N);

RandomFill(c, M * N);

// 输出A、B、C 三个矩阵

std::cout << "A" << std::endl;

Print(a, M, K);

std::cout << "B" << std::endl;

Print(b, K, N);

std::cout << "C" << std::endl;

Print(c, M, N);

/*

我们的目标是计算这么一个公式：

// C := alpha * A * B + beta * C

// 如果只想做两个矩阵的乘法，beta设成0就好了，变为如下式子：

// C := alpha * A * B

*/

/*

函数释义：

sgemm:矩阵间的单精度乘法

s:single 单精度

ge：general 一般的，普通的

m：matix 矩阵

m：multiplication 乘法

*/

/**

参数释义：

param CblasRowMajor 行主序还是列主序，默认行主序

param CblasNoTrans A需不需要转置，不需要

param CblasNoTrans B需不需要转置，不需要

param M 系数M

param N 系数N

param K 系数K

param 1.0f alpha 设为1

param a.data a的缓冲区首地址

param lda a的列数

param b.data b的缓冲区首地址

param ldb b的列数

param 1.0f beta 设为1

param c.data c的缓冲区首地址

param ldc c的列数

*/

int lda = K;

int ldb = N;

int ldc = N;

cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.0f, a.data(), lda, b.data(), ldb, 1.0f, c.data(), ldc);

std::cout << "Result C" << std::endl;

// 三行四列的矩阵 * 四行二列的矩阵 + 三行二列的矩阵，结果为一个三行二列的矩阵

Print(c, M, N);

}

/*

Prints:

A

7 9 13 18

10 12 4 18

3 9 0 5

B

12 2

7 3

7 9

0 12

C

3 9

9 17

0 13

Result C

241 383

241 325

99 106

*/

编译运行：

$ g++ -o openblastest1 openblastest1.cpp `pkg-config --cflags --libs blas-netlib`

$ ./openblastest1

3 6 17 15

2

17

-------------------------------------

A

13 15

6 12

9 1

x

2 7

y

10 19 3

result y

151 134 31

-------------------------------------

A

6 0 6 12

16 11 8 7

9 2 10 2

B

3 7

15 9

2 2

18 9

C

7 13

16 11

2 9

Result C

253 175

371 301

115 128