Fortran调用mkl计算逆矩阵_MATLAB-仿真秀干货文章

Fortran调用mkl计算逆矩阵

零

概述

fortran，编程语言的老爷车，但在科学计算领域依然具有不可撼动的地位。

使用的fortran的时候，常常有矩阵操作的需求，但是fortan很少有内置的矩阵操作函数，这时候需要调用第三方计算库进行矩阵计算。

在vs2013平台，以求解逆矩阵为例讲解如何在vs里面调用mkl计算库。

壹

链接mkl计算库

第一步：在上部菜单栏中选择项目->属性

在“属性页”窗口，依次选择Fortran ->Libraries-> Use Intel Kernel Library->parallel（/Qmkl:parallel）点击确定。

第二步：在“解决方案管理器”中右键解决方案，选择属性

在“属性页”中依次选择Linker->Input->Additional Denendencies，输入：mkl_lapack95.lib

下面的代码，测试Lapack95是否链接成功

      program main            use lapack95            write(*,*)"Lpack95链接成功"            pause      end program main

当正常运行，输出一下内容，说明链接成功

贰

利用方程组求解逆矩阵

求解逆矩阵可以直接根据线性代数中的求解方法，但是过于复杂，在大规模计算中基本是不适用的。fortan采用LU分解的方式求解，这里以方程组AX=B为例，当B为单位矩阵时候，方程组的解X就是A矩阵的逆矩阵。

调用了mkl计算库中的dgesv子程序，该子程序用于求解线性方程组，该子程序的函数头为

