加载中…
个人资料
jackiewu
jackiewu
  • 博客等级:
  • 博客积分:0
  • 博客访问:107
  • 关注人气:11
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

C语言矩阵求逆

(2009-09-21 12:06:25)
标签:

c语言

杂谈

   根据代数里面的知识,可以使用伴随矩阵也可以使用初等行变换来解求解,但是这样如果矩阵的维数较大的时候,使用这种方法,矩阵的维数变大时,计算量急剧的变大,计算时间和使用内存也会按着指数急剧上升,这样的算法的生命力不行。

  使用以下这种算法的计算量和使用内存不会发生急剧的变化,特别是矩阵在维数大的时候。


高斯-约旦法(全选主元)求逆的步骤如下:

首先,对于 k 从 0 到 n - 1 作如下几步:

从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素,并记住次元素所在的行号和列号,在通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。
m(k, k) = 1 / m(k, k)
m(k, j) = m(k, j) * m(k, k),j = 0, 1, ..., n-1;j != k
m(i, j) = m(i, j) - m(i, k) * m(k, j),i, j = 0, 1, ..., n-1;i, j != k
m(i, k) = -m(i, k) * m(k, k),i = 0, 1, ..., n-1;i != k
最后,根据在全选主元过程中所记录的行、列交换的信息进行恢复,恢复的原则如下:在全选主元过程中,先交换的行(列)后进行恢复;原来的行(列)交换用列(行)交换来恢复。

 

#include"stdio.h"
#include"malloc.h"
#include"math.h"    //数学函数
void main()
{ int inv(double *p,int n);
 double a[4][4]={{1,2,0,0},{2,5,0,0},{0,0,3,0},{0,0,0,1}},*ab;
 ab=a[0];
 int n=4,i=0,j;
 i=inv(ab,n);  //调用矩阵求逆
 if(i!=0)      //如果返回值不是0
  for(i=0;i<n;i++)   //输出结果
  { putchar('\n');
    for(j=0;j<n;j++)
     printf("%f  ",a[i][j]);
  }
}
int inv(double *p,int n)
{
  void swap(double *a,double *b);
  int *is,*js,i,j,k,l;
  for(i=0;i<n;i++)
   { putchar('\n');
     for(j=0;j<n;j++)
     printf("%f  ",*(p+i*n+j));
     }
  puts("\n\n\n\n");
  double temp,fmax;
  is=(int *)malloc(n*sizeof(int));
  js=(int *)malloc(n*sizeof(int));
  for(k=0;k<n;k++)
   {
      fmax=0.0;
      for(i=k;i<n;i++)
       for(j=k;j<n;j++)
        { temp=fabs(*(p+i*n+j));//找最大值
          if(temp>fmax)
           { fmax=temp;
             is[k]=i;js[k]=j;
            }
         }
     if((fmax+1.0)==1.0)
       free(is);free(js);
          printf("no inv");
          return(0);
       }
    if((i=is[k])!=k)
      for(j=0;j<n;j++)
        swap(p(k*n+j),p(i*n+j));//交换指针
   if((j=js[k])!=k)
     for(i=0;i<n;i++)
        swap(p(i*n+k),p(i*n+j));  //交换指针
   p[k*n+k]=1.0/p[k*n+k];
   for(j=0;j<n;j++)
     if(j!=k)
        p[k*n+j]*=p[k*n+k];
   for(i=0;i<n;i++)
      if(i!=k)
        for(j=0;j<n;j++)
          if(j!=k)
             p[i*n+j]=p[i*n+j]-p[i*n+k]*p[k*n+j];
   for(i=0;i<n;i++)
     if(i!=k)
       p[i*n+k]*=-p[k*n+k];
 }
 for(k=n-1;k>=0;k--)
   {
     if((j=js[k])!=k)
        for(i=0;i<n;i++)
          swap((p+j*n+i),(p+k*n+i));
     if((i=is[k])!=k)
        for(j=0;j<n;j++)
          swap((p+j*n+i),(p+j*n+k);
  }
  free(is);
  free(js);
  return 1;

}
void swap(double *a,double *b)
double c;
   c=*a;
   *a=*b;
   *b=c;
}

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有