哪位大侠写过关于蚂蚁算法的程序啊

[_xiaolifeidao]

救命啊：）高分求教！

[milpas]

什么是蚂蚁算法???

[_xiaolifeidao]

哎呀，这个说来话长啦
我也不怎么清楚的
我导师让我写程序实现这个算法
天那！

[Arter]

蚁群/蚂蚁算法好多论文上有,代码自己写,不要怕,开始都是很难得!

[wanbaocheng]

我这里有蚂蚁算法的C++源码，是求解旅行售货员（TSP)问题的。不知怎么样？

[wanbaocheng]

把以下6个文件存在同一目录下，然后对anttsp.cpp编译即可。

文件1：anttsp.h

算法源码如下：
anttsp.h
/* C++源程序
* 本文件是求解旅行售货员（TSP）问题的蚂蚁算法的头文件
* 算法的几个主要参数：
* arfa——轨迹的相对重要性（arfa>=0）
* beta——能见度的相对重要性（beta>=0）
* rou ——轨迹的持久性（0<=rou<1），可将1-rou理解为迹衰减度
*/
#ifndef AntTspH
#define AntTspH
class AntTsp{
private:
//double *p;
//int *allow;
matrix d;
int n;
int randchoice(double *pp,int lengthp);
void antk_path(double *t,double *y,double r,double b,double Q,
int *path1,double *z,double *dt);
public:
double arfa,beta,c,rou;
int m,nc;
AntTsp(void){arfa=1;beta=1;c=0.01;rou=0.8;m=0;nc=0;};
~AntTsp(){};
double Ant(int *pathmin);
};
#endif //AntTspH

文件2：anttsp.cpp

#include<iostream.h>
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include"anttsp.h"
#include"matrix.cpp"
double AntTsp::rand0(void){
return (rand()%1000)/1000.0;
}

int AntTsp::randchoice(double *p,int lengthp)
{

double *t,b;
int high=lengthp,low=0,mid,i;
t=(double *)malloc((lengthp+1)*sizeof(double));
if(t==NULL) { printf("assign fail!\n"); exit(0); }
t[0]=0.0;
for(i=1;i<=lengthp;i++)
t[i]=p[i-1];
for(i=2;i<=lengthp;i++)
t[i]+=t[i-1];
b=rand0()*t[lengthp];
while(high-low>=2)
{
mid=floor((high+low)/2.0);
if(b>t[mid]) low=mid;
else high=mid;
}
free(t);
return(low);
}
void AntTsp::antk_path(double *t,double *y,double rr,double bb,double Q,
int *path1,double *z,double *dt)
{
double a;
int sup,i,i0,ii,allowlength;
*z=0;
sup=0;
allowlength=n-1;
a=floor((double)(rand()%1000)/1000.0*n);
path1[0]=(int)a;
for(i=0;i<=path1[0]-1;i++)
allow[i]=i;
for(i=path1[0];i<=n-2;i++)
allow[i]=i+1;

while(sup<=n-2)
{
ii=path1[sup];
for(i=0;i<=allowlength-1;i++)
p[i]=pow(t[n*ii+allow[i]],rr)*pow(y[n*ii+allow[i]],bb);

i0=randchoice(p,allowlength);
path1[sup+1]=allow[i0];
(*z)+=d(ii,allow[i0]);
for(i=i0+1;i<=allowlength-1;i++)
allow[i-1]=allow[i];
allowlength--;
sup++;
}
(*z)+=d(path1[n-1],path1[0]);
for(i=0;i<=n-2;i++)
dt[n*path1[i]+path1[i+1]]+=Q/(*z);
dt[n*path1[n-1]+path1[0]]+=Q/(*z);
}


double AntTsp::Ant(int *pathmin)
{
p=new double[n];
allow=new int[n];
double tspmin;
double *t,*t1,*y,*dt,z,Q,Qmin=1.0e+300,Qmax=0.0;
int *path1;
int i,j,ii;
t=new double[n*n];
if(!t){ printf("assign fail \n"); exit(0); }
t1=new double[n*n];
y=new double[n*n];
if(y==NULL){ printf("assign fail \n"); exit(0); }
dt=new double[n*n];
if(dt==NULL){ printf("assign fail\n"); exit(0); }
path1=new int[n];
if(path1==NULL){ printf("asign fail\n"); exit(0); }

for(i=0;i<=n-1;i++)
{
t[n*i+i]=c;
y[n*i+i]=0.0;
for(j=i+1;j<=n-1;j++)
{
t[n*i+j]=t[n*j+i]=c;
y[n*i+j]=y[n*j+i]=1.0/d(i,j);
}
}
for(i=0;i<=n-1;i++)
for(j=i+1;j<=n-1;j++)
{
if(Qmax<d(i,j)) Qmax=d(i,j);
if(Qmin>d(i,j)) Qmin=d(i,j);
}
Q=c*n*(Qmin+Qmax)/2.0;
tspmin=1.0e+300;
ii=0;
for(i=0;i<=n-1;i++)
for(j=i+1;j<=n-1;j++)
dt[n*i+j]=dt[n*j+i]=0.0;
for(i=0;i<=n-1;i++)
dt[n*i+i]=0.0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
t1[n*i+j]=t[n*i+j];
while(ii<=nc-1){
for(j=0;j<=m-1;j++){
antk_path(t,y,r,b,Q,path1,&z,dt);
if(tspmin>z){
tspmin=z;
for(i=0;i<=n-1;i++)
pathmin[i]=path1[i]+1;
}
for(int i1=0;i1<=n-1;i1++)
for(int j1=0;j1<=n-1;j1++)
t1[n*i1+j1]=pp*t1[n*i1+j1]+dt[n*i1+j1];
}

for(i=0;i<=n-1;i++)
for(j=0;j<=n-1;j++)
t[n*i+j]=t1[n*i+j];
ii++;
}
delete []t;
delete []t1;
delete []y;
delete []dt;
delete []path1;
delete []p;
delete []allow;
return tspmin;
}
void main(void)
{
int i,j;
int *pathmin;
/* double dd[36]={100.0,5.0 ,4.0, 5.0, 10.0, 2.5,
5.0, 100.0,3.0, 4.0, 5.0, 3.5,
4.0, 3.0, 100.0,7.0, 4.5, 5.0,
5.0, 4.0, 7.0, 100.0,2.5, 4.0,
10.0, 5.0, 4.5, 2.5, 100.0,5.5,
2.5, 3.5, 5.0, 4.0, 5.5, 100.0
},
*/
double *dd;
double tspmin;
int n=60;
dd=new double[n*n];
for(i=0;i<n;i++){
for(j=i+1;j<n;j++)
dd[n*i+j]=dd[n*j+i]=(double)(rand()%1000)/1000.0;
dd[n*i+i]=1.0e+200;
}
pathmin=new int[n];
matrix d(dd,n,n);
AntTsp tt;
tt.d=d;
tt.n=n;
tt.m=50;
tt.nc=1;
tt.r=0.06;
tt.b=1.2;
tt.pp=0.9;
tt.c=1.0;
while(tt.nc>0){
tspmin=tt.Ant(pathmin);
printf("tspmin=%lf\n",tspmin);
printf("pathmin:\n ");
for(i=0;i<=n-1;i++)
printf("%d-->",pathmin[i]);
printf("%d\n",pathmin[0]);
cout<<"Input nc: ";
cin>>tt.nc;
}
delete []pathmin;
delete [] dd;
}


文件3：matrix.cpp

#ifndef MATRIX_CPP
#define MATRIX_CPP

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <fstream.h>
// 本程序实现matrix类
// 对matrix类的定义
#include "matrix.h"

matrix::matrix(buffer * b): // 缺省构造函数，产生0行0列空矩阵
rownum(0),colnum(0),istrans(0),isneg(0),issym(0)
{
if(b){ // 如用户给出b，则使用b作为数据缓存
buf=b;
buf->alloc(0);
}
else // 否则，产生一个新的缓存，类别是当前缺省的缓存类
buf = getnewbuffer(0);
}

matrix::matrix(size_t n, buffer * b): // 产生n阶单位方阵
rownum(n),colnum(1),istrans(0),isneg(0),issym(0)
{
if(b){ // 如用户给出b，则使用b作为数据缓存
buf=b;
buf->alloc(n);
}
else // 否则，产生一个新的缓存，类别是当前缺省的缓存类
buf = getnewbuffer(n);
}

matrix unit(size_t n)
{
if(n==0) throw TMESSAGE("n must larger than 0\n");
matrix m(n,n);
for(size_t i=0; i<n; i++)
for(size_t j=0; j<n; j++)
if(i==j) m.set(i,j,1.0);
else m.set(i,j,0.0);
m.issym = 1;
return m;
}


