论坛首页· 友情链接申请·申请版主· 广告投放· 道具中心· 设为首页· 收藏本站
发新话题
打印

熟悉雅可比法的帮忙了!!!!

熟悉雅可比法的帮忙了!!!!

程序如下:(运行出现死循环,不知哪儿出了毛病?请指点!)

#define EPSL 1.0e-12                                                            //迭代的精度水平
#define PI 3.1415926


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


void main()                                                                      //主程序开始;
{
        int size,//协方差矩阵的阶数
                row,col,k,
                p,q;//记录maxa的行与列
        float a[99][99],//存放协方差矩阵
                R[99][99]={0},
                RT[99][99],//矩阵R的转置
                S[99][99],//存放特征向量,赋初值为单位矩阵
                tan,//旋转角度的正切
                alfa,//旋转角度
                maxa=0;//非主对角线上绝对值最大元素
        FILE *fp;

        fp=fopen("sj.txt","r");
        fscanf(fp,"%d",&size);//读入协方差矩阵的阶数

        for(row=0;row<size;row++)
                for(col=0;col<size;col++)
                        fscanf(fp,"%f",&a[row][col]);

        fclose(fp);

        for(row=0;row<size;row++)//给矩阵S赋初值为单位矩阵
                for(col=0;col<size;col++)
                        if(row==col)
                                S[row][col]=1;
                        else
                                S[row][col]=0;


        while(maxa==0||maxa>EPSL)
        {
                for(row=0;row<size;row++)
                        for(col=row+1;col<size;col++)
                                if(fabs(a[row][col])>maxa)
                                {
                                        maxa=a[row][col];//找出非主对角线上绝对值最大的元素
                                        p=row;q=col;//记录maxa的行与列数
                                }

                if(a[p][p]=a[q][q])
                        tan=PI/2.0;
                else
                        tan=2*a[p][q]/(a[p][p]-a[q][q]);

                alfa=atan(tan);//求出旋转角度

                R[p][p]=cos(alfa);
                R[q][q]=cos(alfa);
                R[p][q]=-sin(alfa);
                R[q][p]=sin(alfa);//形成矩阵R

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                RT[col][row]=R[row][col];//形成R矩阵的转置矩阵RT

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                for(k=0;k<size;k++)
                                        a[row][col]+=RT[row][k]*a[k][col];      //计算RT与a的乘积,并赋给a

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                for(k=0;k<size;k++)
                                        a[row][col]+=a[row][k]*R[k][col];      //计算a与R的乘积,并赋给a

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                for(k=0;k<size;k++)
                                        S[row][col]+=S[row][k]*R[k][col];      //计算S与R的乘积,并赋给S
        }

        fp=fopen("sc1.txt","w");

        for(row=0;row<size;row++)
                for(col=0;col<size;col++)
                        fscanf(fp,"%f",a[row][col]);
        fclose(fp);
        fp=fopen("sc2.txt","w");

        for(row=0;row<size;row++)
                for(col=0;col<size;col++)
                        fscanf(fp,"%f",S[row][col]);
                fclose(fp);
}

TOP

你能不能等等,现在没时间仔细看这些代码。暑假好好讨论
因为梦想而努力,因为有你而精彩!
mail: qianzongming@gmail.com

TOP

发新话题