/*求beta2*/ beta2=0; for(i=1;i<=n;i++) { beta2+=y[i]*u[i]; } }while(fabs((beta2-beta1)/beta2) > 1e-12 || k==1); free(u); return beta2; }
/*------------------------------------------------------ ***函数名：fanmifa(double **a,int r,int s,int n,double *y) ***功 能：计算矩阵A按模最大的特征值和对应的特征向量 ***参 数：矩阵A,y为所求特征值对应的特征向量
***返回值：double型特征值
***说 明：~~针对带状存储的矩阵~~ ~~使用的是2范数~~ by lidezheng 日期：2013-10-18
------------------------------------------------------*/ double fanmifa(double **a,int r,int s,int n,double *y) { double beta1,beta2,max; int k,i; double *u; u=(double *)malloc((n+1)*sizeof(double)); /*选取非零向量u*/ for(i=1;i<=n;i++) { u[i]=1; /***--->###这里要讨论初始向量对解的影响###<---***/ } /*反幂法迭代*/ k=0; beta2=0; do{ k++; beta1=beta2; max=0; for(i=1;i<=n;i++)/*求出向量u的2范数*/ { max+=pow((*(u+i)),2); } max=sqrt(max); /*求y(k-1)*/ for(i=1;i<=n;i++) { y[i]=u[i]/max; } /*求Uk,调用LU方法*/ if(LU(a,y,r,s,n,u) == 0)/****注意a,y,u的下标******/ { printf(\调用LU计算方程组出错，ERROR!\\n\); exit(0); } /*求beta2*/ beta2=0;
for(i=1;i<=n;i++) { beta2+=y[i]*u[i]; } }while(fabs((beta2-beta1)/beta2) > 1e-12 || k==1); free(u); return 1/beta2; }
/*------------------------------------------------------
***函数名：LU(double **A,double *B,int r,int s,int n,double *x) ***功 能：Doolittle分解法，解线性方程组 ***参 数：~~解带状存储的矩阵~~
***返回值：0-失败，1-成功 ***说 明：~~针对带状矩阵~~ by lidezheng 日期：2013-10-18 修改：2013-10-19
------------------------------------------------------*/ int LU(double **A,double *B,int r,int s,int n,double *x) { int k,i,j,t,up,down; double temp,**c,*b; c=dmatrix(r+s+1,n); b=(double *)malloc((n+1)*sizeof(double)); /*给局部变量c和b赋初值*/ for(i=1;i<=(r+s+1);i++) for(j=1;j<=n;j++) c[i][j]=A[i][j]; for(i=1;i<=n;i++) *(b+i)=*(B+i); /*lu分解过程*/ for(k=1;k<=n;k++) { if(c[s+1][k] == 0) /*判断Ukk是否为0*/ { return 0; } up=(k+r)>n?n:(k+r); for(j=k;j<=up;j++)/*求U*/ { temp=0; down=((k-r)>(j-s))?(k-r):(j-s); if(down < 1) down=1; for(t=down;t<=k-1;t++) temp+=c[k-t+s+1][t]*c[t-j+s+1][j]; c[k-j+s+1][j]-=temp; } up=((k+r)>n)?n:(k+r); for(i=k+1;i<=up && k
for(t=down;t<=k-1;t++) temp+=c[i-t+s+1][t]*c[t-k+s+1][k]; c[i-k+s+1][k]-=temp; c[i-k+s+1][k]/=c[s+1][k]; } } /*求Ly=b中的y,并将y直接存储到b中*/ for(i=2;i<=n;i++) { temp=0; down=i-r; if(down < 1) down=1; for(t=down;t<=i-1;t++) temp+=c[i-t+s+1][t]*b[t]; b[t]-=temp; } /*求Ux=y中的x*/ x[n]=b[n]/c[s+1][n]; for(i=n-1;i>=1;i--) { temp=0; up=((i+s)>n)?n:(i+s); for(t=i+1;t<=up;t++) temp+=c[i-t+s+1][t]*x[t]; x[i]=(b[i]-temp)/c[s+1][i]; } free_dmatrix(c); free(b); return 1; }
/*------------------------------------------------------ ***函数名：det(double **A,int r,int s,int n)
***功 能：采用Doolittle分解法，用U求矩阵的行列式的值 ***参 数：~~带状存储的矩阵~~
***返回值：0-失败，1-成功 ***说 明：~~针对带状矩阵~~ by lidezheng 日期：2013-10-19 修改：2013-10-19
------------------------------------------------------*/ int det(double **A,int r,int s,int n,double *pdet) { int k,i,j,t,up,down; double temp,**c; c=dmatrix(r+s+1,n); /*给局部变量c和b赋初值*/ for(i=1;i<=(r+s+1);i++) for(j=1;j<=n;j++) c[i][j]=A[i][j]; /*lu分解过程*/ for(k=1;k<=n;k++) { if(c[s+1][k] == 0) /*判断Ukk是否为0*/ { return 0;
} up=(k+r)>n?n:(k+r); for(j=k;j<=up;j++)/*求U*/ { temp=0; down=((k-r)>(j-s))?(k-r):(j-s); if(down < 1) down=1; for(t=down;t<=k-1;t++) temp+=c[k-t+s+1][t]*c[t-j+s+1][j]; c[k-j+s+1][j]-=temp; } up=((k+r)>n)?n:(k+r); for(i=k+1;i<=up && k
/*------------------------------------------------------ ***函数名：dmatrix(int m,int n)
***功 能：生成坐标从1开始的m*n的矩阵 ***参 数：m-行数，n-列数
***返回值：矩阵的首地址
***说 明：要与free_dmatrix(double *a)配套使用 free_dmatrix(double *a)为释放分配的空间的函数 by lidezheng 日期：2013-10-20
------------------------------------------------------*/ double **dmatrix(int m,int n) { int i; double **a; a=(double **)malloc((m+1)*sizeof(double)); if(!a) printf(\); a[1]=(double *)malloc((m*n+1)*sizeof(double)); if(!a[1]) printf(\); for(i=2;i<=m;i++) a[i]=a[i-1]+n; return a; }
void free_dmatrix(double **a) {
free(a[1]); free(a); }
/******完******/