Fortran Coder

标题: 基于mex文件的Matlab和Fortran混合编程问题 [打印本页]
作者: AsIlll    时间: 2023-11-3 12:14
标题: 基于mex文件的Matlab和Fortran混合编程问题
最近正在学习基于mex文件的Matlab和Fortran的混合编程，编写了下面这样一个f90文件，在matlab中编译成mex_function也没有报错，但是一调用生成的function，Matlab就直接崩溃，我是初学者不知道是什么问题，求助大家看看代码有没有什么问题，感激不尽。
[Fortran] 纯文本查看 复制代码
#include "fintrf.h"





subroutine mexFunction(OutSum,OutVar,InSum,InVar)         !函数接口名称必须为mexFunction,



    !OutSum:输出参数个数

    !OutVar:输出参数数组指针

    !InSum:输入参数个数

    !InVar:输入参数数组指针



    !参数顺序不能随意更改



    Integer InSum,OutSum

    

    mwPointer InVar(*),OutVar(*)                           !mwPointer专门用于表示指针变量,这个不能随意用Integer代替

    

    mwPointer mxGetPr, mxGetDimensions, mxCreateNumericArray  !这个对返回指针函数的再次申明,

    integer, parameter :: dp = 8

    Integer , parameter :: myINT  = SELECTED_INT_KIND(8)



    Real(dp),Allocatable :: P(:,:,:,:), F1(:,:,:,:), F2(:,:,:,:), Gamma(:,:,:,:), rc(:,:,:,:), uind(:,:,:,:) ! input and output uind

    Integer,Allocatable :: Dim_P(:), Dim_F(:), Dim_G(:), Dim_rc(:) ! input array dimension

    Real(dp) :: n,co   ! input n and cut_off distance

    Real(dp) :: uindx,uindy,uindz,temuindx,temuindy,temuindz



    Integer I1,J1,K1,I2,J2,K2,NBP,NTP,NSP,NBF,NTF,NSF,PD,FD,GD,RD



    If(InSum/=7)Then

    

        call mexErrMsgIdAndTxt('MATLAB:InputTooBig','输入参数个数必须为7个')

    

        Return

    

    EndIf



    PD = mxGetNumberOfDimensions(P)

    FD = mxGetNumberOfDimensions(F1)

    GD = mxGetNumberOfDimensions(Gamma)

    RD = mxGetNumberOfDimensions(rc)





    Allocate(Dim_P(PD), Dim_F(FD), Dim_G(GD), Dim_rc(RD))



    Dim_P = mxGetDimensions(InVar(1)) !获取第1个输入参数P的维度

    Dim_F = mxGetDimensions(InVar(2)) !获取第2 (3)个输入参数F1 (F2)的维度

    Dim_G = mxGetDimensions(InVar(4)) !获取第4个输入参数Gamma的维度

    Dim_rc = mxGetDimensions(InVar(5)) !获取第5个输入参数rc的维度  

    

    

    Allocate(P(Dim_P(1),Dim_P(2),Dim_P(3),Dim_P(4)))

    Allocate(F1(Dim_F(1),Dim_F(2),Dim_F(3),Dim_F(4)))

    Allocate(F2(Dim_F(1),Dim_F(2),Dim_F(3),Dim_F(4)))

    Allocate(Gamma(Dim_G(1),Dim_G(2),Dim_G(3),Dim_G(4)))

    Allocate(rc(Dim_rc(1),Dim_rc(2),Dim_rc(3),Dim_rc(4)))



    Allocate(Uind(Dim_P(1),Dim_P(2),Dim_P(3),Dim_P(4)))

    



    Call mxCopyPtrToReal8(mxGetPr(InVar(1)),P,Dim_P(1)*Dim_P(2)*Dim_P(3)*Dim_P(4)) !将第1个参数数组赋值给P变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(2)),F1,Dim_F(1)*Dim_F(2)*Dim_F(3)*Dim_F(4)) !将第2个参数数组赋值给F1变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(3)),F2,Dim_F(1)*Dim_F(2)*Dim_F(3)*Dim_F(4)) !将第3个参数数组赋值给F2变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(4)),Gamma,Dim_G(1)*Dim_G(2)*Dim_G(3)*Dim_G(4)) !将第4个参数数组赋值给Gamma变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(5)),rc,Dim_rc(1)*Dim_rc(2)*Dim_rc(3)*Dim_rc(4)) !将第5个参数数组赋值给rc变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(6)),n,1)!将第6个整数变量赋值给n

    Call mxCopyPtrToReal8(mxGetPr(InVar(7)),co,1)!将第7个整数变量赋值给co cut_off_distance

    

    NBP = Dim_P(1)

    NTP = Dim_P(3)

    NSP = Dim_P(4)

    NBF = Dim_F(1)

    NTF = Dim_F(3)

    NSF = Dim_F(4)



    DO K1 = 1,NBP

        DO J1 = 1,NTP

            DO I1 = 1,NSP



                temuindx = 0.0

                temuindy = 0.0

                temuindz = 0.0



                DO K2 = 1,NBF

                    DO J2 = 1,NTF

                        DO I2 = 1,NSF



                            call Biot_p2p(P(I1,1,J1,K1),P(I1,2,J1,K1),P(I1,3,J1,K1),F1(I2,1,J2,K2),F1(I2,2,J2,K2),F1(I2,3,J2,K2),F2(I2,1,J2,K2),F2(I2,2,J2,K2),F2(I2,3,J2,K2),&

                                    Gamma(I1,1,J1,K1),rc(I1,1,J1,K1),co,n,temuindx,temuindy,temuindz)



                            

                            temuindx = temuindx + temuindx

                            temuindy = temuindy + temuindy

                            temuindy = temuindy + temuindy







                        ENDDO

                    ENDDO

                ENDDO





                Uind(I1,1,J1,K1) = temuindx

                Uind(I1,2,J1,K1) = temuindy

                Uind(I1,3,J1,K1) = temuindz









            ENDDO

        ENDDO

    ENDDO



    OutVar(1)=mxCreateNumericArray(PD,Dim_P,REAL*8,0)!给返回参数分配内存



    Call mxCopyReal8ToPtr(Uind,mxGetPr(OutVar(1)),Dim_P(1)*Dim_P(2)*Dim_P(3)*Dim_P(4))!将返回参数赋值给分配的内存

    

    DeAllocate(P,F1,F2,Gamma,rc)!释放临时分配的内存



    Return



