-
2016-12-04.png (25.22 KB, 下载次数: 293)
结果
subroutine wrapper_ODR(FCN,N,M,NP,NQ,BETA,Y,X,&
DELTA,WE,WD,IFIXB,IFIXX,JOB,NDIGIT,TAUFAC,&
SSTOL,PARTOL,MAXIT,IPRINT,LUNERR,LUNRPT,&
STPB,STPD,SCLB,SCLD,WORK,IWORK,INFO,LOWER,UPPER) bind(C, name='wrapper_ODR')
!DEC$ ATTRIBUTES DLLEXPORT :: wrapper_ODR
use iso_c_binding
use ODRPACK95
implicit none
!自訂意凾式介面宣告
interface
subroutine FCN(N,M,NP,NQ,LDN,LDM,LDNP,BETA,XPLUSD,IFIXB,IFIXX,LDIFX,&
IDEVAL,F,FJACB,FJACD,ISTOP) bind(C)
use, intrinsic :: iso_c_binding
implicit none
integer(c_int) :: IDEVAL,ISTOP,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ
real (c_double) :: BETA(1:NP),F(1:LDN,1:NQ),FJACB(1:LDN,1:LDNP,1:NQ), &
FJACD(1:LDN,1:LDM,1:NQ),XPLUSD(1:LDN,1:M)
integer(c_int) :: IFIXB(1:NP),IFIXX(1:LDIFX,1:M)
end subroutine
end interface
integer(c_int),value :: N,M,NP,NQ
real(c_double) :: BETA(1:NP),Y(1:N,1:NQ),X(1:N,1:M)
!!!!!Optional variable
integer(c_int), intent(in), optional :: IFIXB(:),IFIXX(:,:),JOB,NDIGIT,MAXIT&
,IPRINT,LUNERR,LUNRPT,IWORK(:),INFO
real(c_double), intent(in), optional :: DELTA(:,:),&
WE(:,:,:),WD(:,:,:),TAUFAC,SSTOL,PARTOL,&
STPB(:),STPD(:,:),SCLB(:),SCLD(:,:),&
WORK(:),LOWER(:),UPPER(:)
!!!!!Call ODR
call ODR(inter_func,N,M,NP,NQ,BETA,Y,X)
contains
!!!!!Local subroutine
subroutine inter_func(N,M,NP,NQ,LDN,LDM,LDNP,BETA,XPLUSD,IFIXB,IFIXX,LDIFX,&
IDEVAL,F,FJACB,FJACD,ISTOP)
use REAL_PRECISION
integer :: IDEVAL,ISTOP,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ
real (kind=R8) :: BETA(1:NP),F(1:LDN,1:NQ),FJACB(1:LDN,1:LDNP,1:NQ), &
FJACD(1:LDN,1:LDM,1:NQ),XPLUSD(1:LDN,1:M)
integer :: IFIXB(1:NP),IFIXX(1:LDIFX,1:M)
integer(c_int) :: inter_IDEVAL,inter_ISTOP,inter_LDIFX,inter_LDM,&
inter_LDN,inter_LDNP,inter_M,inter_N,inter_NP,inter_NQ
real (c_double) :: inter_BETA(1:NP),inter_F(1:LDN,1:NQ),&
inter_FJACB(1:LDN,1:LDNP,1:NQ),&
inter_FJACD(1:LDN,1:LDM,1:NQ),inter_XPLUSD(1:LDN,1:M)
integer(c_int) :: inter_IFIXB(1:NP),inter_IFIXX(1:LDIFX,1:M)
inter_IDEVAL = IDEVAL
inter_ISTOP = ISTOP
inter_LDIFX = LDIFX
inter_LDM = LDM
inter_LDN = LDN
inter_LDNP = LDNP
inter_M = M
inter_N = N
inter_NP = NP
inter_NQ = NQ
!!!!REAL array
inter_BETA = BETA
inter_F = F
inter_FJACB = FJACB
inter_FJACD = FJACD
inter_XPLUSD = XPLUSD
!!!!INTEGER array
inter_IFIXB = IFIXB
inter_IFIXX = IFIXX
end subroutine inter_func
end subroutine wrapper_ODR
#include "stdafx.h"
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
//宣告ODR呼叫接口
extern "C" {
void wrapper_ODR(void(*)(int*, int*, int*, int*, int*, int*, int*, \
double [2] ,double [1][4],int [2], int [1][4],int* ,\
int*, double [1][4],double [1][2][4],double [1][1][4] ,int* ),\
int N,int M,int NP,int NQ,double BETA[],double Y[][4],double X[][4],\
double DELTA[][4],double WE[],double WD[],int IFIXB[],int IFIXX[],\
int *JOB,int *NDIGIT,double *TAUFAC,double *SSTOL, double *PARTOL,\
int *MAXIT, int *IPRINT, int *LUNERR, int *LUNRPT,double STPB[],\
double STPD[], double SCLB[], double SCLD[], double WORK[], double IWORK[],\
int *INFO, double LOWER[], double UPPER[]);
}
//自定义函数
void FCN(int *N, int *M, int *NP, int *NQ, int *LDN, int *LDM, int *LDNP,\
double BETA[2], double XPLUSD[1][4], int IXIFB[], int IFIXX[1][4],int *LDIFX,\
int *IDEVAL, double F[1][4], double FJACB[1][2][4], double FJACD[1][1][4],int *ISTOP){
*ISTOP = 2;
if (fmod(*IDEVAL,10)!=0) {
for (int i = 0; i < 4; i++) {
F[0] = BETA[0] * XPLUSD[0] + BETA[1]; //Fitting model BETA[0]*x+BETA[1]
}
}
}
int main(){
int NP = 2, N = 4, M = 1, NQ = 1;
double BETA[2] = { 2.0, 0.5 };
double X[1][4] = { 0.0, 1.0, 2.0, 3.0 };
double Y[1][4] = { 2.0, 5.0, 8.0, 11.0 };
wrapper_ODR(FCN, N, M, NP, NQ, BETA, Y, X, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, \
NULL, NULL, NULL, NULL, NULL, NULL, NULL,NULL, NULL, NULL, NULL, NULL, NULL, \
NULL, NULL);
system("pause");
return 0;
}
2016-12-04.png (25.22 KB, 下载次数: 293)
结果
pasuka 发表于 2016-12-6 08:32
为啥不考虑Python，直接吃现成饭
Orthogonal distance regression (scipy.odr) — SciPy v0.18.1 Reference ...
#include "stdafx.h"
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
//宣告ODR呼叫接口
extern "C" {
void wrapper_ODR( void f(int*,int*,double*) , int , int , double * , double * );
}
//自定义函数
void FCN( int *N , int *M , double *delta ){
for (int i = 0; i < *M; i++) {
for (int j = 0; j < *N; j++) {
*((double*)delta + *N*i + j) = 12.0;
}
}
}
int main(){
const int N=4 ;
const int M=3;
double delta[M][N] , arr[N][M];
for (int i = 0; i < N; i++) {
for (int j = 0; j < M; j++) {
delta[j] = 0.0;
arr[j] = 0.0;
}
}
wrapper_ODR( FCN , N , M , &delta[0][0] , &arr[0][0] );
for (int i = 0; i < N; i++) {
for (int j = 0; j < M; j++) {
cout << delta[j] << endl ;
cout << arr[j] << endl;
}
}
system("pause");
return 0;
}
subroutine wrapper_ODR(FCN,N,M,Delta,Arr) bind(C, name='wrapper_ODR')
!DEC$ ATTRIBUTES DLLEXPORT :: wrapper_ODR
use , intrinsic :: iso_c_binding
implicit none
integer(c_int),value :: N,M
type(C_PTR) , value :: arr
real(c_double) , pointer :: fp(:,:) !//传递数组可以用 C_PTR 转成fortran指针
real(c_double) :: Delta(N,M) !//也可以用自动数组
!如果 fcn 是 C 语言提供的，则需要写以下interface，否则不需要
interface
subroutine FCN(N,M,delta) bind(C)
import
integer(c_int) :: n,m
real (c_double) :: delta(n,m)
end subroutine
end interface
call c_f_pointer( arr , fp , shape=[m,n]) !//c指针转成fortran指针
call ODR(fcn,delta,fp)
fp = 100.0d0;
contains
Subroutine ODR( f , d , a )
Procedure(FCN) :: f
real(kind=8) :: d(:,:) , a(:,:)
call f( size(d,dim=1) , size(d,dim=2) , d )
end subroutine ODR
end subroutine wrapper_ODR
Radiosan 发表于 2016-12-6 12:37
您好! 在一年前我老师丢给我ODR要我能让他在VC++下可以用，在这期间我发现可以用Python，我也进入Python ...
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) | Powered by Discuz! X3.2 |