matrix::matrix(size_t r, size_t c, buffer * b):
rownum(r),colnum(c),istrans(0),isneg(0),issym(0)
{
if(b){ // 如用户给出b，则使用b作为数据缓存
buf=b;
buf->alloc(r*c);
}
else // 否则，产生一个新的缓存，类别是当前缺省的缓存类
buf = getnewbuffer(r*c);
}

matrix::matrix(matrix& m) // 拷贝构造函数
{
rownum = m.rownum; // 拷贝行数与列数
colnum = m.colnum;
istrans = m.istrans; // 拷贝转置标记
isneg = m.isneg; // 拷贝负值标记
issym = m.issym; // 拷贝对称标记
buf = m.buf; // 设置指向相同缓存的指针
buf->refnum++; // 缓存的引用数增1
}

matrix::matrix(const char * filename, buffer * b): // 从数据文件构造矩阵
istrans(0), isneg(0), issym(0){
buf=b;
loadfile(filename);
}

matrix::matrix(void * data, size_t r, size_t c, buffer * b):
rownum(r),colnum(c),istrans(0),isneg(0),issym(0) // 数据构造函数
{
if(b){
buf=b;
buf->alloc(r*c);
}
else
buf = getnewbuffer(r*c);
DOUBLE * d = (DOUBLE *)data;
for(size_t i=0; i<r*c; i++) // 这里进行数据拷贝，因此原数据的内存是可以释放的
buf->set(i,d[i]);
checksym(); // 检查是否为对称阵
}

void matrix::loadfile(const char * filename){// 从数据文件构造矩阵
size_t rownum1,colnum1;
int g;
istrans=0;
isneg=0;
issym=0;
char label[10];
ifstream fin(filename); // 打开文件流
fin>>g; // 读取数据文件格式数
fin >> label; // 文件开始必须有matrix关键词
if(strcmp(label, "matrix")!=0) throw -1;
fin >> rownum1 >> colnum1; // 读取行数和列数
if(!fin.good()) throw TMESSAGE("read file error!");
if(buf){
if((rownum1*colnum1)!=(rownum*colnum)){
buf->alloc(rownum1*colnum1);
rownum=rownum1;
colnum=colnum1;
}
else{
rownum=rownum1;
colnum=colnum1;
}
}
else{
buf = getnewbuffer(rownum1*colnum1);
rownum=rownum1;
colnum=colnum1;
}
if(g==2){
size_t line;
for(size_t i=0; i<rownum; i++) { // 依次按行读取
fin >> line; // 读行号
if(line != i+1) throw -1;
fin.width(sizeof(label));
fin >> label; // 行号后跟冒号
if(label[0] != ':') throw TMESSAGE("format error!");
DOUBLE a;
for(size_t j=0; j<colnum; j++) { // 读取一行数据
fin >> a;
set(i,j,a);
}
if(!fin.good()) throw TMESSAGE("read file error!");
}
}
else if(g==1){
size_t nn;
fin>>nn;
for(int i=0;i<(int)rownum;i++)
for(int j=0;j<(int)colnum;j++)
set(i,j,0);
for(size_t i=0;i<(size_t)nn;i++){
int ii,jj;
double val;
fin>>ii>>jj>>val;
set(ii-1,jj-1,val);
}
}
else{
throw -1;
}
checksym(); // 检查是否为对称阵
}

void matrix::savefile(char * filename,int g){
size_t rownum1,colnum1;
ofstream fout(filename);
rownum1=getrownum();
colnum1=getcolnum();
if(g==2){
fout<<2<<"\n";
fout<<"matrix "<<rownum1<<" "<<colnum1<<endl;
for(size_t i=1;i<=rownum1;i++){
fout<<i<<":";
size_t j;
for(j=0;j<colnum1;j++){
fout<<" "<<value(i-1,j);
if((j%10)==9) fout<<endl;
}
if(((j-1)%10)!=9) fout<<endl;
}
}
else if(g==1){
fout<<1<<"\n";
fout<<"matrix "<<rownum1<<" "<<colnum1<<endl;
size_t nn=0;
for(int i=0;i<(int)rownum1;i++)
for(int j=0;j<(int)colnum1;j++)
if(value(i,j)!=0.0) nn++;
fout<<nn<<"\n";
for(int i=0;i<(int)rownum1;i++)
for(int j=0;j<(int)colnum1;j++){
if(value(i,j)!=0.0)
fout<<(i+1)<<" "<<(j+1)<<" "<<value(i,j);
}
}
else{
throw -1;
}
}

DOUBLE & matrix::operator()(size_t r, size_t c)
{
if(r>= rownum || c>= colnum)
throw TMESSAGE("Out range!");
return value(r,c);
}

matrix matrix::operator=(matrix m) // 赋值重载
{
rownum = m.rownum; // 行数和列数的拷贝
colnum = m.colnum;
istrans = m.istrans; // 转置标志的拷贝
isneg = m.isneg; // 取负标志的拷贝
issym = m.issym; // 对称标志的拷贝
if(buf == m.buf) // 如果原缓存与m的缓存一样，则返回
return (*this);
buf->refnum--; // 原缓存不同，则原缓存的引用数减1
if(!buf->refnum)delete buf; // 减1后的原缓存如果引用数为0，则删除原缓存
buf = m.buf; // 将原缓存指针指向m的缓存
buf->refnum++; // 新缓存的引用数增1
checksym(); // 检查是否为对称阵
return (*this); // 返回自己的引用
}

matrix matrix::operator=(DOUBLE a) // 通过赋值运算符将矩阵所有元素设为a
{
if(rownum == 0 || colnum == 0) return (*this);
for(size_t i=0; i<rownum; i++)
for(size_t j=0; j<colnum; j++)
set(i,j,a);
if(rownum == colnum)issym = 1;
return (*this);
}

matrix matrix::operator-() // 矩阵求负，产生负矩阵
{
matrix mm(*this);
mm.neg();
return mm;
}

ostream& operator<<(ostream& o, matrix& m) // 流输出运算符
{
// 先输出关键字matrix,然后是行号，制表符，列号，换行
o << "matrix " << m.rownum << '\t' << m.colnum << endl;
for(size_t i=0; i<m.rownum; i++) { // 依次输出各行
o<< (i+1) <<':'; // 行号后跟冒号。注意输出的行号是从1开始
// 是内部编程的行号加1
size_t k=8; // 后面输入一行的内容，每次一行数字超过八个时
// 就换一行
for(size_t j=0; j<m.colnum; j++) {
o<<'\t'<<m.value(i,j);
if(--k==0) {
k=8;
o<<endl;
}
}
o<<endl; // 每行结束时也换行
}
return o;
}

istream& operator>>(istream& in, matrix& m) // 流输入运算符
{
char label[10];
in.width(sizeof(label));
in >> label; // 输入matrix关键字
if(strcmp(label, "matrix")!=0)
throw TMESSAGE("format error!\n");
in >> m.rownum >> m.colnum; // 输入行号和列号
if(!in.good()) throw TMESSAGE("read file error!");
m.buf->refnum--; // 原缓存引用数减1
if(!m.buf->refnum) delete m.buf; // 如原缓存引用数为0，则删去原缓存
m.isneg = m.istrans = 0; // 转置和取负标志清零
m.buf = getnewbuffer(m.rownum*m.colnum); // 按缺省子类产生新缓存
size_t line; // 按矩阵行输入
for(size_t i=0; i<m.rownum; i++) {
in >> line; // 先输入行号
if(line != i+1) throw TMESSAGE("format error!\n");
in.width(sizeof(label)); // 行号后应跟冒号
in >> label;
if(label[0] != ':') throw TMESSAGE("format error!\n");
DOUBLE a; // 随后是本行的各个数值
for(size_t j=0; j<m.colnum; j++) {
in >> a;
m.set(i,j,a);
}
if(!in.good()) throw TMESSAGE("read file error!");
}
m.checksym(); // 检查是否为对称阵
return in;
}

matrix matrix::operator*=(DOUBLE a) // 矩阵数乘常数a，结果放在原矩阵
{
for(size_t i=0; i<rownum; i++)
for(size_t j=0; j<colnum; j++)
set(i,j,a*value(i,j));
return (*this);
}

matrix matrix::operator*(DOUBLE a) // 矩阵数乘常数a，原矩阵内容不变，返回一新矩阵
{
matrix m(rownum, colnum);
for(size_t i=0; i<rownum; i++)
for(size_t j=0; j<colnum; j++)
m.set(i,j,a*value(i,j));
return m;
}

matrix matrix::operator+(matrix& m) // 矩阵相加，产生一新的矩阵并返回它
{
if(rownum != m.rownum || colnum != m.colnum) // 对应行列必须相同
throw TMESSAGE("can not do add of matrix\n");
matrix mm(rownum, colnum); // 产生一同自己同形的矩阵
DOUBLE a;
for(size_t i=0; i<rownum; i++) // 求和
for(size_t j=0; j<colnum; j++)
{
a = value(i,j)+m.value(i,j);
mm.set(i,j,a);
}
mm.checksym(); // 检查是否为对称阵
return mm;
}

