?? qrmethod.cpp
字號:
#include "math.h"
#include "iostream.h"
#include "stdio.h"
#include <iomanip.h>
#define BOOL int
#define FALSE 0
#define TRUE 1
double ee = 0.00000000001; //精度水平ee=1e-12;
//顯示矩陣
void ShowMatrix(double * P_Matrix,int N)
{
cout<<endl;
for(int i = 0; i < N; i++)
{
cout<<endl<<i+1<<" row numbers are:"<<endl;
for(int j = 0; j < N; j++)
{
cout<<setiosflags(ios::scientific)<<setprecision(11)<<P_Matrix[i*N+j]<<" ";
}
}
}
//顯示列向量
void ShowVector(double * P_Vector,int N)
{
cout<<endl;
for(int i = 0; i < N; i++)
{
cout<<setiosflags(ios::scientific)<<setprecision(11)<<P_Vector[i]<<endl;
}
}
//顯示所有特征值Lmuda
void ShowVectorLumda(double * pf_LmudaRe,double * pf_LmudaIm,int N)
{
cout<<endl;
for(int i = 0; i < N; i++)
{
if(pf_LmudaIm[i] == 0)
{
cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[i]<<endl;
}
if(pf_LmudaIm[i] < 0)
{
cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[i]<<" - "<<-pf_LmudaIm[i]<<" i"<<endl;
}
if(pf_LmudaIm[i] > 0)
{
cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[i]<<" + "<<pf_LmudaIm[i]<<" i"<<endl;
}
}
}
//矩陣乘法——相同維數矩陣相乘
double * Multiply(double * pf_MatraixA, double * pf_MatraixB, int N)
{
int i,j,k;
double * pf_Result = new double[N*N];
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
pf_Result[i*N+j] = 0;
for(k = 0; k < N; k++)
{
pf_Result[i*N+j] += pf_MatraixA[i*N+k]*pf_MatraixB[k*N+j];
}
}
}
return pf_Result;
}
//矩陣乘法——矩陣乘以列向量
double * MultiplyMV(double * pf_MatraixA, double * pf_Vector, int N)
{
int i,k;
double * pf_Result = new double[N];
for(i = 0; i < N; i++)
{
pf_Result[i] = 0;
for(k = 0; k < N; k++)
{
pf_Result[i] += pf_MatraixA[i*N+k]*pf_Vector[k];
}
}
return pf_Result;
}
//矩陣乘法——數乘以矩陣
double * MultiplyNumM(double * pf_MatraixA, double num, int N)
{
int i,j;
double * pf_Result = new double[N*N];
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
pf_Result[i*N+j] = pf_MatraixA[i*N+j] * num;
}
}
return pf_Result;
}
//行向量乘以列向量
double MultiplyRC(double * pf_VectorA, double * pf_VectorB, int N)
{
int i;
double temp = 0;
for(i = 0; i < N; i++)
{
temp += pf_VectorA[i]*pf_VectorB[i];
}
return temp;
}
//列向量乘以行向量
double * MultiplyCR(double * pf_VectorA, double * pf_VectorB, int N)
{
int i,j;
double temp = 0;
double * pf_Result = new double[N*N];
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
pf_Result[i*N+j] = pf_VectorA[i]*pf_VectorB[j];
}
}
return pf_Result;
}
//矩陣轉置
double * MatrixT(double * pf_Matraix, int N)
{
int i,j;
double * pf_MatraixT = new double[N*N];
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
pf_MatraixT[i*N+j] = pf_Matraix[j*N+i];
}
}
return pf_MatraixT;
}
//向量乘以數
double * VectorMultiNum(double * pf_Vector, double num, int N)
{
int i;
double * pf_Result = new double[N];
for(i = 0; i < N; i++)
{
pf_Result[i] = pf_Vector[i]*num;
}
return pf_Result;
}
//向量相減
double * VectorSub(double * pf_VectorA, double * pf_VectorB, int N)
{
int i;
double temp;
double * pf_Result = new double[N];
for(i = 0; i < N; i++)
{
pf_Result[i] = pf_VectorA[i] - pf_VectorB[i];
}
return pf_Result;
}
//矩陣相減
BOOL MatrixSub(double * pf_MatrixA, double * pf_MatrixB, double * pf_Result, int N)
{
int i,j;
for(i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
pf_Result[i*N+j] = pf_MatrixA[i*N+j] - pf_MatrixB[i*N+j];
}
}
return TRUE;
}
//對矩陣第r列進行上三角變換
BOOL HessOneStep(double * pf_Matraix,int r, int N)
{
int i;
double dr,cr,hr,temp;
temp = 0;
for(i = r+1; i < N; i++)
{
temp += pf_Matraix[i*N+r]*pf_Matraix[i*N+r];
}
dr = sqrt(temp);
if (pf_Matraix[(r+1)*N+r] <= 0)
{
cr = dr;
}
else
{
cr = -dr;
}
hr = cr*cr - cr*pf_Matraix[(r+1)*N+r];
double * ur = new double [N];
double * pr;
double * qr;
double tr = 0;
double * wr;
for(i = 0; i < N; i++)
{
if (i <= r)
{
ur[i] = 0;
}
if (i == r+1)
{
ur[i] = pf_Matraix[i*N+r] - cr;
}
if( i > r+1)
{
ur[i] = pf_Matraix[i*N+r];
}
}
double * pf_MatraixT;
pf_MatraixT = MatrixT(pf_Matraix,N);
pr = VectorMultiNum(MultiplyMV(pf_MatraixT,ur,N),1/hr,N);
qr = VectorMultiNum(MultiplyMV(pf_Matraix,ur,N),1/hr,N);
tr = MultiplyRC(pr,ur,N)/hr;
wr = VectorSub(qr,VectorMultiNum(ur,tr,N),N);
MatrixSub(pf_Matraix,MultiplyCR(wr,ur,N),pf_Matraix,N);
MatrixSub(pf_Matraix,MultiplyCR(ur,pr,N),pf_Matraix,N);
delete ur;
delete pr;
delete qr;
delete wr;
delete pf_MatraixT;
return TRUE;
}
//矩陣擬上三角化
double * Hessenberg(double * pf_Matraix,int N)
{
double * pf_HessMatraix = new double[N*N];
int i,j;
//矩陣：A1 = A
for(i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
pf_HessMatraix[i*N+j] = pf_Matraix[i*N+j];
}
}
int r;
BOOL flag;
for(r = 0; r < N-2; r++)
{
flag = FALSE;
for(i = r+2; i < N; i++)
{
if (fabs(pf_HessMatraix[i*N+r]) > ee)
{
flag = TRUE;
break;
}
}
HessOneStep(pf_HessMatraix,r,N);
}
for(i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
if (fabs(pf_HessMatraix[i*N+j]) <= ee)
{
pf_HessMatraix[i*N+j] = 0;
}
}
}
return pf_HessMatraix;
}
void QRDDStep4(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count);
void QRDDStep567(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count);
void QRDDStep8910(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count);
void QRDDStep11(int& m);
void QROneStep(double * pf_Matraix,double * pf_MatraixM,int r,int N);
//step(3)
void QRDDStep3(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
if (fabs(pf_Matraix[(m-1)*N+m-2]) <= ee)
{
pf_LmudaRe[count] = pf_Matraix[(m-1)*N+m-1];
count++;
m = m-1;
QRDDStep4(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
}
else
{
QRDDStep567(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
}
}
//step(4)
void QRDDStep4(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
if (m == 1)
{
pf_LmudaRe[count] = pf_Matraix[0];
count++;
QRDDStep11(m);
}
if(m == 0)
{
QRDDStep11(m);
}
if (m > 1)
{
QRDDStep3(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
}
}
void GetTwoLumda(double * pf_LmudaRe,double * pf_LmudaIm,double b,double c,int& count)
{
double Delta; //delta = b*b - 4ac
Delta = b*b-4*c;
if ( Delta >= 0)
{
pf_LmudaRe[count] = (-b + sqrt(Delta))/2;
count++;
pf_LmudaRe[count] = (-b - sqrt(Delta))/2;
count++;
}
else
{
pf_LmudaRe[count] = -b/2;
pf_LmudaIm[count] = sqrt(-Delta)/2;
count++;
pf_LmudaRe[count] = -b/2;
pf_LmudaIm[count] = -sqrt(-Delta)/2;
count++;
}
}
//step(5)&&step(6)&&step(7)
void QRDDStep567(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
double s1,s2;
int i;
double * pf_MatraixD;
pf_MatraixD = new double[4];
pf_MatraixD[0] = pf_Matraix[(m-2)*N+m-2];
pf_MatraixD[1] = pf_Matraix[(m-2)*N+m-1];
pf_MatraixD[2] = pf_Matraix[(m-1)*N+m-2];
pf_MatraixD[3] = pf_Matraix[(m-1)*N+m-1];
double b,c; //二階方程的系數b,c
b = -(pf_MatraixD[0]+pf_MatraixD[3]);
c = ((pf_MatraixD[0] * pf_MatraixD[3]) - (pf_MatraixD[1] * pf_MatraixD[2]));
if (m == 2)
{
GetTwoLumda(pf_LmudaRe,pf_LmudaIm,b,c,count);
QRDDStep11(m);
}
else
{
if (fabs(pf_Matraix[(m-2)*N+m-3]) <= ee)
{
GetTwoLumda(pf_LmudaRe,pf_LmudaIm,b,c,count);
m = m - 2;
QRDDStep4(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
}
else
{
QRDDStep8910(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
}
}
}
//step(8)&&step(9)&&step(10)
void QRDDStep8910(double * pf_Matraix, double * pf_LmudaRe,double * pf_LmudaIm, int N, int k,int& m,int& count)
{
int i,j,r;
double * pf_MatraixA = new double[m*m];
double * pf_MatraixM = new double[m*m];
for(i = 0; i < m; i++)
{
for(j = 0; j < m; j++)
{
pf_MatraixA[i*m+j] = pf_Matraix[i*N+j];
}
}
double s,t;
s = pf_Matraix[(m-2)*N+m-2] + pf_Matraix[(m-1)*N+m-1];
t = (pf_Matraix[(m-2)*N+m-2] * pf_Matraix[(m-1)*N+m-1])
- (pf_Matraix[(m-1)*N+m-2] * pf_Matraix[(m-2)*N+m-1]);
double * pf_MatraixI= new double[m*m]; //單位矩陣
//單位矩陣賦值
for(i = 0; i < m; i++)
{
for(j = 0; j < m; j++)
{
if (i == j)
{
pf_MatraixI[i*m+j] = 1;
}
else
{
pf_MatraixI[i*m+j] = 0;
}
}
}
MatrixSub(Multiply(pf_MatraixA,pf_MatraixA,m),MultiplyNumM(pf_MatraixA,s,m),pf_MatraixM,m); //M1 = A^2-sA
MatrixSub(pf_MatraixM,MultiplyNumM(pf_MatraixI,-t,m),pf_MatraixM,m); //M = M1-(-tI)
BOOL flag;
for(r = 0; r < m - 1; r++)
{
flag = FALSE;
for(i = r+1; i < m; i++)
{
if (fabs(pf_MatraixM[i*m+r]) > ee)
{
flag = TRUE;
break;
}
}
QROneStep(pf_MatraixM,pf_MatraixA,r,m);
}
//矩陣A的前m*m元素賦值為QR分解后的矩陣
for(i = 0; i < m; i++)
{
for(j = 0; j < m; j++)
{
pf_Matraix[i*N+j] = pf_MatraixA[i*m+j];
}
}
delete pf_MatraixM;
delete pf_MatraixI;
delete pf_MatraixA;
}
//step(11)
void QRDDStep11(int& m)
{
m = 0;
}
//用具體簡化算法進行QR分解
void QROneStep(double * pf_MatraixB,double * pf_MatraixC,int r,int m)
{
int i;
double dr,cr,hr,temp;
temp = 0;
for(i = r; i < m; i++)
{
temp += pf_MatraixB[i*m+r]*pf_MatraixB[i*m+r];
}
dr = sqrt(temp);
if (pf_MatraixB[r*m+r] <= 0)
{
cr = dr;
}
else
{
cr = -dr;
}
hr = cr*cr - cr*pf_MatraixB[r*m+r];
double * ur = new double [m];
double * pr;
double * qr;
double * vr;
double tr = 0;
double * wr;
for(i = 0; i < m; i++)
{
if (i < r)
{
ur[i] = 0;
}
if (i == r)
{
ur[i] = pf_MatraixB[i*m+r] - cr;
}
if( i > r)
{
ur[i] = pf_MatraixB[i*m+r];
}
}
vr = VectorMultiNum(MultiplyMV(MatrixT(pf_MatraixB,m),ur,m),1/hr,m);
MatrixSub(pf_MatraixB,MultiplyCR(ur,vr,m),pf_MatraixB,m); //求B(r+1)
pr = VectorMultiNum(MultiplyMV(MatrixT(pf_MatraixC,m),ur,m),1/hr,m);
qr = VectorMultiNum(MultiplyMV(pf_MatraixC,ur,m),1/hr,m);
tr = MultiplyRC(pr,ur,m) / hr;
wr = VectorSub(qr,VectorMultiNum(ur,tr,m),m);
MatrixSub(pf_MatraixC,MultiplyCR(wr,ur,m),pf_MatraixC,m); //求C(r+1)
MatrixSub(pf_MatraixC,MultiplyCR(ur,pr,m),pf_MatraixC,m); //求C(r+1)
delete ur;
delete pr;
delete qr;
delete wr;
delete vr;
}
//對擬上三角化的矩陣進行QR分解求特征值
BOOL QRDoubleDispace(double * pf_Matraix,int N, double * pf_LmudaRe,double * pf_LmudaIm)
{
int i,k,m,count;
count = 0;
m = N;
while (m != 0)
{
QRDDStep3(pf_Matraix,pf_LmudaRe,pf_LmudaIm,N,k,m,count);
}
for(i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
if (fabs(pf_Matraix[i*N+j]) <= ee)
{
pf_Matraix[i*N+j] = 0;
}
}
}
return TRUE;
}
//對矩陣進行選主元素高斯消去
BOOL GaussRe(double * P_Matrix,double * P_MatrixG, int N)
{
double m,temp;
int i,j,k,l;
int flag = 0;
if (P_Matrix[0] == 0)
{
return FALSE;
}
//復制原始矩陣和向量
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
P_MatrixG[i*N+j] = P_Matrix[i*N+j];
}
}
//開始消去法
for(k = 0; k < N-1; k++)
{
temp = P_MatrixG[k];
flag = k;
//找出主元素
for(l = k; l < N; l++ )
{
if (fabs(P_MatrixG[l*N+k]) > fabs(temp))
{
temp = P_MatrixG[l*N+k];
flag = l;
}
}
//交換行
if (flag != k)
{
for(l = k; l < N; l++)//交換系數矩陣元素
{
temp = P_MatrixG[flag*N+l];
P_MatrixG[flag*N+l] = P_MatrixG[k*N+l];
P_MatrixG[k*N+l] = temp;
}
}
for(i = k+1; i < N; i++)
{
m = P_MatrixG[i*N+k]/P_MatrixG[k*N+k];
P_MatrixG[i*N+k] = 0; //消去的元素為0
for(j = k+1; j < N; j++)
{
P_MatrixG[i*N+j] = P_MatrixG[i*N+j]-m*P_MatrixG[k*N+j];
}
}
}
return TRUE;
}
//回代求解
BOOL Calculate(double * P_MatrixG, double * P_VectorX, int N)
{
int j,k;
double temp;
P_VectorX[N-1] = 1; //因為是齊次線性方程組，將最后一個x直接賦值為1
for(k = N-2; k >= 0; k--)
{
temp = 0;
for(j = k+1; j < N; j++)
{
temp += P_MatrixG[k*N+j]*P_VectorX[j];
}
P_VectorX[k] = -temp/P_MatrixG[k*N+k];
}
return TRUE;
}
//求出每個實特征值的特征向量
void CalVectors(double * pf_Matrix,int N, double * pf_LmudaRe,double * pf_LmudaIm)
{
int i,j,k;
double * pf_MatrixASubLmuda = new double[N*N];
double * pf_VectorX = new double[N];
for (k = 0; k < N; k++)
{
if (pf_LmudaIm[k] != 0) //若特征值不是實數，不進行計算
{
continue;
}
cout<<endl<<"This is the vector of "<<k+1<<" (";
cout<<setiosflags(ios::scientific)<<setprecision(11)<<pf_LmudaRe[k]<<"):";
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
if (i == j)
{
pf_MatrixASubLmuda[i*N+j] = pf_Matrix[i*N+j] - pf_LmudaRe[k];
}
else
{
pf_MatrixASubLmuda[i*N+j] = pf_Matrix[i*N+j];
}
}
}
double * pf_MatrixG = new double[N*N];
GaussRe(pf_MatrixASubLmuda,pf_MatrixG,N);
Calculate(pf_MatrixG,pf_VectorX,N);
ShowVector(pf_VectorX,N);
}
}
void main()
{
int N;
int i,j;
N = 10;
double * pf_MatrixA = new double[N*N];
double * pf_Result;
double * pf_LmudaIm = new double[N];
double * pf_LmudaRe = new double[N];
//Lumda虛部為0
for(i = 0; i < N; i++)
{
pf_LmudaIm[i] = 0;
}
for(i = 0; i < N; i++)
{
for(j = 0; j < N; j++)
{
if (i == j)
{
pf_MatrixA[i*N+j] = 1.5*cos(i+1+1.2*(j+1));
}
else
{
pf_MatrixA[i*N+j] = sin(0.5*(i+1)+0.2*(j+1));
}
}
}
pf_Result = Hessenberg(pf_MatrixA,N);
cout<<"The hessenberg matrix is:"<<endl;
ShowMatrix(pf_Result,N);
QRDoubleDispace(pf_Result,N,pf_LmudaRe,pf_LmudaIm);
cout<<endl<<"The matrix after QR is :";
ShowMatrix(pf_Result,N);
cout<<endl<<"The Lmudas are :";
ShowVectorLumda(pf_LmudaRe,pf_LmudaIm,N);
cout<<endl<<"The X vectors are :";
CalVectors(pf_MatrixA,N,pf_LmudaRe,pf_LmudaIm);
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -