Question

谁有编好的求矩阵特征值的c++程序？

最好是用并行计算采用OpenMP的

Answer 1

#include<iostream.h>
#include<stdlib.h>
#include<math.h>
#include<iomanip.h>

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//class Matrix定义矩阵类

const int Max_xy = 20; //矩阵的最大维数
class Matrix
{
private:
double data[Max_xy][Max_xy];
unsigned x, y; //x,y;
public:
Matrix(); //默认构造函数
Matrix(const Matrix &source); //拷贝构造函数
void creat(); //输入矩阵
void init();
void transpose(); //矩阵转置
void show(); //输入此矩阵
double mode() const; //求一维矩阵的长度
void check_shiduichen(); //检查是否为是对称矩阵
void creat_unit(unsigned i); //生成i行单位矩阵
void set_x(unsigned xx); //设置行数
void set_y(unsigned yy); //设置列数
unsigned get_x(); //得到行数
unsigned get_y(); //得到列数
void shucheng(double changshu); //数乘运算
void setdata(unsigned i, unsigned j, double source); //定位输入数据
double getdata(unsigned i, unsigned j); //定位得到数据
void sturm(); //求特征值
unsigned sturm_s(double m); //计算sturm系列的同好数
Matrix operator = (const Matrix &right);
friend Matrix &operator + (const Matrix &left, const Matrix &right); //重载+号
friend Matrix &operator - (const Matrix &left, const Matrix &right); //重载-号
friend Matrix &operator * (const Matrix &left, const Matrix &right); //重载乘号
friend ostream &operator <<(ostream &os, const Matrix &source); //重载输出
friend void Householder(Matrix &source); //用Householde矩阵将实对称矩阵化为三对角矩阵

};

Matrix temp_Matrix; //全局变量Matrix

//===================================================================

//--------------------默认构造函数
Matrix::Matrix()
{
init();
}

//----------------------------拷贝构造函数
Matrix::Matrix(const Matrix &source)
{
init();
x = source.x;
y = source.y;
for(unsigned i = 0; i < x; i++)
for(unsigned j = 0; j < y; j++)
data[i][j] = source.data[i][j];
}

//------------------------------------------初始化矩阵元素
void Matrix::init()
{
x = y = 0;
for(unsigned i = 0; i < Max_xy; i++)
for(unsigned j = 0; j < Max_xy; j++)
data[i][j] = 0;

}

//------------------------------矩阵转置
void Matrix::transpose()
{
double temp;
int t;
for(unsigned i = 0; i < Max_xy; i++)
for(unsigned j = 0; j <= i; j++)
{
temp = data[i][j];
data[i][j] = data[j][i];
data[j][i] = temp;
}
t = x;
x = y;
y = t;
}

//--------------------------------------求一维矩阵的长度
double Matrix::mode() const
{
double s = 0;
unsigned i, j;
if(x == 1)
for(i = 0, j = 0; j < y; j++)
s += data[i][j] * data[i][j];
else if(y == 1)
for(i = 0, j = 0; i < x; i++)
s += data[i][j] * data[i][j];
else
{
cout << "\n不是一维的!";
exit(0);
}
s = sqrt(s);
return (s);
}

//----------------------------------------重载=号
Matrix Matrix::operator = (const Matrix &source)
{
x = source.x;
y = source.y;
for(unsigned i = 0; i < x; i++)
for(unsigned j = 0; j < y; j++)
data[i][j] = source.data[i][j];
return *this;

}

//-------------------------------------------重载+号
Matrix &operator + (const Matrix &left, const Matrix &right)
{
if(left.x != right.x || left.y != right.y)
{
cout << "\n维数不相等,不能相加!";
exit(0);
}
for(unsigned i = 0; i < left.x; i++)
for(unsigned j = 0; j < left.y; j++)
temp_Matrix.data[i][j] = right.data[i][j] + left.data[i][j];
temp_Matrix.x = right.x;
temp_Matrix.y = right.y;
return temp_Matrix;
}