matrix matrix::operator+=(matrix m) // 矩阵求和，自己内容改变为和
{
DOUBLE a;
for(size_t i=0; i<rownum; i++)
for(size_t j=0; j<colnum; j++)
{
a = value(i,j)+m.value(i,j);
set(i,j,a);
}
checksym(); // 检查是否为对称阵
return (*this);
}

matrix matrix::operator+(DOUBLE a) // 矩阵加常数，指每一元素加一固定的常数，产生
// 新矩阵，原矩阵不变
{
matrix m(rownum, colnum);
for(size_t i=0; i<rownum; i++)
for(size_t j=0; j<colnum; j++)
m.set(i,j,a+value(i,j));
return m;
}

matrix matrix::operator+=(DOUBLE a) // 矩阵自身加常数，自身内容改变
{
for(size_t i=0; i<rownum; i++)
for(size_t j=0; j<colnum; j++)
set(i,j,a+value(i,j));
return (*this);
}

matrix operator-(DOUBLE a, matrix& m) { // 常数减矩阵，产生新的矩阵
return (-m)+a;
};


matrix matrix::operator-(matrix m) // 矩阵相减，产生新的矩阵
{
matrix mm(*this); // 产生一同自己同形的矩阵
//mm += -m; // 加上相应的负矩阵
matrix mm1;
mm1 = mm.operator-();
mm+= mm1;
return mm;
}

matrix matrix::operator-=(matrix m) // 矩阵相减，结果修改原矩阵
{
(*this) += (-m);
return (*this);
}

matrix matrix::operator*(matrix m) // 矩阵相乘，原矩阵内容不变，产生一新矩阵
{
if(colnum != m.rownum) // 必须满足相乘条件
throw TMESSAGE("can not multiply!");
matrix mm(rownum,m.colnum); // 计算并产生一合要求的矩阵放乘积
DOUBLE a;
for(size_t i=0; i<rownum; i++) // 计算乘积
for(size_t j=0; j<m.colnum; j++){
a = 0.0;
for(size_t k=0; k<colnum; k++)
a += value(i,k)*m.value(k,j);
mm.set(i,j,a);
}
mm.checksym(); // 检查是否为对称阵
return mm; // 返回乘积
}

matrix matrix::operator*=(matrix m) // 矩阵相乘，自己修改成积矩阵
{
(*this) = (*this)*m;
return (*this);
}

matrix matrix::t() // 矩阵转置，产生新的矩阵
{
matrix mm(*this);
mm.trans();
return mm;
}

int matrix::isnear(matrix m, double e) // 检查二矩阵是否近似相等
{
if(rownum != m.rownum || colnum != m.colnum) return 0;
for(size_t i=0; i< rownum; i++)
for(size_t j=0; j< colnum; j++)
if(fabs(value(i,j)-m.value(i,j)) > e) return 0;
return 1;
}

int matrix::isnearunit(double e) // 检查矩阵是否近似为单位矩阵
{
if(rownum != colnum) return 0;
return isnear(unit(rownum), e);
}

matrix matrix::row(size_t r) // 提取第r行行向量
{
matrix mm(1, colnum);
for(size_t i=0; i< colnum; i++)
mm.set(0, i, value(r,i));
return mm;
}

matrix matrix::col(size_t c) // 提取第c列列向量
{
matrix mm(rownum, 1);
for(size_t i=0; i< rownum; i++)
mm.set(i, value(i, c));
return mm;
}

void matrix::swapr(size_t r1, size_t r2, size_t k) // 交换矩阵r1和r2两行
{
DOUBLE a;
for(size_t i=k; i<colnum; i++) {
a = value(r1, i);
set(r1, i, value(r2, i));
set(r2, i, a);
}
}

void matrix::swapc(size_t c1, size_t c2, size_t k) // 交换c1和c2两列
{
DOUBLE a;
for(size_t i=k; i<colnum; i++) {
a = value(i, c1);
set(i, c1, value(i, c2));
set(i, c2, a);
}
}

DOUBLE matrix::maxabs(size_t &r, size_t &c, size_t k) // 求第k行和第k列后的主元及位置
{
DOUBLE a=0.0;
for(size_t i=k;i<rownum;i++)
for(size_t j=k;j<colnum;j++)
if(a < fabs(value(i,j))) {
r=i;c=j;a=fabs(value(i,j));
}
return a;
}

size_t matrix::zgsxy(matrix & m, int fn) // 进行主高斯消元运算，fn为参数，缺省为0
/* 本矩阵其实是常数阵，而矩阵m必须是方阵
运算过程其实是对本矩阵和m同时作行初等变换，
运算结果m的对角线相乘将是行列式，而本矩阵则变换成
自己的原矩阵被m的逆阵左乘，m的秩被返回，如果秩等于阶数
则本矩阵中的内容已经是唯一解
*/
{
if(rownum != m.rownum || m.rownum != m.colnum) // 本矩阵行数必须与m相等
// 且m必须是方阵
throw TMESSAGE("can not divide!");
lbuffer * bb = getnewlbuffer(rownum); // 产生一维数为行数的长整数缓存区
lbuffer & b = (*bb); // 用引用的办法使下面的程序容易懂
size_t is;
DOUBLE a;
size_t i,j,rank=0;
size_t k;
for(k=0; k<rownum; k++) { // 从第0行到第k行进行主高斯消元
if(m.maxabs(is, i, k)==0) // 求m中第k级主元，主元所在的行，列存在is,i中
break; // 如果主元为零，则m不可逆，运算失败
rank = k+1; // rank存放当前的阶数
b.retrieve(k) = i; // 将长整数缓存区的第k个值设为i
if(i != k)
m.swapc(k, i); // 交换m中i,k两列
if(is != k) {
m.swapr(k, is, k); // 交换m中i,k两行,从k列以后交换
swapr(k, is); // 交换本矩阵中i,k两行
}
a = m.value(k,k); // 取出主元元素
for (j=k+1;j<rownum;j++) // 本意是将m的第k行除以主元
// 但只需把第k行的第k+1列以上除以主元即可
// 这样还保留了主元作行列式运算用
m.set(k,j,m.value(k,j)/a);
for (j=0;j<colnum;j++) // 将本矩阵的第k行除以主元
set(k,j,value(k,j)/a);
// 上面两步相当于将m和本矩阵构成的增广矩阵第k行除以主元
// 下面对增广矩阵作行基本初等变换使第k行的其余列均为零
// 但0值无必要计算，因此从第k+1列开始计算
for(j=k+1;j<rownum;j++) // j代表列，本矩阵的行数就是m的列数
for(i=0;i<rownum;i++) //i代表行，依次对各行计算，k行除外
if(i!=k)
m.set(i,j,m.value(i,j)-m.value(i,k)*m.value(k,j));
// 对本矩阵亦作同样的计算
for(j=0;j<colnum;j++)
for(i=0;i<rownum;i++)
if(i!=k)
set(i,j,value(i,j)-m.value(i,k)*value(k,j));
} // 主高斯消元循环k结束
if(fn == 1) {
for(j=0; j<rank; j++)
for(i=0; i<rownum; i++)
if(i==j) m.set(i,i,1.0);
else
m.set(i,j,0.0);
for(k = rank; k>0; k--)
m.swapc(k-1,(size_t)b[k-1]);
}
for(k = rank; k>0; k--) // 将本矩阵中的各行按b中内容进行调节
if(b[k-1] != k-1)
swapr(k-1,(size_t)b[k-1]); // 行交换
delete bb; // 释放长整数缓存
return rank; // 返回mm的秩
}

matrix matrix::operator/=(matrix m) // 利用重载的除法符号/=来解方程
// 本矩阵设为d,则方程为mx=d,考虑解写成x=d/m的形式，
// 而方程的解也存放在d中，则实际编程时写d/=m
{
if(zgsxy(m)!=rownum) // 如秩不等于m的阶数，则方程无解
throw TMESSAGE("can not divide!");
return *this;
}

matrix matrix::operator/(matrix m) // 左乘m的逆矩阵产生新矩阵
{
m.inv(); // m的逆矩阵
return (*this)*m;
}

