//
for (r = 0; r < N-1; r++){ //判断R[i][r](i=r+1,r+2,...,N-1)是否全为零 sum = 0;//变量sum使用前清零 for (i = r+1; i < N; i++){ sum = sum || R[i][r]; } if (sum){ //计算dr sum = 0; for (i = r; i < N; i++){ sum += R[i][r] * R[i][r]; } dr = sqrt(sum);
//计算cr if (R[r][r] > 0) cr = -dr; else cr = dr;
//计算hr
hr = cr * cr - cr * R[r][r]; //取向量ur[]
for (i = 0; i < N; i++){ if (i < r) ur[i] = 0; else if (i == r) ur[i] = R[i][r] - cr; else ur[i] = R[i][r]; }
//计算向量wr[] for (i = 0; i < N; i++){ sum = 0; for (j = r; j < N; j++) sum += Q[i][j] * ur[j]; wr[i] = sum; }
//计算新的矩阵Qr+1[][],存储在Q[][]里面 for (i = 0; i < N; i++) for (j = 0; j < N; j++) Q[i][j] = Q[i][j] - wr[i] * ur[j] / hr; //计算向量pr[] for (i = 0; i < N; i++){ sum = 0; for (j = r; j < N; j++) sum += R[j][i] * ur[j]; pr[i] = sum / hr; }
}
}
} //计算新产生的R[][] for (i = 0; i < N; i++) for (j = 0; j < N; j++) R[i][j] = R[i][j] - ur[i] * pr[j];
/*************************************
Mk[][]矩阵的QR分解及求Ak+1[][]的计算
形参m为每次降阶之后的值 无返回值 **************************************/
void QR_decompositionMk(double mk[][N],int m, double a[][N]){ intr,i,j; double tr,hr,cr,dr,sum; double ur[N],pr[N],qr[N],wr[N],vr[N],br[N][N],Cr[N][N];
//取矩阵br[][]
for (i = 0; i <= m; i++) for (j = 0; j <= m; j++) br[i][j] = mk[i][j]; //取矩阵Cr[][]
for (i = 0; i <= m; i++) for (j = 0; j <= m; j++) Cr[i][j] = a[i][j]; //
for (r = 0; r <= m-1; r++){
//判断br[i][r](i=r+1,r+2,...,m)是否全为零 sum = 0;//变量sum使用前清零 for (i = r+1; i <= m; i++){ sum = sum || br[i][r]; }
if (sum){ //计算dr sum = 0; for (i = r; i <= m; i++){ sum += br[i][r] * br[i][r]; } dr = sqrt(sum); //计算cr if (br[r][r] > 0) cr = -dr; else cr = dr; //计算hr hr = cr * cr - cr * br[r][r]; //取向量ur[] for (i = 0; i <= m; i++){ if (i < r)
}
ur[i] = 0; else if (i == r) ur[i] = br[i][r] - cr; else ur[i] = br[i][r]; }
//计算向量vr[]
for (i = 0; i <= m; i++){ sum = 0; for (j = r; j <= m; j++) sum += br[j][i] * ur[j]; vr[i] = sum / hr; }
//计算新的矩阵Br+1[][],存储在Br[][]里面 for (i = 0; i <= m; i++) for (j = 0; j <= m; j++) br[i][j] = br[i][j] - ur[i] * vr[j]; //计算向量pr[]
for (i = 0; i <= m; i++){ sum = 0; for (j = r; j <= m; j++) sum += Cr[j][i] * ur[j]; pr[i] = sum / hr; }
//计算向量qr[]
for (i = 0; i <= m; i++){ sum = 0; for (j = r; j <= m; j++) sum += Cr[i][j] * ur[j]; qr[i] = sum / hr; }
//计算tr sum = 0;
for (i = r; i <= m; i++) sum += pr[i] * ur[i]; tr = sum / hr;
//计算向量wr[] for (i = 0; i <= m; i++){ if (i < r) wr[i] = qr[i]; else wr[i] = qr[i] - tr * ur[i]; }
//计算新产生的矩阵Cr[][] for (i = 0; i <= m; i++) for (j = 0; j <= m; j++) Cr[i][j] = Cr[i][j] - wr[i] * ur[j] - ur[i] * pr[j];
} //将计算出的Cr[][]矩阵中的值赋给矩阵A[][]，得到新的矩阵。 for (i = 0; i <= m; i++) for (j = 0; j <= m; j++) a[i][j] = Cr[i][j]; } /******************************************* 高斯列主元素消元法求矩阵实特征值的特征向量 等到的特征向量放在X[]里面 返回值为整型 ********************************************/ intGauss_column(double a[][N],double b[],double X[]) { inti,j,k,row_no; double temp,a_ik,m_ik,sum = 0;
//消元过程
for (k = 0; k < N - 1; k++){
//选行号，并记录 row_no = k;
a_ik = abs(a[k][k]); for(i = k + 1; i < N; i++){ if (abs(a[i][k]) >a_ik){ a_ik = abs(a[i][k]); row_no = i; } }
//交换刚刚选择的最大行与第K行的所有元素
if (row_no != k){//如果不是第K行的元素最大，则交换行 for (j = k; j < N; j++){ temp = a[k][j]; a[k][j] = a[row_no][j]; a[row_no][j] = temp; } temp = b[k]; b[k] = b[row_no]; b[row_no] = temp; } if(a[k][k] == 0){ printf(\ return(0); } for(i = k+1; i < N; i++){ m_ik = a[i][k] / a[k][k]; a[i][k] = 0; for(j = k+1; j < N; j++){ a[i][j] = a[i][j] - m_ik * a[k][j]; } b[i] = b[i] - m_ik * b[k];
}
} }
//回代过程，由于a的行列式值为零，所以有无穷多个解。所以不妨设x[N-1] := 1 X[N-1] = 1;//b[N-1] / a[N-1][N-1]; for(k = N-2; k >= 0; k--){ sum = 0; for(j = k+1; j < N; j++) sum += X[j] * a[k][j]; X[k] = (b[k] - sum) / a[k][k]; }
return 0;
/******************************************* 打印特征向量的值到result文件中
无返回值 ********************************************/ void print_lambda(complex lambda[]){ int i; fprintf(pReadFile,\已经求出所有的%d个特征值λ：\\n\ //printf(\已经求出所有的%d个特征值λ：\\n\ for (i = 0; i < N; i++){ if (lambda[i].ipart != 0){ if (lambda[i].ipart> 0) fprintf(pReadFile,\ //printf(\ else fprintf(pReadFile,\ //printf(\ } else fprintf(pReadFile,\ //printf(\ } fprintf(pReadFile,\ //printf(\} /******************************************* 打印参量传入矩阵到result文件中 无返回值 ********************************************/ void print_matrix(double mat[][N]){ inti,j; for (i = 0; i < N; i++){ for (j = 0; j < N; j++)