//---------------------------------------------重载+号
Matrix &operator - (const Matrix &left, const Matrix &right)
{
if(left.x != right.x || left.y != right.y)
{
cout << "\n维数不相等,不能相减!";
exit(0);
}
for(unsigned i = 0; i < left.x; i++)
for(unsigned j = 0; j < left.y; j++)
temp_Matrix.data[i][j] = left.data[i][j] - right.data[i][j];
temp_Matrix.x = right.x;
temp_Matrix.y = right.y;
return temp_Matrix;
}

//----------------------------------------重载乘号
Matrix &operator * (const Matrix &left, const Matrix &right)
{
if(left.y != right.x)
{
cout << "\n两个矩阵相乘错误.";
exit(0);
}
temp_Matrix.init();
unsigned i, j, k;
for(i = 0; i < left.x; i++)
for(j = 0; j < right.y; j++)
for(k = 0; k < left.y; k++)
temp_Matrix.data[i][j] += left.data[i][k] * right.data[k][j];

temp_Matrix.x = left.x;
temp_Matrix.y = right.y;
return temp_Matrix;
}

//-------------------------------------------输入矩阵
void Matrix::creat()
{
cout << "输如行列式:";
cout << "\n行数:";
cin >> x;
cout << "列数:";
cin >> y;
for(unsigned i = 0; i < x; i++)
{
cout << "输入第" << i + 1 << "行:";
for(unsigned j = 0; j < y; j++)
cin >> data[i][j];
}

}

//----------------------------------------------输出矩阵
void Matrix::show()
{
unsigned i, j;
cout << "\n\n矩阵表示如下:";
for(i = 0; i < x; i++)
{
cout << endl;
for(j = 0; j < y; j++)
cout << setw(7) << setiosflags(ios::left) << data[i][j];
}
}

//----------------------------------用Householder矩阵化为实对称矩阵
void Householder(Matrix &source)
{
unsigned i, lenth, k, m, n, flag;
double temp_lie_x[Max_xy], temp_lie_y[Max_xy], temp[Max_xy];
double s;
for(i = 0; i < source.x - 2; i++)
{
for(k = 0; k < Max_xy; k++)
{
//初始化为0
temp_lie_x[k] = 0;
temp_lie_y[k] = 0;
temp[k] = 0;
}

for(lenth = 0; lenth + i + 1 < source.x; lenth++) //提取第data[i+1][i]到data[x-1][i]的数据存到temp_lie[Max_xy];
temp_lie_x[lenth] = source.data[lenth+i+1][i];
for(k = 0, s = 0; k < lenth; k++) //j为临时变量的个数.
s += temp_lie_x[k] * temp_lie_x[k];
s = sqrt(s);
temp_lie_y[0] = -s;
for(k = 1; k < lenth; k++)
temp_lie_y[k] = 0;

for(k = 0; k < lenth; k++)
temp[k] = temp_lie_x[k] - temp_lie_y[k];

for(k = 0, flag = 0; k < lenth; k++) //假如以上两个向量相等则退出,则跳出以进行下一次变换
if(temp[k] != 0)
{
flag = 1;
break;
}
if(flag == 0)
continue;

Matrix part_h, I, x_y; //定义Matrix变量

I.creat_unit(lenth);
x_y.x = lenth, x_y.y = 1;
//对x_y赋值
for(k = 0; k < lenth; k++)
x_y.data[k][0] = temp[k];
//求x_y的转置
Matrix zhuanzhi_x_y( x_y ); //拷贝构造函数
zhuanzhi_x_y.transpose();
s = 2.0 / ( x_y.mode() * x_y.mode() );
x_y.shucheng(s);
temp_Matrix = x_y * zhuanzhi_x_y;
part_h = I - x_y * zhuanzhi_x_y;
Matrix H;
H.creat_unit(source.x);
for(m = i + 1; m < source.x; m++)
for(n = i + 1; n < source.y; n++)
H.data[m][n] = part_h.data[m-i-1][n-i-1]; //得到最后的Householder矩阵
source = source * H;
source = H * source;

}

for(i = 0; i < source.x; i++)
for(k = 0; k < source.y; k++)
if(fabs(source.data[i][k]) < 1e-13)
source.data[i][k] = 0;

}