matrix matrix::inv() // 用全选主元高斯-约当法求逆矩阵
{
if(rownum != colnum || rownum == 0)
throw TMESSAGE("Can not calculate inverse");
size_t i,j,k;
DOUBLE d,p;
lbuffer * isp = getnewlbuffer(rownum); // 产生一维数为行数的长整数缓存区
lbuffer * jsp = getnewlbuffer(rownum); // 产生一维数为行数的长整数缓存区
lbuffer& is = *isp; // 使用引用使程序看起来方便
lbuffer& js = *jsp;
for(k=0; k<rownum; k++)
{
d = maxabs(i, j, k); // 全主元的位置和值
is[k] = i;
js[k] = j;
if(d==0.0) {
delete isp;
delete jsp;
throw TMESSAGE("can not inverse");
}
if (is[k] != k) swapr(k,(size_t)is[k]);
if (js[k] != k) swapc(k,(size_t)js[k]);
p = 1.0/value(k,k);
set(k,k,p);
for (j=0; j<rownum; j++)
if (j!=k) set(k,j,value(k,j)*p);
for (i=0; i<rownum; i++)
if (i!=k)
for (j=0; j<rownum; j++)
if (j!=k) set(i,j,value(i,j)-value(i,k)*value(k,j));
for (i=0; i<rownum; i++)
if (i!=k) set(i,k,-value(i,k)*p);
} // end for k
for (k=rownum; k>0; k--)
{ if (js[k-1]!=(k-1)) swapr((size_t)js[k-1], k-1);
if (is[k-1]!=k-1) swapc((size_t)is[k-1], k-1);
}
delete isp;
delete jsp;
checksym(); // 检查是否为对称阵
return (*this);
}

matrix matrix::operator~() // 求逆矩阵，但产生新矩阵
{
matrix m(*this);
m.inv();
return m;
}

matrix operator/(DOUBLE a, matrix& m) // 求逆矩阵再乘常数
{
matrix mm(m);
mm.inv();
if(a != 1.0) mm*=a;
return mm;
}

matrix matrix::operator/=(DOUBLE a) // 所有元素乘a的倒数，自身改变
{
return operator*=(1/a);
}

matrix matrix::operator/(DOUBLE a) // 所有元素乘a的倒数，产生新的矩阵
{
matrix m(*this);
m/=a;
return m;
}

DOUBLE matrix::det(DOUBLE err) // 求行列式的值
{
if(rownum != colnum || rownum == 0)
throw TMESSAGE("Can not calculate det");
matrix m(*this);
size_t rank;
return m.detandrank(rank, err);
}

size_t matrix::rank(DOUBLE err) // 求矩阵的秩
{
if(rownum==0 || colnum==0)return 0;
size_t k;
k = rownum > colnum ? colnum : rownum;
matrix m(k,k); // 产生k阶方阵
for(size_t i=0; i<k; i++)
for(size_t j=0; j<k; j++)
m.set(i,j,value(i,j));
size_t r;
m.detandrank(r, err);
return r;
}

DOUBLE matrix::detandrank(size_t & rank, DOUBLE err) // 求方阵的行列式和秩
{
if(rownum != colnum || rownum == 0)
throw TMESSAGE("calculate det and rank error!");
size_t i,j,k,is,js;
double f,detv,q,d;
f=1.0; detv=1.0;
rank = 0;
for (k=0; k<rownum-1; k++)
{
q = maxabs(is, js, k);
if(q<err) return 0.0; // 如主元太小，则行列式的值被认为是0
rank++; // 秩增1
if(is!=k) { f=-f; swapr(k,is,k); }
if(js!=k) { f=-f; swapc(k,js,k); }
q = value(k,k);
detv *= q;
for (i=k+1; i<rownum; i++)
{
d=value(i,k)/q;
for (j=k+1; j<rownum; j++)
set(i,j,value(i,j)-d*value(k,j));
}
} // end loop k
q = value(rownum-1,rownum-1);
if(q != 0.0 ) rank++;
return f*detv*q;
}

void matrix::checksym() // 检查和尝试调整矩阵到对称矩阵
{
issym = 0; // 先假设矩阵非对称
if(rownum != colnum) return; // 行列不等当然不是对称矩阵
DOUBLE a,b;
for(size_t i=1; i<rownum; i++) // 从第二行开始检查
for(size_t j=0; j<i; j++) // 从第一列到第i-1列
{
a = value(i,j);
b = value(j,i);
if( fabs(a-b) >= defaulterr ) return; // 发现不对称，返回
if(a!=b)set(i,j,b); // 差别很小就进行微调
}
issym = 1; // 符合对称阵标准
}

void matrix::house(buffer & b, buffer & c)// 用豪斯荷尔德变换将对称阵变为对称三对
// 角阵，b返回主对角线元素，c返回次对角线元素
{
if(!issym) throw TMESSAGE("not symatry");
size_t i,j,k;
DOUBLE h,f,g,h2,a;
for (i=rownum-1; i>0; i--)
{ h=0.0;
if (i>1)
for (k=0; k<i; k++)
{ a = value(i,k); h += a*a;}
if (h == 0.0)
{ c[i] = 0.0;
if (i==1) c[i] = value(i,i-1);
b[i] = 0.0;
}
else
{ c[i] = sqrt(h);
a = value(i,i-1);
if (a > 0.0) c[i] = -c[i];
h -= a*c[i];
set(i,i-1,a-c[i]);
f=0.0;
for (j=0; j<i; j++)
{ set(j,i,value(i,j)/h);
g=0.0;
for (k=0; k<=j; k++)
g += value(j,k)*value(i,k);
if(j<=i-2)
for (k=j+1; k<i; k++)
g += value(k,j)*value(i,k);
c[j] = g/h;
f += g*value(j,i);
}
h2=f/(2*h);
for (j=0; j<i; j++)
{ f=value(i,j);
g=c[j] - h2*f;
c[j] = g;
for (k=0; k<=j; k++)
set(j,k, value(j,k)-f*c[k]-g*value(i,k) );
}
b[i] = h;
}
}
for (i=0; i<=rownum-2; i++) c[i] = c[i+1];
c[rownum-1] = 0.0;
b[0] = 0.0;
for (i=0; i<rownum; i++)
{ if ((b[i]!=0.0)&&(i>=1))
for (j=0; j<i; j++)
{ g=0.0;
for (k=0; k<i; k++)
g=g+value(i,k)*value(k,j);
for (k=0; k<i; k++)
set(k,j,value(k,j)-g*value(k,i));
}
b[i] = value(i,i);
set(i,i,1.0);
if (i>=1)
for (j=0; j<=i-1; j++)
{ set(i,j,0.0);
set(j,i,0.0); }
}
return;
}

void matrix::trieigen(buffer& b, buffer& c, size_t l, DOUBLE eps)
// 计算三对角阵的全部特征值与特征向量
{ size_t i,j,k,m,it;
double d,f,h,g,p,r,e,s;
c[rownum-1]=0.0; d=0.0; f=0.0;
for (j=0; j<rownum; j++)
{ it=0;
h=eps*(fabs(b[j])+fabs(c[j]));
if (h>d) d=h;
m=j;
while ((m<rownum)&&(fabs(c[m])>d)) m+=1;
if (m!=j)
{ do
{ if (it==l) throw TMESSAGE("fial to calculate eigen value");
it += 1;
g=b[j];
p=(b[j+1]-g)/(2.0*c[j]);
r=sqrt(p*p+1.0);
if (p>=0.0) b[j]=c[j]/(p+r);
else b[j]=c[j]/(p-r);
h=g-b[j];
for (i=j+1; i<rownum; i++)
b[i]-=h;
f=f+h; p=b[m]; e=1.0; s=0.0;
for (i=m-1; i>=j; i--)
{ g=e*c[i]; h=e*p;
if (fabs(p)>=fabs(c[i]))
{ e=c[i]/p; r=sqrt(e*e+1.0);
c[i+1]=s*p*r; s=e/r; e=1.0/r;
}
else
{ e=p/c[i]; r=sqrt(e*e+1.0);
c[i+1]=s*c[i]*r;
s=1.0/r; e=e/r;
}
p=e*b[i]-s*g;
b[i+1]=h+s*(e*g+s*b[i]);
for (k=0; k<rownum; k++)
{
h=value(k,i+1);
set(k,i+1, s*value(k,i)+e*h);;
set(k,i,e*value(k,i)-s*h);
}
if(i==0) break;
}
c[j]=s*p; b[j]=e*p;
}
while (fabs(c[j])>d);
}
b[j]+=f;
}
for (i=0; i<=rownum; i++)
{ k=i; p=b[i];
if (i+1<rownum)
{ j=i+1;
while ((j<rownum)&&(b[j]<=p))
{ k=j; p=b[j]; j++;}
}
if (k!=i)
{ b[k]=b[i]; b[i]=p;
for (j=0; j<rownum; j++)
{ p=value(j,i);
set(j,i,value(j,k));
set(j,k,p);
}
}
}
}

