[Fortran] 纯文本查看 复制代码
!Fortran的算法，Guass_Jordan法解线性方程
subroutine gauss_jordan(pa,pb,pc,n) Bind(C,Name="gauss_jordan")
use , intrinsic :: ISO_C_Binding
type(C_PTR) , value :: pa , pb , pc
integer , value :: n
real , pointer :: a(:,:),b(:),c(:)
call c_f_pointer( pa , a , [n,n])
call c_f_pointer( pb , b , [n])
call c_f_pointer( pc , c , [n])
print *,"a=",a
print *,"b=",b
print *,"c=",c
call up(a,b,n)
call up(a,b,n)
call low(a,b,n)
forall(i=1:n)c(i)=b(i)/a(i,i)
end subroutine gauss_jordan
subroutine up(a,b,n)
real a(n,n),b(n)
do i=1,n-1
do j=i+1,n
p=a(j,i)/a(i,i)
a(j,i:n)=a(j,i:n)-a(i,i:n)*p
b(j)=b(j)-b(i)*p
end do
end do
end subroutine up
subroutine low(a,b,n)
real a(n,n),b(n)
do i=n,2,-1
do j=i-1,1,-1
p=a(j,i)/a(i,i)
a(j,i:n)=a(j,i:n)-a(i,i:n)*p
b(j)=b(j)-b(i)*p
end do
end do
end subroutine low
[C++] 纯文本查看 复制代码
extern "C" void gauss_jordan( float *a , float *b , float *c , int n );
const int N = 5;
int main(int argc, char *argv[]){
float *a , *b , *c;
a = new float[N*N];
b = new float[N];
c = new float[N];
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
a[i*N+j] = (i+1)*j;
}
b = i*10-3;
c = i*-30+23;
}
gauss_jordan( a , b , c , N );
delete[] a;
delete[] b;
delete[] c;
return 0;
}