//------------------------------------检查是否为实对称矩阵
void Matrix::check_shiduichen()
{
if(x != y)
{
cout << "\n\n不是是对称矩阵(行列不相等)\n\n";
exit(0);
}

for(unsigned i = 0; i < x; i++)
for(unsigned j = 0; j < y; j++)
{
if(data[i][j] != data[j][i])
{
cout << "\n\n不是实对称矩阵!(不对称!)\n\n";
exit(0);
}
}
}

ostream &operator <<(ostream &os, const Matrix &source)
{
unsigned i, j;
for(i = 0; i < source.x; i++)
{
os << "\n";
for(j = 0; j < source.y; j++)
os << setw(10) << setiosflags(ios::left) << source.data[i][j] << "\t";
}
os << endl;
return os;
}

void Matrix::set_x(unsigned xx) //设置行数
{
x = xx;
}

void Matrix::set_y(unsigned yy) //设置列数
{
y = yy;
}

unsigned Matrix::get_x() //得到行数
{
return x;
}

unsigned Matrix::get_y() //得到列数
{
return y;
}

double Matrix::getdata(unsigned i, unsigned j) //定位得到数据
{
return data[i][j];
}

void Matrix::setdata(unsigned i, unsigned j, double source)
{
data[i][j] = source;
if(x < i) x = i;
if(y < j) y = j;
}

void Matrix::creat_unit(unsigned lenth) //生成i行单位矩阵
{
init();
x = y = lenth;
for(unsigned i = 0; i < lenth; i++)
data[i][i] = 1;
}

void Matrix::shucheng(double changshu) //数乘运算
{
unsigned i, j;
for(i = 0; i < x; i++)
for(j = 0; j < y; j++)
data[i][j] *= changshu;
}

void Matrix::sturm()
{
double s, r, maxhang, a, b;
unsigned i, j, m;
for(i = 0, maxhang = 0; i < x; i++)
{
for( s = 0, j = 0; j < y; j++)
s += fabs(data[i][j]);
if(s > maxhang)
maxhang = s;
}
a = -maxhang;
b = maxhang;
m = sturm_s(a) - sturm_s(b);
for(i = 1; i <= m; i++)
{

a = -maxhang;
b = maxhang;
do
{
r = 0.5 * (a + b);
if(sturm_s(r) >= i)
a = r;
else
b = r;
}
while(fabs(a - b) > 1e-11);
cout << "\n特征值" << i << ": " << setiosflags(ios::left) << setw(10) << setprecision(10)
<< setiosflags(ios::fixed) << 0.5 * (a + b);
}
}

unsigned Matrix::sturm_s(double r)
{
double p[Max_xy+1];
int temp[Max_xy+1];
unsigned m, i;
p[0] = 1;
p[1] = data[0][0] - r;
for(i = 2; i <= x; i++)
p[i] = (data[i-1][i-1] - r) * p[i-1] - data[i-2][i-1] * data[i-2][i-1] * p[i-2];
temp[0] = 1;
for(i = 1; i <= x; i++)
if(p[i] > 1e-14)
temp[i] = 1;
else if(p[i] < -1e-14)
temp[i] = -1;
else
temp[i] = temp[i-1];

for(i = 1, m = 0; i <= x; i++)
if(temp[i] + temp[i-1] != 0)
m++;
return m;

}

void main()
{
Matrix a;
a.creat();
cout << "输入的矩阵为:";
cout << a;

a.sturm();
cout << endl;
}

Answer 2

这个矩阵输入上有问题