论坛 › 程序语言综合 › 求计算方阵行列式的C代码!

bobnie 发表于 2008-6-8 21:15:06

求计算方阵行列式的C代码!

目前考虑的算法有两种:
1.按行或者按列展开法;
2.对角化法(化为上三角阵即可);
哪位大侠能赐教?如果使用第二种算法,请同时输出上三角矩阵.谢谢......

albertshu 发表于 2008-6-9 11:40:35

这是以前刚学c时候写的结构不怎么好给楼主参考吧

这里采用的计算方法是对角化法 源程序：

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define SIZE sizeof(double)
typedef double real;
int Nn;
real *p;

void main()
{
        void inputm(real *,int);
        void printm(real *,int);
        real jisuan(int);
        real *definem(int);

        int n;
        real answer;

        printf("****determinant solver****\n");
        printf("                     by dragon\n");
        do {
                printf("input n please:");
                scanf("%d",&n);
           }while (n<=0);
        Nn=n;

        p=definem(n);
        inputm(p,n);
        answer=jisuan(n);

        printf("the answer is %.2f\n",answer);

        system("pause");
}

real *definem(int n)
{
        real *pt;
        pt=(real *)malloc(n*n*SIZE);
        return(pt);
}

void inputm(real *pt,int n)
{
        int i,j; double k;

        for(i=1;i<=n;i++)
        {
                for(j=1;j<=n;j++)
                {
                        printf("a%d%d=",i,j);
                        scanf("%lf",&k);
                        *pt=(real)k;
                        pt++;
                        printf("");
                }
                printf("\n");
        }
}

real jisuan(int n)
{
        real a(int,int);
        int swap(int);
        int i,j;
        real r,re,k=1;

        if(n==1)return (*p);
        if(a(n,n)==0)
                {
                        k=-1;
                        if(swap(n)==0) return 0;
                }

        for(i=1;i<=(n-1);i++)
        {
                r=a(i,n)/a(n,n);
                for(j=1;j<=(n-1);j++)
                        *(p+(i-1)*Nn+j-1)=a(i,j)-r*a(n,j);

        }

        k=k*a(n,n);
        re=k*jisuan(n-1);
        return re;
}

int swap(int n)
{
        real a(int,int);
        int i,j;
        real t;
        for(i=0;i<n;i++)
        {
                if(*(p+i*Nn+n-1)!=0)
                {
                        for(j=0;j<=n;j++)
                        {
                                t=a(i+1,j+1);
                                *(p+i*Nn+j)=a(n,j+1);
                                *(p+(n-1)*Nn+j)=t;
                        }
                        return 1;
                }
        }
        return 0;
}

real a(int x,int y)
{
        real t;
        t=*(p+(x-1)*Nn+y-1);
        return(t);
}

dragon_ 发表于 2008-6-9 14:07:35

哈哈，这个源代码我有。数值分析课的各类算法从头写到尾。:lol
最好使用LU分解法，数值稳定性强，最后输出下三角矩阵的对角线元素乘积即可。

bobnie 发表于 2008-6-9 16:23:38

跪求二楼大侠给出关键算法的注释,小弟数学能力有些差,

bobnie 发表于 2008-6-9 16:38:10

回复 2# 的帖子

跪求二楼大侠给出关键算法的注释,小弟数学能力有些差,

dragon_ 发表于 2008-6-16 01:32:02

对矩阵进行LU分解的子程序如下：void LUDecomposition(double matrix, int rank)
{
int i, j, k, t;
double tmp;
   
    for(k = 0; k < rank; k++)
    {
      for(i = k; i < rank; i++)
      {
            tmp = 0;
            for(t = 0; t < k; t++)tmp+= matrix * matrix;
            matrix-= tmp;
      }

      if(k < rank - 1)
      {
            for(j = k + 1; j < rank; j++)
            {
                tmp = 0;
                for(t = 0; t < k; t++)tmp+= matrix * matrix;
                matrix = (matrix - tmp) / matrix;
            }
      }
      else{}
    }
   
}matrix为所求矩阵，rank为矩阵的阶。
最后计算matrix主对角线上元素的积就是矩阵的行列式了。算这个应该没难度吧。
就是这么简单，不像2楼的代码那样繁琐。并且速度快，数值稳定性强。但是使用这个算法时要注意，矩阵matrix的所有顺序主子阵必须都是非奇异的。否则还是要使用Gauss消去法。
关于矩阵的LU分解，请参看任意一本数值分析教材。

[ 本帖最后由 dragon_ 于 2008-6-16 01:51 编辑 ]

dragon_ 发表于 2008-6-16 01:43:53

突然发现2楼的代码里有这么两行……
      printf("****determinant solver****\n");
      printf("                     by dragon\n");
by dragon……
囧死…………………:L
那不是我写的……:L

bobnie 发表于 2008-6-16 17:59:31

回复 6# 的帖子

狂谢大哥,以后多多向你学习......

bobnie 发表于 2008-6-24 14:31:29

高斯消去法的程序已经调试成功，有需要的请留下Email

dragon_ 发表于 2008-6-24 16:35:04

说实话，估计EM学数学的人不多，学过数值分析的人就更少。
把用雅各比法求实对称矩阵的特征值和特征向量的C代码
和求矩阵的协方差矩阵的C代码
给我一份吧。谢谢。:loveliness:

[ 本帖最后由 dragon_ 于 2011-1-14 19:54 编辑 ]

twl859588 发表于 2008-7-5 22:15:48

C用来干这个有点晕，好像MATLAB更好用...

查看完整版本: 求计算方阵行列式的C代码!