>下面是一个定义数组的例子,我估计你定义不了是因为在栈里面无法定义,而在堆里面是可以定义的!</P>
>#define T 100000</P>
>#include "stdlib.h"<BR>#include "math.h"<BR>#include "stdio.h"<BR>#include <string.h><BR>struct position<BR>{<BR> int x;<BR> int y;<BR>};<BR>//////////////////////////////////////////////////////////////////////////////<BR>//摘自《C常用算法程序集(第二版)》,实矩阵求逆的全选主高斯-约当法<BR>//<BR>//int brinv(a,n)<BR>//<BR>//本函数返回一个整形标志值。若返回的标志值为0,则表示矩阵A奇异,还输出信息<BR>//“err * * not inv”;若返回的标志值不为0,则表示正常返回。<BR>//形参说明:<BR>//a——双精度实型二维数组,体积为n*n。存放原矩阵A;返回时存放其逆矩阵A逆。<BR>//n——整形变量。矩阵的阶数。<BR>//////////////////////////////////////////////////////////////////////////////<BR>int brinv(double a[],int n)<BR>{ <BR> int *is,*js,i,j,k,l,u,v;<BR> double d,p;<BR> is=(int *)malloc((unsigned)n*sizeof(int));<BR> js=(int *)malloc((unsigned)n*sizeof(int));<BR> for (k=0; k<=n-1; k++)<BR> {<BR> d=0.0;<BR> for (i=k; i<=n-1; i++)<BR> for (j=k; j<=n-1; j++)<BR> {<BR> l=i*n+j; <BR> p=fabs(a[l]);<BR> if (p>d)<BR> {<BR> d=p;<BR> is[k]=i;<BR> js[k]=j;<BR> }<BR> }<BR> if (fabs(d)<0.0000001)<BR> {<BR> free(is);<BR> free(js);<BR> printf("err**not inv\n");<BR> return(0);<BR> }<BR> if (is[k]!=k)<BR> {<BR> for (j=0; j<=n-1; j++)<BR> { <BR> u=k*n+j;<BR> v=is[k]*n+j;<BR> p=a;<BR> a=a[v]; <BR> a[v]=p;<BR> }<BR> }<BR> if (js[k]!=k)<BR> {<BR> for (i=0; i<=n-1; i++)<BR> {<BR> u=i*n+k; <BR> v=i*n+js[k];<BR> p=a; <BR> a=a[v];<BR> a[v]=p;<BR> }<BR> }<BR> l=k*n+k;<BR> a[l]=1.0/a[l];<BR> for (j=0; j<=n-1; j++)<BR> {<BR> if (j!=k)<BR> {<BR> u=k*n+j;<BR> a=a*a[l];<BR> }<BR> }<BR> for (i=0; i<=n-1; i++)<BR> {<BR> if (i!=k)<BR> {<BR> for (j=0; j<=n-1; j++)<BR> {<BR> if (j!=k)<BR> {<BR> u=i*n+j;<BR> a=a-a[i*n+k]*a[k*n+j];<BR> }<BR> }<BR> }<BR> }<BR> for (i=0; i<=n-1; i++)<BR> {<BR> if (i!=k)<BR> {<BR> u=i*n+k;<BR> a=-a*a[l];<BR> }<BR> }<BR> }<BR> for (k=n-1; k>=0; k--)<BR> {<BR> if (js[k]!=k)<BR> {<BR> for (j=0; j<=n-1; j++)<BR> {<BR> u=k*n+j; <BR> v=js[k]*n+j;<BR> p=a; <BR> a=a[v];<BR> a[v]=p;<BR> }<BR> }<BR> if (is[k]!=k)<BR> {<BR> for (i=0; i<=n-1; i++)<BR> {<BR> u=i*n+k;<BR> v=i*n+is[k];<BR> p=a; <BR> a=a[v];<BR> a[v]=p;<BR> }<BR> }<BR> }<BR> free(is); <BR> free(js);<BR> return(1);<BR>}</P>" {5 C! c7 l, w4 R: q- B
>//////////////////////////////////////////////////////////////////////////////<BR>//摘自《C常用算法程序集(第二版)》,实矩阵相乘<BR>//<BR>//void brmul(a,b,m,n,k,c)<BR>//<BR>//形参说明:<BR>//a——双精度实型二维数组,体积为m*n。存放矩阵A的元素。<BR>//b——双精度实型二维数组,体积为n*k。存放矩阵B的元素。<BR>//m——整形变量。矩阵A的行数,也是乘积矩阵C=AB的行数。<BR>//n——整形变量。矩阵A的列数,也是矩阵B的行数。<BR>//k——整形变量。矩阵B的列数,也是乘积矩阵C=AB的列数。<BR>//c——双精度实型二维数组,体积为m*k。返回乘积矩阵C=AB的元素。<BR>//////////////////////////////////////////////////////////////////////////////<BR>void brmul(double a[],double b[],int m,int n,int k,double c[])<BR>{<BR> int i,j,l,u;<BR> for (i=0; i<=m-1; i++)<BR> {<BR> for (j=0; j<=k-1; j++)<BR> {<BR> u=i*k+j;<BR> c=0.0;<BR> for (l=0; l<=n-1; l++) c=c+a[i*n+l]*b[l*k+j];<BR> }<BR> }<BR> return;<BR>}</P>
>//////////////////////////////////////////////////////////////////////////////<BR>void main()<BR>{<BR> FILE *para_fp,*in_fp,*out_fp,*fp;<BR> int n;// 特征空间的大小,即共n个特征<BR> int lines,columns;//对应数据场的行列值<BR> int m;//分类数<BR> int *nl;//一维数组,存放每类的样本点数量<BR> struct position **style;//二维数组,存放每类的样本点坐标<BR> double *xp,**x;//二维数组x,存放训练样本数据的数组;<BR> double *xtp,**xt;//数组x的逆,存放训练样本数据的数组;<BR> double *up,**u;//二维数组,存放每一类的均值<BR> double ***sigma,*sigmap,**sigmapp;//三维数组,存放每类的散度矩阵<BR> double ***anti_sigma,*anti_sigmap,**anti_sigmapp;//三维数组,存放每类的散度矩阵的逆矩阵<BR> double **sigma_temp,*sigma_tempp,**anti_sigma_temp,*anti_sigma_tempp;//临时用的矩阵<BR> int return_brinv;//矩阵求逆返回值,标示矩阵是否奇异<BR> double *g;//数组,存放像元x 与每一类的广义距离<BR> double min_g;<BR> int min_i;//存放最大g的i<BR> int datatype;//对应每一数据类型所占的字节数<BR> int i,j,k;//循环变量<BR> long count;//循环变量,扫描整个数据场<BR> char d[10];//原数据场的数据类型<BR> //以下为临时变量<BR> struct position style_temp;<BR> char pels_temp;<BR> char temp;<BR> short temp_0;<BR> int temp_1;<BR> long temp_2;<BR> float temp_3;<BR> double temp_4;</P>
> double g_temp0[1],*g_temp1p,**g_temp1;<BR>//////////////////////////////////////////////////////////////////////////////////////////////<BR>// 变量说明<BR>//<BR>// nl每类样本点数量 style每类的样本点坐标<BR>// | |<BR>// nl[0] style[0][0]---style[0][(nl[0])]<BR>// nl[1] style[1][0]---style[0][(nl[1])]<BR>//<BR>// …… ……<BR>//<BR>// nl[m] style[m][0]---style[0][(nl[m])] <BR>//<BR>//sigma[m][n][n]:存放所有类的sigma矩阵;<BR>//anti_sigma[m][n][n]:存放所有类的sigma矩阵的逆矩阵;<BR>//<BR>////////////////////////////////////////////////////////////////////////////////////////////////</P>! {7 O1 _: }( S; I' M
><BR>////////////读取特征空间的大小//////////////////////////////////<BR> if((para_fp=fopen("guangyi.par","rt"))==NULL)<BR> { <BR> return;<BR> }<BR> fscanf(para_fp,"%d",&n);//读取特征空间的大小<BR> fscanf(para_fp,"%d",&m);//预分类数<BR> //band=(int *)malloc((unsigned)n*sizeof(int));//存取n个波段号<BR> //for(i=0;i<n;i++) band=0;<BR> //for(i=0;i<n;i++)<BR> //{<BR> // fscanf(para_fp,"%d",&(band));<BR> //}<BR> fclose(para_fp);</P>
> //for(i=0;i<n;i++)<BR> //{<BR> //}<BR>////////////读取数据场的行列值//////////////////////////////////<BR> if((para_fp=fopen("sourceimg.par","rt"))==NULL)<BR> { <BR> return;<BR> }<BR> fscanf(para_fp,"%d,%d",&lines,&columns);<BR> fclose(para_fp);</P>2 L4 {, }* q9 z
> printf("lines=%d,columns=%d\n",lines,columns);<BR> ///////////////////////////////////<BR> if((fp=fopen("leixing.par","rt"))==NULL)<BR> { <BR> return;<BR> }<BR> i=0;<BR> while(!feof(fp))<BR> { <BR> d=fgetc(fp);<BR> // printf("%c\n",d);<BR> i=i+1;<BR> }<BR> d[i-1]='\0';<BR> //fscanf(fp,"%s",d);</P>
> if(strcmp(d,"char")==0) <BR> {<BR> datatype=1;<BR> }<BR> if(strcmp(d,"short")==0)<BR> {<BR> datatype=2;<BR> }<BR> if(strcmp(d,"int")==0)<BR> {<BR> datatype=4;<BR> }<BR> if(strcmp(d,"long")==0)<BR> {<BR> datatype=4;<BR> }<BR> if(strcmp(d,"float")==0)<BR> {<BR> datatype=4;<BR> }<BR> if(strcmp(d,"double")==0)<BR> {<BR> datatype=8;<BR> }<BR> printf("datatype=%d\n",datatype);<BR> printf("%s",d);<BR>////////////读取样本数和样本区//////////////////////////////////<BR> style=(struct position **)malloc((unsigned)m*sizeof(struct position *));<BR> nl=(int*)malloc((unsigned)m*sizeof(int));</P>
> if((para_fp=fopen("style.par","rt"))==NULL)<BR> { <BR> free(style);<BR> free(nl);<BR> return;<BR> }<BR> for(i=0;i<m;i++)<BR> {<BR> fscanf(para_fp,"%d",&(nl));<BR> style=(struct position *)malloc((unsigned)(nl)*sizeof(struct position));<BR> for(j=0;j<(nl);j++) <BR> {<BR> style_temp.x=0;<BR> style_temp.y=0;<BR> style[j]=style_temp;<BR> }<BR> for(j=0;j<(nl);j++)<BR> {<BR> // 读入每一类的样本点的坐标<BR> fscanf(para_fp,"%d,%d",&(style_temp.x),&(style_temp.y));<BR> style[j]=style_temp;<BR> }<BR> }<BR> fclose(para_fp);</P>
> for(i=0;i<m;i++)<BR> {<BR> printf("\nClass %3d:%3d\n",i,nl);<BR> for(j=0;j<(nl);j++)<BR> {<BR> style_temp=style[j];<BR> // printf("(%d,%d)",style_temp.x,style_temp.y);<BR> }<BR> }<BR>////////////计算m个类别训练组的u,sigma,sigma逆//////////////////////////////////<BR> if((in_fp=fopen("sourceimg.dat","rb"))==NULL)<BR> { <BR> for(k=0;k<m;k++)<BR> {<BR> free(style[k]);<BR> }<BR> free(style);<BR> free(nl);<BR> return;<BR> }</P># e# ~5 g& U8 ^7 j9 x
> up=(double *)malloc((unsigned)n*m*sizeof(double));<BR> u=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> u[k]=&(up[k*m]);<BR> }<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<m;j++)<BR> {<BR> u[k][j]=0.0;<BR> }<BR> }</P>
> sigmap=(double *)malloc((unsigned)m*n*n*sizeof(double));<BR> sigmapp=(double **)malloc((unsigned)m*n*sizeof(double *));<BR> sigma=(double ***)malloc((unsigned)m*sizeof(double **));</P>5 t' ?5 S, ^! P
> for(i=0;i<m;i++) <BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> sigmapp[i*n+j]=&(sigmap[i*n*n+j*n]);<BR> }<BR> sigma=&(sigmapp[i*n]);<BR> }</P>
> anti_sigmap=(double *)malloc((unsigned)m*n*n*sizeof(double));<BR> anti_sigmapp=(double **)malloc((unsigned)m*n*sizeof(double *));<BR> anti_sigma=(double ***)malloc((unsigned)m*sizeof(double **));</P>4 q3 `- x E. g' t8 I! P
> for(i=0;i<m;i++) <BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> anti_sigmapp[i*n+j]=&(anti_sigmap[i*n*n+j*n]);<BR> }<BR> anti_sigma=&(anti_sigmapp[i*n]);<BR> }</P>
> //对每一类进行循环,计算计算m个类别训练组的u,sigma,sigma逆<BR> for(i=0;i<m;i++)<BR> {<BR> sigma_tempp=(double *)malloc((unsigned)n*n*sizeof(double));<BR> sigma_temp=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> sigma_temp[k]=&(sigma_tempp[k*n]);<BR> }</P>* w, o; i" b3 V, K6 l
> anti_sigma_tempp=(double *)malloc((unsigned)n*n*sizeof(double));<BR> anti_sigma_temp=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> anti_sigma_temp[k]=&(anti_sigma_tempp[k*n]);<BR> }</P>. @3 w7 f/ x0 l2 l" j4 U2 y
> xp=(double *)malloc((unsigned)n*nl*sizeof(double));<BR> x=(double **)malloc((unsigned)n*sizeof(double*)); //训练样本<BR> for(k=0;k<n;k++)<BR> {<BR> x[k]=&(xp[k*nl]);<BR> }<BR> xtp=(double *)malloc((unsigned)nl*n*sizeof(double));<BR> xt=(double **)malloc((unsigned)nl*sizeof(double*)); //训练样本的转置<BR> for(j=0;j<nl;j++)<BR> {<BR> xt[j]=&(xtp[j*n]);<BR> }</P>
><BR>///////////////////读取矩阵X///////////////////////////////////<BR> for(j=0;j<nl;j++)<BR> { <BR> style_temp=style[j];<BR> for(k=0;k<n;k++)<BR> {<BR> //style_temp=style[j];<BR> if(k==0) <BR> { <BR> fseek(in_fp,(style_temp.x*lines+style_temp.y)*datatype,SEEK_SET);<BR> }<BR> else<BR> {<BR> fseek(in_fp,(lines*columns-1)*datatype,SEEK_CUR);<BR> }<BR> /////////////////从数据场文件中读入训练样本的数据<BR> if(1==datatype)<BR> {<BR> fread(&temp,datatype,1,in_fp);<BR> x[k][j]=(double)temp;<BR> }<BR> if(2==datatype)<BR> {<BR> fread(&temp_0,datatype,1,in_fp);<BR> x[k][j]=(double)temp_0;<BR> }<BR> if(strcmp(d,"int")==0)<BR> {<BR> fread(&temp_1,datatype,1,in_fp);<BR> x[k][j]=(double)temp_1;<BR> }<BR> if(strcmp(d,"long")==0)<BR> {<BR> fread(&temp_2,datatype,1,in_fp); <BR> x[k][j]=(double)temp_2;<BR> } <BR> if(strcmp(d,"float")==0)<BR> {<BR> fread(&temp_3,datatype,1,in_fp); <BR> x[k][j]=(double)temp_3;<BR> }<BR> if(8==datatype)<BR> {<BR> fread(&temp_4,datatype,1,in_fp); <BR> x[k][j]=(double)temp_4;<BR> }<BR> /////////////////////<BR> // fread(&temp,datatype,1,in_fp);<BR> // x[k][j]=(double)temp;<BR> // if(x[k][j]<0) x[k][j]+=256;<BR> }<BR> }<BR> // printf("\nWhen Class= %3d,X is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<nl;j++)<BR> {<BR> // printf("%10.5f",x[k][j]);<BR> }<BR> // printf("\n");<BR> }<BR>///////////////////求X的转置矩阵Xt////////////////////////////////<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<nl;j++)<BR> {<BR> //计算训练样本的转置矩阵<BR> xt[j][k]=x[k][j];<BR> }<BR> }</P>9 i3 H- P% n1 Q2 P5 ^+ X# O
> printf("\nWhen Class= %3d,Xt is:\n",i);<BR> for(j=0;j<nl;j++)<BR> {<BR> for(k=0;k<n;k++)<BR> {<BR> // printf("%10.5f",xt[j][k]);<BR> }<BR> // printf("\n");<BR> }</P>
>///////////////////求矩阵u/////////////////////////////////////////<BR> printf("\nWhen Class= %3d,U is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<nl;j++)<BR> { <BR> u[k]+=x[k][j];<BR> }<BR> u[k]=(u[k])/(nl);//求每一类的样本均值<BR> // printf("%13.7lf\n",u[k]);<BR> }</P>
> // printf("\n");</P>
>///////////////////计算矩阵Sigma=X*Xt//////////////////////////////<BR> brmul(xp,xtp,n,nl,n,sigma_tempp);<BR> for (k=0;k<n;k++)<BR> for(j=0;j<n;j++)<BR> sigma_temp[k][j]=sigma_temp[k][j]/(nl-1); //求每一类的散度矩阵<BR> printf("\nWhen Class= %3d,Sigma is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> // printf("%13.3lf",sigma_temp[k][j]);<BR> }<BR> // printf("\n");<BR> }</P>
><BR>///////////////////计算矩阵Sigma的逆矩阵anti_sigma////////////////<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> anti_sigma_temp[k][j]=sigma_temp[k][j];<BR> }<BR> }</P>
> return_brinv=brinv(anti_sigma_tempp,n);//求散度矩阵的逆矩阵</P>3 S& v& v# a7 D1 M
> if(return_brinv==0)//矩阵Sigma奇异,退出程序,重新选择参数<BR> {<BR> free(x);<BR> free(xp);<BR> free(xt);<BR> free(xtp);<BR> free(sigma_tempp);<BR> free(sigma_temp);<BR> free(anti_sigma_tempp);<BR> free(anti_sigma_temp);</P>% B: z) p: H {# \: v9 d
> <BR> <BR> // free(p);<BR> free(nl);<BR> for(k=0;k<m;k++)<BR> {<BR> free(style[k]);<BR> }<BR> free(style);<BR> free(u);<BR> free(up);<BR> free(sigmap);<BR> free(sigmapp);<BR> free(sigma);<BR> free(anti_sigmap);<BR> free(anti_sigmapp);<BR> free(anti_sigma);<BR> return;<BR> }</P>0 ] E1 B/ `. g& O: a' ~6 K
> printf("\nWhen Class= %3d,Anti_Sigma is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> // printf("%13.3e ",anti_sigma_temp[k][j]);<BR> }<BR> // printf("\n");<BR> }</P>/ `$ F6 I) E8 L| 欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) | Powered by Discuz! X2.5 |