matrix matrix::eigen(matrix & eigenvalue, DOUBLE eps, size_t l)
// 计算对称阵的全部特征向量和特征值
// 返回按列排放的特征向量，而eigenvalue将返回一维矩阵，为所有特征值
// 组成的列向量
{
if(!issym) throw TMESSAGE("it is not symetic matrix");
eigenvalue = matrix(rownum); // 产生n行1列向量准备放置特征值
matrix m(*this); // 复制自己产生一新矩阵
if(rownum == 1) { // 如果只有1X1矩阵，则特征向量为1，特征值为value(0,0)
m.set(0,0,1.0);
eigenvalue.set(0,value(0,0));
return m;
}
buffer * b, *c;
b = getnewbuffer(rownum);
c = getnewbuffer(rownum);
m.house(*b,*c); // 转换成三对角阵
m.trieigen(*b,*c,l,eps); // 算出特征向量和特征值
for(size_t i=0; i<rownum; i++) // 复制b的内容到eigenvalue中
eigenvalue.set(i,(*b)[i]);
return m;
}

void matrix::hessenberg() // 将一般实矩阵约化为赫申伯格矩阵
{
size_t i,j,k;
double d,t;
for (k=1; k<rownum-1; k++)
{ d=0.0;
for (j=k; j<rownum; j++)
{ t=value(j,k-1);
if (fabs(t)>fabs(d))
{ d=t; i=j;}
}
if (fabs(d)!=0.0)
{ if (i!=k)
{ for (j=k-1; j<rownum; j++)
{
t = value(i,j);
set(i,j,value(k,j));
set(k,j,t);
}
for (j=0; j<rownum; j++)
{
t = value(j,i);
set(j,i,value(j,k));
set(j,k,t);
}
}
for (i=k+1; i<rownum; i++)
{
t = value(i,k-1)/d;
set(i,k-1,0.0);
for (j=k; j<rownum; j++)
set(i,j,value(i,j)-t*value(k,j));
for (j=0; j<rownum; j++)
set(j,k,value(j,k)+t*value(j,i));
}
}
}
}

void matrix::qreigen(matrix & u1, matrix & v1, size_t jt, DOUBLE eps)
// 求一般实矩阵的所有特征根
// a和b均返回rownum行一列的列向量矩阵，返回所有特征根的实部和虚部
{
matrix a(*this);
a.hessenberg(); // 先算出赫申伯格矩阵
u1 = matrix(rownum);
v1 = matrix(rownum);
buffer * uu = getnewbuffer(rownum);
buffer * vv = getnewbuffer(rownum);
buffer &u = *uu;
buffer &v = *vv;
size_t m,it,i,j,k,l;
size_t iir,iic,jjr,jjc,kkr,kkc,llr,llc;
DOUBLE b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
it=0; m=rownum;
while (m!=0)
{ l=m-1;
while ( l>0 && (fabs(a.value(l,l-1))>eps*
(fabs(a.value(l-1,l-1))+fabs(a.value(l,l))))) l--;
iir = m-1; iic = m-1;
jjr = m-1; jjc = m-2;
kkr = m-2; kkc = m-1;
llr = m-2; llc = m-2;
if (l==m-1)
{ u[m-1]=a.value(m-1,m-1); v[m-1]=0.0;
m--; it=0;
}
else if (l==m-2)
{ b=-(a.value(iir,iic)+a.value(llr,llc));
c=a.value(iir,iic)*a.value(llr,llc)-
a.value(jjr,jjc)*a.value(kkr,kkc);
w=b*b-4.0*c;
y=sqrt(fabs(w));
if (w>0.0)
{ xy=1.0;
if (b<0.0) xy=-1.0;
u[m-1]=(-b-xy*y)/2.0;
u[m-2]=c/u[m-1];
v[m-1]=0.0; v[m-2]=0.0;
}
else
{ u[m-1]=-b/2.0; u[m-2]=u[m-1];
v[m-1]=y/2.0; v[m-2]=-v[m-1];
}
m=m-2; it=0;
}
else
{
if (it>=jt) {
delete uu;
delete vv;
throw TMESSAGE("fail to calculate eigenvalue");
}
it++;
for (j=l+2; j<m; j++)
a.set(j,j-2,0.0);
for (j=l+3; j<m; j++)
a.set(j,j-3,0.0);
for (k=l; k+1<m; k++)
{ if (k!=l)
{ p=a.value(k,k-1); q=a.value(k+1,k-1);
r=0.0;
if (k!=m-2) r=a.value(k+2,k-1);
}
else
{
x=a.value(iir,iic)+a.value(llr,llc);
y=a.value(llr,llc)*a.value(iir,iic)-
a.value(kkr,kkc)*a.value(jjr,jjc);
iir = l; iic = l;
jjr = l; jjc = l+1;
kkr = l+1; kkc = l;
llr = l+1; llc = l+1;
p=a.value(iir,iic)*(a.value(iir,iic)-x)
+a.value(jjr,jjc)*a.value(kkr,kkc)+y;
q=a.value(kkr,kkc)*(a.value(iir,iic)+a.value(llr,llc)-x);
r=a.value(kkr,kkc)*a.value(l+2,l+1);
}
if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
{ xy=1.0;
if (p<0.0) xy=-1.0;
s=xy*sqrt(p*p+q*q+r*r);
if (k!=l) a.set(k,k-1,-s);
e=-q/s; f=-r/s; x=-p/s;
y=-x-f*r/(p+s);
g=e*r/(p+s);
z=-x-e*q/(p+s);
for (j=k; j<m; j++)
{
iir = k; iic = j;
jjr = k+1; jjc = j;
p=x*a.value(iir,iic)+e*a.value(jjr,jjc);
q=e*a.value(iir,iic)+y*a.value(jjr,jjc);
r=f*a.value(iir,iic)+g*a.value(jjr,jjc);
if (k!=m-2)
{ kkr = k+2; kkc = j;
p=p+f*a.value(kkr,kkc);
q=q+g*a.value(kkr,kkc);
r=r+z*a.value(kkr,kkc);
a.set(kkr,kkc,r);
}
a.set(jjr,jjc,q);
a.set(iir,iic,p);
}
j=k+3;
if (j>=m-1) j=m-1;
for (i=l; i<=j; i++)
{ iir = i; iic = k;
jjr = i; jjc = k+1;
p=x*a.value(iir,iic)+e*a.value(jjr,jjc);
q=e*a.value(iir,iic)+y*a.value(jjr,jjc);
r=f*a.value(iir,iic)+g*a.value(jjr,jjc);
if (k!=m-2)
{ kkr = i; kkc = k+2;
p=p+f*a.value(kkr,kkc);
q=q+g*a.value(kkr,kkc);
r=r+z*a.value(kkr,kkc);
a.set(kkr,kkc,r);
}
a.set(jjr,jjc,q);
a.set(iir,iic,p);
}
}
}
}
}
for(i=0;i<rownum;i++) {
u1.set(i,u[i]);
v1.set(i,v[i]);
}
delete uu;
delete vv;
}

DOUBLE gassrand(int rr) // 返回一零均值单位方差的正态分布随机数
{
static DOUBLE r=3.0; // 种子
if(rr) r = rr;
int i,m;
DOUBLE s,w,v,t;
s=65536.0; w=2053.0; v=13849.0;
t=0.0;
for (i=1; i<=12; i++)
{ r=r*w+v; m=(int)(r/s);
r-=m*s; t+=r/s;
}
t-=6.0;
return(t);
}

gassvector::gassvector(matrix & r): //r必须是正定对称阵，为正态随机向量的协方差
a(r)
{
if(r.rownum != r.colnum) throw TMESSAGE("must be a sqare matrix");
if(!r.issym) throw TMESSAGE("must be a symetic matrix");
matrix evalue;
a = a.eigen(evalue);
for(size_t i=0; i<a.colnum; i++) {
DOUBLE e = sqrt(evalue(i));
for(size_t j=0; j<a.rownum; j++)
a.set(j,i,a.value(j,i)*e);
}
}

matrix gassvector::operator()(matrix & r) // 返回给定协方差矩阵的正态随机向量
{
a = r;
matrix evalue;
a = a.eigen(evalue);
size_t i;
for(i=0; i<a.colnum; i++) {
DOUBLE e = sqrt(evalue(i));
for(size_t j=0; j<a.rownum; j++)
a.set(j,i,a.value(j,i)*e);
}
matrix rr(a.rownum); // 产生列向量
for(i=0; i<a.rownum; i++)
rr.set(i,gassrand());
return a*rr;
}

matrix gassvector::operator()() // 返回已设定的协方差矩阵的正态随机向量
{
matrix rr(a.rownum);
for(size_t i=0; i<a.rownum; i++)
rr.set(i,gassrand());
return a*rr;
}

void gassvector::setdata(matrix & r) // 根据协方差矩阵设置增益矩阵
{
if(!r.issym) throw TMESSAGE("r must be symetic!");
a = r;
matrix evalue;
a = a.eigen(evalue);
for(size_t i=0; i<a.colnum; i++) {
if(evalue(i)<0.0) throw TMESSAGE("var matrix not positive!");
DOUBLE e = sqrt(evalue(i));
for(size_t j=0; j<a.rownum; j++)
a.set(j,i,a.value(j,i)*e);
}
}