End SubRoutine



Subroutine Biot_p2p(px,py,pz,f1x,f1y,f1z,f2x,f2y,f2z,gamma,rc,co,n,temuindx,temuindy,temuindz)



    implicit none



    integer, parameter :: dp = 8

    Real(dp),Intent(In):: px,py,pz,f1x,f1y,f1z,f2x,f2y,f2z,gamma

    Real(dp),Intent(Out):: temuindx,temuindy,temuindz

    Real(dp) :: ldx, ldy, ldz, pxx1, pxx2, pyy1, pyy2, pzz1, pzz2, r1, r2, ldr, len, r1dr2, r1tr2, cnu, den, ubar,cmn,rc

    Real ::  co,cutoff,cored,n

    Real, parameter :: pi = acos(-1.0)

    Real, parameter :: infinity = 999999.0



    cutoff = co

    cored = n



    

    temuindx = 0.0

    temuindy = 0.0

    temuindz = 0.0



    ldx = f2x-f1x

    ldy = f2y-f1y

    ldz = f2z-f1z



    pxx1 = px-f1x

    pxx2 = px-f2x

    pyy1 = py-f1y

    pyy2 = py-f2y

    pzz1 = pz-f1z

    pzz2 = pz-f2z



    cmn = cored



    r1 = sqrt(pxx1*pxx1+pyy1*pyy1+pzz1*pzz1)

    r2 = sqrt(pxx2*pxx2+pyy2*pyy2+pzz2*pzz2)





    if((r1<cutoff).AND.(r2<cutoff)) then

    

        ldr = (ldx*pxx1 + ldy*pxx2 + ldz*pzz1)

        ldr = ldr*ldr

        len = ldx*ldx + ldy*ldy + ldz*ldz !!L^2

        r1dr2 = pxx1*pxx2+pyy1*pyy2+pzz1*pzz2

        r1tr2 = r1*r2



        cnu = (r1*r1)-(ldr/len)

        cnu = cnu*((rc**cmn+cnu**cmn)**(-1.0/cmn))



        den = 4.0*pi*(r1tr2*(r1tr2 + r1dr2))



        ubar = cnu*gamma*(r1+r2)

        ubar = ubar/den



        if (isnan(ubar)) then

           temuindx = 0.0

           temuindy = 0.0

           temuindz = 0.0 

        else if (ubar.GE.infinity) then

           temuindx = 0.0

           temuindy = 0.0

           temuindz = 0.0         

        else

            temuindx = ubar*(pyy1*pzz2-pzz1*pyy2)

            temuindy = ubar*(pzz1*pxx2-pxx1*pzz2)

            temuindz = ubar*(pxx1*pyy2-pyy1*pxx2)

        endif



    endif





