Chinaunix首页 | 论坛 | 博客
  • 博客访问: 543808
  • 博文数量: 92
  • 博客积分: 2511
  • 博客等级: 少校
  • 技术积分: 932
  • 用 户 组: 普通用户
  • 注册时间: 2008-10-19 10:10
文章分类
文章存档

2011年(6)

2010年(27)

2009年(37)

2008年(22)

我的朋友

分类: LINUX

2009-12-28 23:12:42

    幂法的原理大概每本数值分析的书上都会讲吧,这里参考的是清华大学出版社的《数值分析》第四版,在此就不赘述原理了。


/*Author:RainMan
Date:2009-12-28
请使用符合C99标准的编译器编译以下程序
输入文件由以下格式构成,第一行是矩阵的维数
以下各行是矩阵的值,比如
3 3
1.0 1.0 0.5
1.0 1.0 0.25
0.5 0.25 2.0
*/

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

void multiple(double**A,double *V,int dim_x,int dim_y);
double max(double *V,int dim);
void div_matrix(double *V,int dim,double m);

int main(){
    FILE *file = fopen("dengjin.txt","r");
    int dim_x,dim_y;
    double **A,*V;
    double miu0 = 0,miu1 = 10;/*just make sure to enter the loop*/
    fscanf(file,"%d %d",&dim_x,&dim_y);

/*load in data*/
    A = (double **)malloc(sizeof(double *)*dim_x);
    V = (double *)malloc(sizeof(double)*dim_y);
    for(int i=0;i<dim_x;i++)
        A[i] = (double *)malloc(sizeof(double)*dim_y);

    for(int i=0;i<dim_x;i++)
        for(int j=0;j<dim_y;j++)
            fscanf(file,"%lf",&A[i][j]);

    for(int i=0;i<dim_y;i++)
        V[i] = 1;/*initialing a vector with any value*/

    while(fabs(miu1-miu0) >= 1E-8){
        multiple(A,V,dim_x,dim_y);
        miu0 = miu1;
        miu1 = max(V,dim_y);
        div_matrix(V,dim_y,miu1);
        for(int i=0;i<dim_y;i++)
            printf("%10.8lf ",V[i]);
        puts("");
    }
    printf("Eigenvalue: %10.8lf\n",miu1);
    //delocating


    free(V);
    for(int i=0;i<dim_y;i++)
        free(A[i]);
    free(A);

    return EXIT_SUCCESS;
}

void multiple(double**A,double *V,int dim_x,int dim_y){
    double *tmp = (double *)malloc(sizeof(double)*dim_y);
    for(int i=0;i<dim_y;i++)
        tmp[i] = 0;

    for(int i=0;i<dim_x;i++)
        for(int j=0;j<dim_y;j++)
            tmp[i] += A[i][j]*V[j];

    for(int i=0;i<dim_y;i++)
        V[i] = tmp[i];
    free(tmp);
}

double max(double *V,int dim){
    double tmp = V[0];
    for(int i=1;i<dim;i++)
            if(fabs(V[i]) > fabs(tmp))
                tmp = V[i];
    return tmp;
}

void div_matrix(double *V,int dim,double m){
    for(int i=0;i<dim;i++)
        V[i] /= m;
}


阅读(7087) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~