int matrix::ispositive() // 判定矩阵是否对称非负定阵，如是返回1，否则返回0
{
if(!issym) return 0;
matrix ev;
eigen(ev);
for(size_t i=0; i<rownum; i++)
if(ev(i)<0) return 0;
return 1;
}

matrix matrix::cholesky(matrix& dd) // 用乔里斯基分解法求对称正定阵的线性
// 方程组ax=d返回方程组的解，本身为a，改变为分解阵u,d不变
{
if(!issym) throw TMESSAGE("not symetic!");
if(dd.rownum != colnum) throw TMESSAGE("dd's rownum not right!");
matrix md(dd);
size_t i,j,k,u,v;
if(value(0,0)<=0.0) throw TMESSAGE("not positive");
set(0,0,sqrt(value(0,0))); // a[0]=sqrt(a[0]);
buffer& a = (*buf);
buffer& d = (*(md.buf));
size_t n = rownum;
size_t m = dd.colnum;
for (j=1; j<n; j++) a[j]=a[j]/a[0];
for (i=1; i<n; i++)
{ u=i*n+i;
for (j=1; j<=i; j++)
{ v=(j-1)*n+i;
a[u]=a[u]-a[v]*a[v];
}
if (a[u]<=0.0) throw TMESSAGE("not positive");
a[u]=sqrt(a[u]);
if (i!=(n-1))
{ for (j=i+1; j<n; j++)
{ v=i*n+j;
for (k=1; k<=i; k++)
a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
a[v]=a[v]/a[u];
}
}
}
for (j=0; j<m; j++)
{ d[j]=d[j]/a[0];
for (i=1; i<=n-1; i++)
{ u=i*n+i; v=i*m+j;
for (k=1; k<=i; k++)
d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
d[v]=d[v]/a[u];
}
}
for (j=0; j<=m-1; j++)
{ u=(n-1)*m+j;
d[u]=d[u]/a[n*n-1];
for (k=n-1; k>=1; k--)
{ u=(k-1)*m+j;
for (i=k; i<=n-1; i++)
{ v=(k-1)*n+i;
d[u]=d[u]-a[v]*d[i*m+j];
}
v=(k-1)*n+k-1;
d[u]=d[u]/a[v];
}
}
if(n>1)
for(j=1; j<n; j++)
for(i=0; i<j; i++)
a[i+j*n]=0.0;
return md;
}

DOUBLE lineopt(matrix& aa,matrix& bb, matrix& cc, matrix & xx)
// 线性规划最优点寻找程序，aa为mXn不等式约束条件左端系数矩阵，bb为不等式约束
// 条件的右端项，为m维向量，cc为目标函数系数，n维向量，xx返回极小点，为n维向量
{
if(aa.rownum != bb.rownum || aa.colnum != cc.rownum ||
aa.colnum != xx.rownum) throw TMESSAGE("dimenstion not right!");
size_t n=aa.colnum, m=aa.rownum;
size_t i,mn,k,j;
matrix a(m,n+m);
for(i=0;i<m;i++) {
for(j=0;j<n;j++)
a.set(i,j,aa(i,j));
for(j=n;j<n+m;j++)
if(j-n == i) a.set(i,j,1.0);
else a.set(i,j,0.0);
}
matrix c(m+n);
c = 0.0;
for(i=0;i<m;i++)
c.set(i,cc(i));
lbuffer* jjs = getnewlbuffer(m);
lbuffer& js = (*jjs);
DOUBLE s,z,dd,y; //,*p,*d;

for (i=0; i<m; i++) js[i]=n+i;
matrix p(m,m);
matrix d;
mn=m+n;
matrix x(mn);
while (1)
{ for (i=0; i<m; i++)
for (j=0; j<m; j++)
p.set(i,j,a(i,(size_t)js[j]));
p.inv();
d = p*a;
x = 0.0;
for (i=0; i<m; i++)
{ s=0.0;
for (j=0; j<=m-1; j++)
s+=p(i,j)*bb(j);
x.set((size_t)js[i],s);
}
k=mn; dd=1.0e-35;
for (j=0; j<mn; j++)
{ z=0.0;
for (i=0; i<m; i++)
z+=c((size_t)js[i])*d(i,j);
z-=c(j);
if (z>dd) { dd=z; k=j;}
}
if (k==mn)
{ s=0.0;
for (j=0; j<n; j++) {
xx.set(j,x(j));
s+=c(j)*x(j);
}
delete jjs;
return s;
}
j=m;
dd=1.0e+20;
for (i=0; i<=m-1; i++)
if (d(i,k)>=1.0e-20) // [i*mn+k]>=1.0e-20)
{ y=x(size_t(js[i]))/d(i,k);
if (y<dd) { dd=y; j=i;}
}
if (j==m) { delete jjs;
throw TMESSAGE("lineopt failed!");
}
js[j]=k;
}
}
matrix matrix::mean(void){
matrix mat(1,getcolnum());
double sum,rownum;
rownum=(double)getrownum();
for(int j=0;j<(int)getcolnum();j++){
sum=0.0;
for(int i=0;i<(int)getrownum();i++)
sum+=value(i,j);
mat.set(0,j,sum/rownum);
}
return mat;
}
matrix matrix::cov(void){
matrix mat_cov(getcolnum(),getcolnum());
matrix mat_mean;
mat_mean=mean();
double sum,rownum;
rownum=(double)getrownum();
for(int i=0;i<(int)getcolnum();i++)
for(int j=i;j<(int)getcolnum();j++){
sum=0.0;
for(int k=0;k<(int)getrownum();k++)
sum+=value(k,i)*value(k,j);
sum-=(rownum*mat_mean(0,i)*mat_mean(0,j));
sum/=(rownum-1);
mat_cov.set(i,j,sum);
mat_cov.set(j,i,sum);
}
return mat_cov;
}
#endif //MATRIX_CPP


文件4：matrix.h

#ifndef MATRIX_H

#define MATRIX_H

// 定义适合数学运算的实数matrix类，数据将存放在buffer类中

#include <iostream.h>
#include <iomanip.h>
#include <stdio.h>

// buffer.h包含实数缓存类buffer的定义
#include "buffer.cpp"

