Fortran Coder

查看: 903|回复: 7
打印 上一主题 下一主题

[混编] 基于mex文件的Matlab和Fortran混合编程问题

[复制链接]

5

帖子

1

主题

0

精华

入门

F 币
31 元
贡献
16 点
跳转到指定楼层
楼主
发表于 2023-11-3 12:14:04 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
最近正在学习基于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.

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

798

帖子

2

主题

0

精华

大宗师

F 币
3793 元
贡献
2268 点
沙发
发表于 2023-11-3 20:13:47 | 只看该作者
你把fortran里面所有的计算语句都删除,看看有没有问题。这样可以确定到底是计算问题还是接口问题。

5

帖子

1

主题

0

精华

入门

F 币
31 元
贡献
16 点
板凳
 楼主| 发表于 2023-11-4 12:51:51 | 只看该作者
li913 发表于 2023-11-3 20:13
你把fortran里面所有的计算语句都删除,看看有没有问题。这样可以确定到底是计算问题还是接口问题。 ...

好的,谢谢您回复,我试试

5

帖子

1

主题

0

精华

入门

F 币
31 元
贡献
16 点
地板
 楼主| 发表于 2023-11-4 18:53:34 | 只看该作者
本帖最后由 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

798

帖子

2

主题

0

精华

大宗师

F 币
3793 元
贡献
2268 点
5#
发表于 2023-11-5 19:25:19 | 只看该作者
matlab 默认浮点数是双精度,整数是8字节。

5

帖子

1

主题

0

精华

入门

F 币
31 元
贡献
16 点
6#
 楼主| 发表于 2023-11-12 11:22:04 | 只看该作者
li913 发表于 2023-11-5 19:25
matlab 默认浮点数是双精度,整数是8字节。

抱歉,我没太明白您的意思,是声明的问题吗,我又声明了一下变量的精度,浮点数为双精度,整数为8字节,但是仍然出现了同样的问题

5

帖子

1

主题

0

精华

入门

F 币
31 元
贡献
16 点
7#
 楼主| 发表于 2023-11-29 22:02:04 | 只看该作者
已经解决了。谢谢大家

6

帖子

2

主题

0

精华

入门

F 币
41 元
贡献
12 点
8#
发表于 2024-3-8 16:16:41 | 只看该作者
AsIlll 发表于 2023-11-29 22:02
已经解决了。谢谢大家

楼主怎么解决的?就是改了变量声明的精度吗?
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-29 05:47

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表