lapack_int LAPACKE_sgesv( n,      // (input) 线性方程的个数,n>=0                         nrhs,      // (input) 矩阵B的列数，即线性方程组右端的项个数,nrhs>=0                         a,       // (input/output)系数矩阵A，维度为nxn                         lda,       // (input) A矩阵的第一维                         ipiv,      // (output) 置换矩阵，ipiv[i]表示矩阵A的第i行与第ipiv[i]行进行了交换                         b,      // (input/output)B矩阵                         ldb       // (input) B矩阵的第一维                        info     //储存结算结果有关参数，info为0时计算正确                        )

详细的代码可以参考：

LAPACK: dgesv (netlib.org)

首先采用matlab计算逆矩阵，然后将矩阵的维度、初始矩阵和matlab计算的逆矩阵都输出到外部txt文件中。

其中，维度dim用于fortran定义矩阵，初始矩阵用于fortran读取计算逆矩阵，matlab求解的逆矩阵用于和fortran的计算结果进行对比。

首先是matlab程序

clc;    clear;path='C:\Users\nnn\Desktop\inv_martix\';%%N=100;   %   矩阵大小
m=rand(N,N);      %   随机矩阵 inv_m=inv(m);   %   matlab求解逆矩阵
%%  输出矩阵维度fid=fopen(path+"dim.txt",'w');fprintf(fid,"%d",N);fclose(fid);
% 输出随机矩阵供fortran读求解逆矩阵OutputMartix(path+"m.txt",m);
% 输出matlab求解的逆矩阵OutputMartix(path+"inv_m.txt",inv_m);
%%  定义输出矩阵的函数function OutputMartix(path,m)% 打开文件fid1=fopen(path,'w');
for i=1:size(m,1)    for ii=1:size(m,1)        if ii==size(m,1)            fprintf(fid1,"%f\n",m(i,ii));        else            fprintf(fid1,"%f\t",m(i,ii));        end            endend%关闭文件fclose(fid1);end%%

上面的程序向我们请求一个矩阵大小参数，然后根据维度生成一个随机矩阵；程序运行结束会生成三个文件，分别为

m.txt      储存原始矩阵dim.txt    储存矩阵维度inv_m.txt  储存matlab计算的逆矩阵

当求解的矩阵为5X5时候matlab生成的初始随机矩阵为

0.238472  0.887758  0.062132  0.503891  0.7421880.158590  0.557903  0.603668  0.408332  0.3658670.038231  0.363472  0.110736  0.632575  0.3978250.611598  0.373223  0.183794  0.758511  0.5724170.426455  0.536710  0.507714  0.643364  0.089276

matlab求解的逆矩阵为

0.282723  -0.567128  -2.133005  1.384150  0.6038691.590926  -0.857213  -0.536741  -1.365810  1.436023-1.019661  2.210717  -0.330703  0.223589  -0.542941-0.697804  -0.878189  2.135785  -0.129090  0.710454-0.087321  1.618724  -0.094985  1.257878  -2.348602

下面是完整的fortran程序

      program main          use lapack95            integer info,dim      C     声明动态数组            double precision, dimension (:,:), allocatable::a,b,matlab_inv_a,     1 error      integer, dimension (:), allocatable ::ipivCC     读取方程左端量维度C      open(100, file='C:\Users\nnn\Desktop\inv_martix\dim.txt',      1 status='old')        read(100,*)dim      allocate(a(dim,dim))      allocate(b(dim,dim))      allocate(matlab_inv_a(dim,dim))      allocate(ipiv(dim))      allocate(error(dim,dim))      close(100)      C   C     读取方程左端量AC            open(100, file='C:\Users\nnn\Desktop\inv_martix\m.txt',      1 status='old')        do i=1,dim          read(100,*)(a(i,ii),ii=1,dim)      enddo      close(100)      C      C     读取方程matlab的解答Inv_A C            open(100, file='C:\Users\nnn\Desktop\inv_martix\inv_m.txt',      1 status='old')        do i=1,dim          read(100,*)(matlab_inv_a(i,ii),ii=1,dim)      enddo      close(100)                        C  AX=B的B为单位矩阵            do i=1,dim        do ii=1,dim          if(i.eq.ii) then            b(i,ii)=1.d0          else            b(i,ii)=0.d0            endif        enddo      enddo      C     mkl lapack 求解AX=B            call dgesv(dim,dim,a,dim,ipiv,b,dim,info)
C     检查计算结果            if(info.eq.0)then          write(*,*)"计算结束！"      elseif(info.lt.0)then          write(*,*)"第",info,"个参数出现错误"      endif      C     计算fortan与matlab的误差      error=b-matlab_inv_a            C     输出误差            do i=1,dim                      write(*,*)"================================================"          write(*,*)(error(i,ii),ii=1,dim)      enddo            pause      end program main

fortan程序的矩阵采用了动态数组，涉及到allocate的使用，该参数允许我们在运行时刻定义数组。

文件的读取涉及到open和close的使用，采用逐行读取并赋值的方式。

程序最后的数组error用于比较fortan与matlab的计算误差，当matlab中矩阵行列为5时候，fortran计算的逆矩阵为

  0.282724810474718      -0.567127932917041       -2.13300554826328   1.38414882116874       0.603869463986308   1.59092466221457      -0.857213823820332      -0.536744751112424  -1.36580656716275        1.43602351412321  -1.01965883092666        2.21071417933664      -0.330702302220375  0.223585458186896      -0.542937520517426 -0.697807460769760      -0.878187005828073        2.13578796695868 -0.129088088425678       0.710452275103190 -8.731710238368200E-002   1.61872328492087      -9.498347204982187E-002   1.25787466658899       -2.34860174914708

二者的误差error矩阵为

  1.810474718144661E-006  6.708295874346959E-008 -5.482632801090404E-007 -1.178831259940338E-006  4.639863080413420E-007 -1.337785428878746E-006 -8.238203322852300E-007 -3.751112424343894E-006  3.432837253347643E-006  5.141232108929472E-007  2.169073342717098E-006 -2.820663357105957E-006  6.977796250806634E-007 -3.541813103580260E-006  3.479482574397785E-006 -3.460769759988658E-006  1.994171926700261E-006  2.966958676076104E-006  1.911574322166487E-006 -1.724896809829346E-006  3.897616317999342E-006 -7.150791270227330E-007  1.527950178131787E-006 -3.333411010331133E-006  2.508529246547653E-007

比较发现，fortran的计算结果与matlab一致，事实上，二者内部采用的计算库都是相同的计算库，所有的误差都是数据截断产生的，二者的计算结果是可以互相替代的。

以上就是fortran在vs里面调用mkl计算库计算逆矩阵的过程，通过链接mkl计算库，可以大大增强fortran的计算性能。

来源：有限元先生

C++函数基础概念

零概述为了程序的模块化管理，任何计算机语言都会设计函数，将多次执行的程序封装为函数，以供任何时候调用。下面介绍C++中的函数使用。壹函数简介数据类型就是前文讲解的内容，经常使用的数据类型有：int、float等等，函数名称的命名与变量名称的命名要保持相同的法则，函数名后面开括号里的内容是调用函数的参数。依据输入和输出的类型，C++函数可以大致分为以下几类有输入，有输出有输入，无输出无输入，有输出无输入，无输出一个函数的完整使用过程，不但包括函数的具体定义，还要在主程序main函数之前进行声明，如#include<iostream>usingnamespacestd;\\函数声明输出类型函数名称（参数1，参数2...）intmain(){\\函数调用return0;}\\函额具体定义输出类型函数名称（参数1，参数2...）{程序语句return函数返回值；}下面依次介绍这几种函数贰有输入，有输出设计一个简单的数学函数，这个函数分为三个阶段，x在不同的阶段y有不同对应的函数值，如y=x-1x<0y=x^2x>=0&&x<10y=x^3x>+10对应的C++程序为#include<iostream>#include<string>usingnamespacestd;doublesubroutine(doublex);intmain(){doublex=0;cout<<"PleaseinputX:\t"<<endl;cin>>x;cout<<"TheYis:\t"<<subroutine(x)<<endl;system("pause");return0;}doublesubroutine(doublex){if(x<0){returnx-1;}elseif(x>=0&&x<10){returnx*x;}else{returnx*x*x;}}上面的程序首先向屏幕请求一个输入，然后将用户输入的数值传入函数中，函数判断x的范围并进行相应的计算，当输入的数值为56时，输出的内容为PleaseinputX:56TheYis:175616叁有输入，无输出这里设计一个函数，这个函向屏幕输出许多行字符，输出的行数是通过参数传进来，这个函数不向主程序返回任何数值，只是向屏幕输出以字符，如#include<iostream>#include<string>usingnamespacestd;voidsubroutine(intx);intmain(){intx=0;cout<<"输入行数:\t"<<endl;cin>>x;subroutine(x);system("pause");return0;}voidsubroutine(intx){for(inti=1;i<=x;i++){cout<<"第"<<i<<"行，共计"<<x<<"行。\n"<<endl;}}上面的程序接受一个整型参数，这个参数数值就是函向屏幕输入的字符串行数，当输入的数值为8时，输出的结果为输入行数:8第1行，共计8行。第2行，共计8行。第3行，共计8行。第4行，共计8行。第5行，共计8行。第6行，共计8行。第7行，共计8行。第8行，共计8行。肆无输入，有输出下面的函数没有接接受输入，但是向主程序返回了一个整型的数值，并将该数值赋值给一个整型变量，并且向屏幕输出了结果#include<iostream>#include<string>usingnamespacestd;\\definefunctionintsubroutine();intmain(){intoutput=0;\\callsubroutineoutput=subroutine();cout<<output<<endl;system("pause");return0;}\\definefunctionintsubroutine(){cout<<"这是一个无输入但是有数值输出的函数！"<<endl;cout<<"这个函数输出的数值内容为："<<endl;return23;}上面的函数输出结果为这是一个无输入但是有数值输出的函数！这个函数输出的数值内容为：23伍无输入，无输出无输入也无输出的函数定义为voidfunctionname(){\\codes}这种函数声明方式，不需要在函数中添加代码return，只需要设置输出类型为void即可。设置无输入也无输出函数也可以采用如下方式intfunctionname(){\\codesreturn0;}上面代码将函数输出类型设置为int，而且在函数程序中添加了return0，这意味着程序会返回一个0数值，这也是可行的。同样的道理，也可以stringfunctionname(){return"";}这就返回了一个空的字符串，也是可行的。下面给出了一个无输入也无输出的函数示例#include<iostream>#include<string>usingnamespacestd;voidsubroutine();intmain(){subroutine();system("pause");return0;}voidsubroutine(){cout<<"这是一个没有输入，也没有数值输出的函数。"<<endl;cout<<"这个函数仅仅向屏幕输出了这两句话。"<<endl;}运行上述程序，只会向屏幕输出两句话，而没有任何的数值输出，输出结果为这是一个没有输入，也没有数值输出的函数。这个函数仅仅向屏幕输出了这两句话。陆函数与数组的联合使用函数不但可以用来接收一些简单的数据，还可以接收数组和指针类型的数据，这里介绍一个接收数组类型的函数例子假设要计算一个班级里面所有人一学期的签到次数，就需要创建一个数组，数组的索引就是第几个人，数组的数值就是一个人签到的次数，这里假设一个班级里面有8个人，设计如下程序#include<iostream>usingnamespacestd;constintpeople=8;intsum(intarr[],intn);//第六行代码改成这样也是可以运行的,intsum(int[],int);intmain(){intdata;//班级里面每个人额签到次数intarr[8]={1,2,3,4,5,6,7,8};//调用函数data=sum(arr,people);cout<<"班级里所有人签到的总次数为：\t"<<data<<endl;system("pause");return0;}intsum(intarr[],intn){intm=0;for(inti=0;i<n;i++){m+=arr[i];}returnm;}注意第七行的注释，将intsum(intarrr[],intn)修改为intsum(int[],int)也是可以正确运行的，这是因为函数原型中的参数不需要准确的变量名称，只需要知道数据类型就行。关于函数、数组和指针的应用，会更加的有意思，下一部分帖子会重点讲解指针与数组的内容。来源：有限元先生