// 实数矩阵类，每个矩阵有rownum行，colnum列，总共rownum乘colnum个实数
// 下标一律从0开始起算，即存在第0行第0列
class matrix {
public:
buffer * buf; // 缓存区指针，指向存放实数的地方，可以是任何媒体
// buffer类是一抽象类，具体的子类可以将数据存放在内存或者
// 临时文件中，用户还可以设计自己的缓存区子类
size_t rownum, colnum; // 矩阵的行数和列数
unsigned istrans:1; // 转置标志，0为缺省，1为转置
unsigned isneg:1; // 为负标志，0为非负，1为取负
unsigned issym:1; // 对称标志，0为非对称，1为对称，不具有约束力
// 在下面的所有构造函数中，缓存区指针b通常不设，缺省为0，这样程序将根据
// 自动申请缺省的缓存区类型。如果用户定义了特殊的缓顾虑区类型，可以
// 先申请后在构造函数中指明
matrix(buffer * b=0); // 缺省构造函数，产生0行0列空矩阵
matrix(size_t n, buffer * b=0); // 产生nX1列向量
matrix(size_t r, size_t c, buffer * b=0); // 构造函数，产生r行c列矩阵
matrix(matrix& m); //拷贝构造函数，产生的矩阵内容同m相同，
// 新产生的矩阵并不申请新的缓存，而是指向与m同样的缓存，
// 只是它们共同指向的缓存的引用数增1，这样可以节省存储资源
matrix(const char * filename, buffer * b=0); /* 文件构造函数，
使用数据文件对矩阵初使化，数据文件为文本文件，格式：
以小写的matrix关键词打头，然后是空格，再跟矩阵的行数，再跟
空格，再跟列数，再跟空格，接着是内容，
内容按行分，首先是行号，从1开始计数，行号后面紧跟着冒号：，
然后是这行的数据，各个数据之间都用空格或者换行隔开 */
matrix(void * data, size_t r, size_t c=1, buffer * b=0); /* 数据构造函数，
行数为r,列数为c,data必须指向一存放实数的内存，按先列后行的
顺序进行存放。构造函数进行数据拷贝，因此在构造函数完成后,
data指向的数据区可以释放 */
size_t getrownum(){return istrans?colnum:rownum; };//获得矩阵的行数
size_t getcolnum(){return istrans?rownum:colnum; };//获得矩阵的列数
void loadfile(const char * filename);//从指定数据文件中获得矩阵数据
void savefile(char * filename,int g=2);//存矩阵数据到指定数据文件中
DOUBLE & operator()(size_t r, size_t c); // 重载函数运算符()返回第r行第c列的实数，
DOUBLE & operator()(size_t r){return (*this)(r,0);}; // 取回第r行0列的实数
virtual void set(size_t r, size_t c, DOUBLE v) {
size_t l;
v = isneg ? -v : v;
l = istrans ? r+c*rownum : r*colnum+c;
buf = buf->set(l, v); }; // 设置第r行第c列的值为v，
// 在设置时如果缓存的引用数大于1，产生新的缓存
// 因此每次设置返回的缓存区指针都要再赋给buf。
void set(size_t r, DOUBLE v) {
set(r,0,v);}; // 设置第r行0列的值为v,用作列向量
virtual ~matrix(){ // 析构函数，如果缓存的引用数大于1，说明还有其它的
// 矩阵在用此缓存，则不释放缓存，只是将引用数减1
buf->refnum--;
if(!buf->refnum) delete buf;};
matrix operator=(matrix m); // 重载赋值运算符，使矩阵的内容等于另一个矩阵
matrix operator=(DOUBLE a); // 通过赋值运算符将矩阵所有元素设为a
matrix operator-();
matrix operator*(DOUBLE a); // 重载乘法运算符，实现矩阵与数相乘，矩阵在前
// 数在后
matrix operator*=(DOUBLE a); // 重载*=运算符，功能同上，但*产生了一个新
// 的矩阵，原矩阵不变，而本运算符修改了原矩阵
friend matrix operator*(DOUBLE a, matrix& m);
// 数乘矩阵，但数在前矩阵在后
matrix operator+(matrix& m); // 矩阵相加，原二矩阵不变并产生一个新的和矩阵
matrix operator+=(matrix m); // 矩阵相加，不产生新矩阵，有一矩阵的内容改变为和

matrix operator+(DOUBLE a); // 矩阵加常数，指每一元素加一固定的常数，产生
// 新矩阵，原矩阵不变
friend matrix operator+(DOUBLE a, matrix& m); // 常数加矩阵，产生新的和矩阵
matrix operator+=(DOUBLE a); // 矩阵自身加常数，自身内容改变

matrix operator*(matrix m); // 矩阵相乘，产生出新的矩阵
matrix operator*=(matrix m); // 矩阵相乘，结果修改原矩阵

//matrix operator-(); // 矩阵求负，产生新的负矩阵
matrix trans(){
size_t r = rownum; rownum = colnum; colnum = r;
istrans = !istrans; return (*this);}; // 矩阵自身转置
matrix t(); // 矩阵转置，产生新的矩阵
matrix neg(){isneg = !isneg; return (*this);}; // 将自己求负

matrix operator-(matrix m); // 矩阵相减，产生新的矩阵
matrix operator-=(matrix m); // 矩阵相减，结果修改原矩阵
matrix operator-(DOUBLE a) // 矩阵减常数，指每一元素减一固定的常数，产生
// 新矩阵，原矩阵不变
{ return (*this)+(-a); };
friend matrix operator-(DOUBLE a, matrix m); // 常数减矩阵，产生新的矩阵
matrix operator-=(DOUBLE a) // 矩阵自身减常数，自身内容改变
{ return operator+=(-a);};

int isnear(matrix m, double e = defaulterr); // 检查两个矩阵是否近似相等,
// 如近似相等则返回1，否则返回0，e为允许误差
int isnearunit(double e = defaulterr); // 检查矩阵是否近似为单位矩阵
// 如是则返回1，否则返回0，e为允许误差

matrix row(size_t r); // 提取第r行行向量
matrix col(size_t c); // 提取第c列列向量

matrix mean(void);
matrix cov(void);

void checksym(); // 检查和尝试调整矩阵到对称矩阵

void swapr(size_t r1, size_t r2, size_t k=0); // 交换两行,只交换k列以上
void swapc(size_t c1, size_t c2, size_t k=0); // 交换两列,只交换k行以上
void swap(size_t r1, size_t c1, size_t r2, size_t c2){ // 交换两个元素
DOUBLE a; a=(*this)(r1,c1);set(r1,c1,(*this)(r2,c2));
set(r2,c2,a);};
DOUBLE maxabs(size_t &r, size_t &c, size_t k=0);
// 求主元值和位置，即k行k列之后的最大元素，
// 元素所在的行列将放在r,c中，返回此元素值
size_t zgsxy(matrix & m, int fn=0); // 主高斯消元运算
matrix operator/=(matrix m); // 用主高斯消元法求解线性方程的解
// 矩阵本身为系数矩阵，m为常数向量，返回解向量
matrix operator/(matrix m); // 左乘m的逆矩阵

matrix inv(); // 用全选主元高斯-约当法求逆矩阵
matrix operator~(); // 求逆矩阵，但产生新矩阵
friend matrix operator/(DOUBLE a, matrix m); // 求逆矩阵再乘常数
matrix operator/=(DOUBLE a); // 所有元素乘a的倒数，自身改变
matrix operator/(DOUBLE a); // 所有元素乘a的倒数，产生新的矩阵
matrix cholesky(matrix &d); // 用乔里斯基分解法求对称正定阵的线性
// 方程组ax=d返回方程组的解，本身为a，改变为分解阵u,d不变

DOUBLE det(DOUBLE err=defaulterr); // 求行列式的值
size_t rank(DOUBLE err=defaulterr); // 求矩阵的秩

void house(buffer & b, buffer & c); // 用豪斯荷尔德变换将对称阵变为对称三对
// 角阵，b返回主对角线元素，c返回次对角线元素
void trieigen(buffer& b, buffer& c, size_t l=600, DOUBLE eps=defaulterr);
// 计算三对角阵的全部特征值与特征向量
matrix eigen(matrix & eigenvalue, DOUBLE eps=defaulterr, size_t l=600);
// 计算对称阵的全部特征向量和特征值
// 返回按列排放的特征向量，而eigenvalue将返回一维矩阵，为所有特征值
// 组成的列向量
int ispositive(); // 判定矩阵是否对称非负定阵，如是返回1，否则返回0

void hessenberg(); // 将一般实矩阵约化为赫申伯格矩阵
void qreigen(matrix & a, matrix & b, size_t l=600, DOUBLE eps=defaulterr);
// 求一般实矩阵的所有特征根
// a和b均返回rownum行一列的列向量矩阵，返回所有特征根的实部和虚部

friend ostream& operator<<(ostream& o, matrix& m);
friend istream& operator>>(istream& in, matrix& m);
DOUBLE & value(size_t r, size_t c){ // 返回第r行c列的数值，不做行列检查
DOUBLE &a = istrans ? (*buf)[r+c*rownum] : (*buf)[r*colnum + c];
a=isneg ? -a : a;
return a;
};
protected:
DOUBLE detandrank(size_t & r, DOUBLE err); // 求方阵的行列式和秩
};

inline matrix operator*(DOUBLE a, matrix& m) {
return m*a;
};

inline matrix operator+(DOUBLE a, matrix& m) { // 常数加矩阵，产生新的和矩阵
return m+a;
};

matrix operator-(DOUBLE a, matrix& m); // 常数减矩阵，产生新的矩阵
matrix operator/(DOUBLE a, matrix& m); // 求逆矩阵再乘常数


ostream& operator<<(ostream& o, matrix& m);
istream& operator>>(istream& in, matrix& m);

matrix unit(size_t n); // 产生n阶单位矩阵

DOUBLE gassrand(int rr=0); // 返回一零均值单位方差的正态分布随机数
// rr为种子

class gassvector //返回零均值给定协方差阵的正态随机向量的类
{
public:
matrix a; // a是增益矩阵，由协方差矩阵算出
gassvector(matrix & r); //r必须是正定对称阵，为正态随机向量的协方差
matrix operator()(matrix & r); // 返回给定协方差矩阵的正态随机向量
matrix operator()(); // 返回已设定的协方差矩阵的正态随机向量
void setdata(matrix & r); // 根据协方差矩阵设置增益矩阵
};

DOUBLE lineopt(matrix& a,matrix& b, matrix& c, matrix & x);
// 线性规划最优点寻找程序，a为mXn不等式约束条件左端系数矩阵，b为不等式约束
// 条件的右端项，为m维向量，c为目标函数系数，n维向量，x返回极小点，为n维向量

#endif // MATRIX_H

文件5：buffer.cpp

#ifndef BUFFER_CPP
#define BUFFER_CPP
#include <stdio.h>
#include <string.h>
#include <fstream.h>
#include <math.h>

// 功能：实现buffer.h中的各个缓冲类的函数

#include "buffer.h"

DOUBLE defaulterr = 1e-8; // 缺省的误差调整值
int doadjust = 1; // 决定数据是否根据误差值自动向整数靠扰，
// 其实是向小数点后三位靠扰

void membuffer::Set(size_t n, DOUBLE d)
{
// 如调整标志设置，则调整
if(doadjust) adjust(d);
if(refnum == 1) { // 如果引用数为1，则将第n个实数的值设为d，并返回自己
// 的指针
retrieve(n) = d;
return this;
}
refnum --; // 引用数大于1，因此将本缓冲区的引用数减1
buffer * b;
b = clone(); // 克隆一个内容同自己一样的新缓冲区
(b->retrieve(n)) = d; // 将新缓冲区的第n个实数的值设为d
return b; // 返回新缓冲区的指针
}