End SubRoutine Biot_p2p


崩溃信息：This error was detected while a MEX-file was running. If the MEX-file
is not an official MathWorks function, please examine its source code
for errors. Please consult the External Interfaces Guide for information
on debugging MEX-files.

作者: li913    时间: 2023-11-3 20:13
你把fortran里面所有的计算语句都删除，看看有没有问题。这样可以确定到底是计算问题还是接口问题。
作者: AsIlll    时间: 2023-11-4 12:51
li913 发表于 2023-11-3 20:13
你把fortran里面所有的计算语句都删除，看看有没有问题。这样可以确定到底是计算问题还是接口问题。 ...

好的，谢谢您回复，我试试
作者: AsIlll    时间: 2023-11-4 18:53
本帖最后由 AsIlll 于 2023-11-4 18:54 编辑

我将计算语句全部删除之后，仍然出现了程序崩溃的问题，应该是接口的问题，但是我没看出来是什么问题，请问有大佬帮我看看吗
[Fortran] 纯文本查看 复制代码
#include "fintrf.h"





subroutine mexFunction(OutSum,OutVar,InSum,InVar)         !函数接口名称必须为mexFunction,



    !OutSum:输出参数个数

    !OutVar:输出参数数组指针

    !InSum:输入参数个数

    !InVar:输入参数数组指针



    !参数顺序不能随意更改



    Integer InSum,OutSum

    

    mwPointer InVar(*),OutVar(*)                           !mwPointer专门用于表示指针变量,这个不能随意用Integer代替

    

    mwPointer mxGetPr, mxGetDimensions, mxCreateNumericArray  !这个对返回指针函数的再次申明,

    integer, parameter :: dp = 8

    Integer , parameter :: myINT  = SELECTED_INT_KIND(8)



    Real(dp),Allocatable :: P(:,:,:,:), F1(:,:,:,:), F2(:,:,:,:), Gamma(:,:,:,:), rc(:,:,:,:), uind(:,:,:,:) ! input and output uind

    Integer,Allocatable :: Dim_P(:), Dim_F(:), Dim_G(:), Dim_rc(:) ! input array dimension

    Real(dp) :: n,co   ! input n and cut_off distance

    Real(dp) :: uindx,uindy,uindz,temuindx,temuindy,temuindz



    Integer I1,J1,K1,I2,J2,K2,NBP,NTP,NSP,NBF,NTF,NSF,PD,FD,GD,RD



    If(InSum/=7)Then

    

        call mexErrMsgIdAndTxt('MATLAB:InputTooBig','输入参数个数必须为7个')

    

        Return

    

    EndIf



    PD = mxGetNumberOfDimensions(P)

    FD = mxGetNumberOfDimensions(F1)

    GD = mxGetNumberOfDimensions(Gamma)

    RD = mxGetNumberOfDimensions(rc)





    Allocate(Dim_P(PD), Dim_F(FD), Dim_G(GD), Dim_rc(RD))



    Dim_P = mxGetDimensions(InVar(1)) !获取第1个输入参数P的维度

    Dim_F = mxGetDimensions(InVar(2)) !获取第2 (3)个输入参数F1 (F2)的维度

    Dim_G = mxGetDimensions(InVar(4)) !获取第4个输入参数Gamma的维度

    Dim_rc = mxGetDimensions(InVar(5)) !获取第5个输入参数rc的维度  

    

    

    Allocate(P(Dim_P(1),Dim_P(2),Dim_P(3),Dim_P(4)))

    Allocate(F1(Dim_F(1),Dim_F(2),Dim_F(3),Dim_F(4)))

    Allocate(F2(Dim_F(1),Dim_F(2),Dim_F(3),Dim_F(4)))

    Allocate(Gamma(Dim_G(1),Dim_G(2),Dim_G(3),Dim_G(4)))

    Allocate(rc(Dim_rc(1),Dim_rc(2),Dim_rc(3),Dim_rc(4)))



    Allocate(Uind(Dim_P(1),Dim_P(2),Dim_P(3),Dim_P(4)))

    



    Call mxCopyPtrToReal8(mxGetPr(InVar(1)),P,Dim_P(1)*Dim_P(2)*Dim_P(3)*Dim_P(4)) !将第1个参数数组赋值给P变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(2)),F1,Dim_F(1)*Dim_F(2)*Dim_F(3)*Dim_F(4)) !将第2个参数数组赋值给F1变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(3)),F2,Dim_F(1)*Dim_F(2)*Dim_F(3)*Dim_F(4)) !将第3个参数数组赋值给F2变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(4)),Gamma,Dim_G(1)*Dim_G(2)*Dim_G(3)*Dim_G(4)) !将第4个参数数组赋值给Gamma变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(5)),rc,Dim_rc(1)*Dim_rc(2)*Dim_rc(3)*Dim_rc(4)) !将第5个参数数组赋值给rc变量

    Call mxCopyPtrToReal8(mxGetPr(InVar(6)),n,1)!将第6个整数变量赋值给n

    Call mxCopyPtrToReal8(mxGetPr(InVar(7)),co,1)!将第7个整数变量赋值给co cut_off_distance

    

    NBP = Dim_P(1)

    NTP = Dim_P(3)

    NSP = Dim_P(4)

    NBF = Dim_F(1)

    NTF = Dim_F(3)

    NSF = Dim_F(4)





    OutVar(1)=mxCreateNumericArray(PD,Dim_P,REAL*8,0)!给返回参数分配内存



    Call mxCopyReal8ToPtr(Uind,mxGetPr(OutVar(1)),Dim_P(1)*Dim_P(2)*Dim_P(3)*Dim_P(4))!将返回参数赋值给分配的内存

    

    DeAllocate(P,F1,F2,Gamma,rc)!释放临时分配的内存



    Return



End SubRoutine
作者: li913    时间: 2023-11-5 19:25
matlab 默认浮点数是双精度，整数是8字节。
作者: AsIlll    时间: 2023-11-12 11:22
li913 发表于 2023-11-5 19:25
matlab 默认浮点数是双精度，整数是8字节。

抱歉，我没太明白您的意思，是声明的问题吗，我又声明了一下变量的精度，浮点数为双精度，整数为8字节，但是仍然出现了同样的问题
作者: AsIlll    时间: 2023-11-29 22:02
已经解决了。谢谢大家
作者: Alohafancy    时间: 2024-3-8 16:16
AsIlll 发表于 2023-11-29 22:02
已经解决了。谢谢大家

楼主怎么解决的？就是改了变量声明的精度吗？



欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2