#include <math.h>
#include <stdio.h>
int brinv(double a[], int n)
{ int *is,*js,i,j,k,l,u,v
double d,p
is=malloc(n*sizeof(int))
js=malloc(n*sizeof(int))
for (k=0k<=n-1k++)
{ d=0.0
for (i=ki<=n-1i++)
for (j=kj<=n-1j++)
{ l=i*n+jp=fabs(a[l])
if (p>d) { d=pis[k]=ijs[k]=j}
}
if (d+1.0==1.0)
{ free(is)free(js)printf("err**not inv\n")
return(0)
}
if (is[k]!=k)
for (j=0j<=n-1j++)
{ u=k*n+jv=is[k]*n+j
p=a[u]a[u]=a[v]a[v]=p
}
if (js[k]!=k)
for (i=0i<=n-1i++)
{ u=i*n+kv=i*n+js[k]
p=a[u]a[u]=a[v]a[v]=p
}
l=k*n+k
a[l]=1.0/a[l]
for (j=0j<=n-1j++)
if (j!=k)
{ u=k*n+ja[u]=a[u]*a[l]}
for (i=0i<=n-1i++)
if (i!=k)
for (j=0j<=n-1j++)
if (j!=k)
{ u=i*n+j
a[u]=a[u]-a[i*n+k]*a[k*n+j]
}
for (i=0i<=n-1i++)
if (i!=k)
{ u=i*n+ka[u]=-a[u]*a[l]}
}
for (k=n-1k>=0k--)
{ if (js[k]!=k)
for (j=0j<=n-1j++)
{ u=k*n+jv=js[k]*n+j
p=a[u]a[u]=a[v]a[v]=p
}
if (is[k]!=k)
for (i=0i<=n-1i++)
{ u=i*n+kv=i*n+is[k]
p=a[u]a[u]=a[v]a[v]=p
}
}
free(is)free(js)
return(1)
}
void brmul(double a[], double b[],int m,int n,int k,double c[])
{ int i,j,l,u
for (i=0i<=m-1i++)
for (j=0j<=k-1j++)
{ u=i*k+jc[u]=0.0
for (l=0l<=n-1l++)
c[u]=c[u]+a[i*n+l]*b[l*k+j]
}
return
}
int main()
{ int i,j
static double a[4][4]={ {0.2368,0.2471,0.2568,1.2671},
{1.1161,0.1254,0.1397,0.1490},
{0.1582,1.1675,0.1768,0.1871},
{0.1968,0.2071,1.2168,0.2271}}
static double b[4][4],c[4][4]
for (i=0i<=3i++)
for (j=0j<=3j++)
b[i][j]=a[i][j]
i=brinv(a,4)
if (i!=0)
{ printf("MAT A IS:\n")
for (i=0i<=3i++)
{ for (j=0j<=3j++)
printf("%13.7e ",b[i][j])
printf("\n")
}
printf("\n")
printf("MAT A- IS:\n")
for (i=0i<=3i++)
{ for (j=0j<=3j++)
printf("%13.7e ",a[i][j])
printf("\n")
}
printf("\n")
printf("MAT AA- IS:\n")
brmul(b,a,4,4,4,c)
for (i=0i<=3i++)
{ for (j=0j<=3j++)
printf("%13.7e ",c[i][j])
printf("\n")
}
}
}
这个行么？不知道是不是你要的，希望对你有用
//源程序如下#include<stdio.h>#include<conio.h>
#include<string.h>
#include<iostream.h>
#include<stdlib.h>
#include<math.h>
#define max 100void inputstyle(int *) //输入函数
void input(int **,int) //输入函数
long danx(int **,int)
int sgnx(int)
void martx(int **,int)int main(void)
{
int style=0,i=0
int matrix[max][max],*p[max]
for(i=0i<maxi++)*(p+i)=matrix[i] //*(p+i)是指针，指向第i个字符串
char exit1=' '
while(exit1!='E'&&exit1!='e'){ printf("求n阶矩阵的逆\n") inputstyle(&style)
input(p,style)
printf("原矩阵为:\n")
for(i=0i<stylei++){
for(int j=0j<stylej++){
printf("%4d",matrix[i][j])
}
printf("\n")
}
martx(p,style)
printf("\n")
printf("Exit=e Continue=Press any key\n")
cin>>exit1
fflush(stdin)
printf("\n\n")}
return(0)
} void input(int **p,int n){
for(int i=0i<ni++){
for(int j=0j<nj++){
printf("输入矩阵(%d行,%d列)元素:",j+1,i+1)
*(*(p+j)+i)=0
scanf("%d",*(p+j)+i)
fflush(stdin)
}
}
}void inputstyle(int *style){
do{
printf("输入矩阵n*n阶数n(0<n<%d):",max)
fflush(stdin)
scanf("%d",style)
fflush(stdin)
}while(*style<=0 &&*style>max)
}long danx(int **p,int n){
int i=0,j1=0,k1=0,j2=0,k2=0
long sum=0
int operate[max][max],*po[max]
for(i=0i<maxi++)*(po+i)=operate[i]if(n==1)return *(*(p+0)+0)
else{
for(i=0i<ni++){
for(j1=1,j2=0j1<nj1++,j2++){
k1=-1k2=-1
while(k2<n-1){
k1++
k2++
if(k1==i)k1++
*(*(po+j2)+k2)=*(*(p+j1)+k1)
}
}
/*for(int i1=0i1<n-1i1++){
for(int h1=0h1<n-1h1++){
printf("(%d,%d)%d ",i1,h1,*(*(po+h1)+i1))
}
printf("\n")
}*/
sum+=*(*(p+0)+i) * sgnx(1+i+1) * danx(po,n-1)
}
return sum
}
}int sgnx(int i){
if(i%2==0)return(1)
else return(-1)
}void martx(int **p,int n){
int i=0,j=0,j1=0,k1=0,j2=0,k2=0,num=0
int tramform[max][max]
int operate[max][max],*po[max]
for(i=0i<maxi++)*(po+i)=operate[i]
num=danx(p,n)
if(num==0)printf("矩阵不可逆\n")
else{
if(n==1)printf("矩阵的逆为: 1/%d\n",num)
else{
printf("矩阵的逆为: 系数 1/%d *\n",num)
for(i=0i<ni++){
for(j=0j<nj++){
j1=-1j2=-1
while(j2<n-1){
j1++j2++
if(j1==j)j1++ k1=-1k2=-1
while(k2<n-1){
k1++
k2++
if(k1==i)k1++
*(*(po+j2)+k2)=*(*(p+j1)+k1)
}
}
tramform[i][j]=sgnx(2+i+j) * danx(po,n-1)
}
}
for(i=0i<ni++){
for(j=0j<nj++){
printf("%4d",tramform[i][j])
}
printf("\n")
}
}
}
}
//运行结果//希望对你有帮助
这是一段f90的代码，使用是时候要加上use inv_mat，然后就可以求某方阵的逆矩阵了。module inv_mat
! Description : 计算逆矩阵
contains
subroutine inv(A,invA,N)
! Purpose : 计算逆矩阵
!-----------------------------------------------------
! Input parameters :
! 1. A: 需要求逆的矩阵
! 2. N: 矩阵的维度
! Output parameters :
! 1. invA: A的逆矩阵
implicit real*8(a-z)
integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)
E=0
!设置E为单位矩阵
do i=1,n
E(i,i)=1
end do
call mateq(A,E,invA,N,N)
end subroutine inv
subroutine mateq(A,B,X,N,M)
! Purpose : 高斯列主元消去法计算矩阵方程
! AX=B
implicit real*8(a-z)
integer::N,M,i
real*8::A(N,N),B(N,M),X(N,M)
real*8::btemp(N),xtemp(N)
do i=1,M
btemp=B(:,i)
call elgauss(A,btemp,xtemp,N)
X(:,i)=xtemp
end do
end subroutine mateq
subroutine elgauss(A,b,x,N)
! Purpose : 高斯列主元消去法
! Ax=b
implicit real*8(a-z)
integer::i,k,N
integer::id_max !主元素标号
real*8::A(N,N),b(N),x(N)
real*8::Aup(N,N),bup(N)
!Ab为增广矩阵 [Ab]
real*8::Ab(N,N+1)
real*8::vtemp1(N+1),vtemp2(N+1)
Ab(1:N,1:N)=A
Ab(:,N+1)=b
do k=1,N-1
elmax=dabs(Ab(k,k))
id_max=k
do i=k+1,n
if (dabs(Ab(i,k))>elmax) then
elmax=Ab(i,k)
id_max=i
end if
end do
vtemp1=Ab(k,:)
vtemp2=Ab(id_max,:)
Ab(k,:)=vtemp2
Ab(id_max,:)=vtemp1
!#########################################################
do i=k+1,N
temp=Ab(i,k)/Ab(k,k)
Ab(i,:)=Ab(i,:)-temp*Ab(k,:)
end do
end do
Aup(:,:)=Ab(1:N,1:N)
bup(:)=Ab(:,N+1)
call uptri(Aup,bup,x,n)
end subroutine elgauss
subroutine uptri(A,b,x,N)
! Purpose : 上三角方程组的回带方法
! Ax=b
implicit real*8(a-z)
integer::i,j,N
real*8::A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
!回带部分
do i=n-1,1,-1
x(i)=b(i)
do j=i+1,N
x(i)=x(i)-a(i,j)*x(j)
end do
x(i)=x(i)/A(i,i)
end do
end subroutine uptri
end module inv_mat
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）