membuffer* membuffer::clone() // 克隆一个内容和自己一样的实数内存缓冲区，并返回其指针
{
buffer* b;
b = new membuffer(length); // 建立一个新的实数内存缓冲区，尺寸和自己的一样
if(length > 0) // 如果自己有内容，将全部拷贝到新产生的缓冲区中
for(size_t i=0; i<length; i++)
(b->retrieve(i)) = (*this)[i];
return b;
}

void membuffer::alloc(size_t num){
DOUBLE *buf1;
buf1=new DOUBLE[num]; //分配新的内存
if(length>num)
for(int i=0;i<(int)num;i++) buf1[i]=buf[i];
else{
int i;
for(i=0;i<(int)length;i++) buf1[i]=buf[i];
for(i=length;i<(int)num;i++) buf1[i]=0;
}
delete []buf; //释放原来的内存
length = num;
buf=buf1; //把新的内存地址赋给 buf
};

// 实数磁盘缓冲区构造函数
diskbuffer::diskbuffer(size_t lth):buffer(),length(lth),n(0)
{
tempfp = tmpfile(); // 打开一新的临时文件，指针赋给tempfp
}

// 实数磁盘缓冲区的析构函数
diskbuffer::~diskbuffer()
{
fclose(tempfp); // 关闭临时文件
}

void diskbuffer::alloc(size_t num) // 将磁盘缓冲区重新定为尺寸为num
{
length = num;
if(!tempfp) // 如果临时文件尚未打开，则打开临时文件
tempfp = tmpfile();
n = 0; // 当前缓冲区指针设在开始，就是0处
}

void diskbuffer::release() // 释放磁盘缓冲区
{
fclose(tempfp); // 关闭临时文件
buf = 0.0;
length = 0;
}

DOUBLE& diskbuffer::retrieve(size_t i) // 返回实数磁盘缓冲区的第i个实数
{
long off;
off = n*sizeof(DOUBLE); // 计算出当前指针n对应的文件的偏移量
fseek(tempfp, off, SEEK_SET); // 移动文件指针到给定的偏移量处
fwrite(&buf, sizeof(DOUBLE), 1, tempfp); // 将当前的buf中值写到文件中
off = i*sizeof(DOUBLE); // 计算第i个实数在文件的什么位置
fseek(tempfp, off, SEEK_SET); // 移动文件指针到给定的偏移量处
fread(&buf, sizeof(DOUBLE), 1, tempfp); // 读回所需的实数值到buf中
n = i; // 并将n改为i
return buf; // 返回读出的实数的引用
}

buffer* diskbuffer::clone() // 克隆一个新的内容与自己一样的实数磁盘缓冲区
{
buffer* b;
b = new diskbuffer(length); // 定位产生一个与自己长度一样的新缓冲区
if(length > 0) // 如果长度大于0，则自己的内容拷贝到新缓冲区
for(size_t i=0; i<length; i++)
(b->retrieve(i)) = (*this)[i];
return b; // 返回新缓冲区的指针
}


// 将第n个长整数的值设为v，如果缓冲区的引用数大于1，
// 则会先克隆出新的缓冲区，然后再进行操作，并返回新缓冲区的指针
lbuffer* lbuffer::set(size_t n, long v)
{
if(refnum == 1) { // 引用数为1
retrieve(n) = v; // 将第n个长整数的值设为v
return this; // 返回自己
}
refnum --; // 引用数大于1，引用数减1
lbuffer * b;
b = clone(); // 克隆出新的内容一样的长整数缓冲区
(b->retrieve(n)) = v; // 将新缓冲区的第n个长整数值设为v
return b; // 返回新缓冲区的指针
}

ostream& operator<<(ostream& o, lbuffer& b)
{
size_t maxitem = 5;
for(size_t i=0; i<b.len(); i++) {
o << b[i] << '\t';
if(((i+1) % maxitem)==0)
o << endl;
}
o << endl;
return o;
}

lbuffer* lmembuffer::clone() // 克隆一个内容一新的长整数内存缓存区
{
lbuffer* b;
b = new lmembuffer(length); // 产生长度与自己一样的长整数缓存区
if(length > 0) // 如果长度大于0，则将自己的内容拷贝过去
for(size_t i=0; i<length; i++)
(b->retrieve(i)) = (*this)[i];
return b; // 返回新缓存区的指针
}

// 长整数磁盘缓存区的构造函数
ldiskbuffer::ldiskbuffer(size_t lth):lbuffer(),length(lth),n(0)
{
tempfp = tmpfile(); // 打开临时文件
}

// 长整数磁盘缓存区的析构函灵敏
ldiskbuffer::~ldiskbuffer()
{
fclose(tempfp); // 关闭临时文件
}

// 将长度定改为num
void ldiskbuffer::alloc(size_t num)
{
length = num;
if(!tempfp) // 如临时文件尚未打开，则打开
tempfp = tmpfile();
n = 0;
}

// 长整数磁盘缓存区释放
void ldiskbuffer::release()
{
fclose(tempfp); // 关闭临时文件
buf = 0;
length = 0;
}

// 检索第i个值的引用
long& ldiskbuffer::retrieve(size_t i)
{
long off;
off = n*sizeof(long); // 先将当前的内容保存
fseek(tempfp, off, SEEK_SET);
fwrite(&buf, sizeof(long), 1, tempfp);
off = i*sizeof(long); // 再取出第i个长整数到buf中
fseek(tempfp, off, SEEK_SET);
fread(&buf, sizeof(long), 1, tempfp);
n = i;
return buf;
}

// 克隆一个和内容和自己一样的长整数磁盘缓存区
lbuffer* ldiskbuffer::clone()
{
lbuffer* b;
b = new ldiskbuffer(length); // 产生长度和自己一样的缓存区
if(length > 0) // 如果自己内容不为空，拷贝内容到新的缓存区
for(size_t i=0; i<length; i++)
(b->retrieve(i)) = (*this)[i];
return b;
}

DOUBLE adjust(DOUBLE & a) // 将实数调整为最靠近小数点后二位的实数
// 准则：如果一个实数距某个二位小数小于十的负八次幂，则调整为这个小数
// 否则不变
{
DOUBLE b = floor(a*100.0+0.5)/100.0; // b是四舍五入到小数点后二位
if(fabs(b-a)<defaulterr) // 如果舍入误差极小，则抛弃
a = b;
return a;
}

char * throwmessage(int l, char * f, char * c)
{
static char a[100];
sprintf(a,"file:%s,line:%d\nmessage:%s\n",f,l,c);
return a;
}
#endif //BUFFER_CPP

文件6：buffer.h
#ifndef MEMBUFFER_H
#define MEMBUFFER_H

// 功能：定义数据缓冲区各抽象类

#ifndef DOUBLE

#define DOUBLE double
// #define DOUBLE long double // 如double精度不够，则可将此行去注释，而将上行注释
// 即可得到更精确的结果
#endif // DOUBLE

//#include <stdio.h>
#ifndef size_t
typedef unsigned int size_t
#endif

extern DOUBLE defaulterr; // 缺省的误差标准，为十的负八次方
extern int doadjust; // 逻辑变量，为1则作调整，否则不作，缺省为1
DOUBLE adjust(DOUBLE & a); // 将实数调整为最靠近的小数点后二位数

class membuffer{
public:
DOUBLE * buf; // 指向存放双精度实数数据的内存缓冲区的指针
size_t length; // 数组的以实数个数计的长度
membuffer():buf(0),length(0){}; // 构造函数，产生一个容量为0的实数缓冲区
membuffer(size_t lth):length(lth){ // 构造函数，产生长度为lth的实数缓冲区
if(lth > 0)
buf = new DOUBLE[lth];
else buf = 0; };
~membuffer(){delete []buf;} // 析构函数，释放缓冲区
void alloc(size_t num);
// 将缓冲区的容量改为num个实数，尽量保留原来的内容
DOUBLE& retrieve(size_t n){ // 返回第n个字节，n的范围从0到length-1
return buf[n]; };
void release() { delete []buf; length = 0; }; // 释放内存，数组的长度变为0
membuffer* clone(); // 克隆自己
membuffer* newbuf(size_t n=0){ // 产生一个长度为n的与自己同类别的新缓冲区并返回
// 相应的指针
return new membuffer(n);};
size_t len(){return length;}; // 返回数组的长度
void Set(size_t n,DOUBLE d);
};
#endif // BUFFER_H

[_xiaolifeidao]

wonderful!great thanks to wanbaocheng(www)　:)
let me try:

[_xiaolifeidao]

thanks!

[wanbaocheng]

不客气，有问题我们以后可继续讨论。