ccrosslist.cpp

一个用十字链表实现的稀疏矩阵类

				temp = 0;
				Element *x,*y;
				x=Rheaders[i];
				y=mx.CHeaders[j];
				while(x && y)
				{
					if(x->j == y->i)
					{
						temp+=x->value * y->value;
						x = x->right;
						y = y->down;
						
					}
					else if(x->j > y->i)
					{
						y = y ->down;
					}
					else
					{
						x = x->right;
					}
					
				}
				if(temp)
				{
					result.InsertElement(i,j,temp);
				}
			}
		}


	return result;
}


CrossList CrossList::operator *(double n)
{
	CrossList result(m_r,m_c);
	Element *x;
	for(int i=0; i<result.m_r; i++)
	{
		x = Rheaders[i];
		while(x)
		{
			result.InsertElement(x->i, x->j, x->value*n);
			x = x->right;
		}
	}
	return result;
}

CrossList CrossList::operator /(double n)
{
	CrossList result(m_r,m_c);
	Element *x;
	for(int i=0; i<result.m_r; i++)
	{
		x = Rheaders[i];
		while(x)
		{
			result.InsertElement(x->i, x->j, x->value/n);
			x = x->right;
		}
	}
	return result;
}

CrossList& CrossList::operator =(const CrossList& mx)
{
	if(&mx != this)
	{
		if(m_elements)
		{
			delete CHeaders;
			CHeaders=NULL;
			delete Rheaders;
			Rheaders=NULL;
			delete m_elements;
			m_elements = NULL;
		}
		m_r = mx.m_r;
		m_c = mx.m_c;
		m_elements_counter = 0;
		m_elements_MaxNum	=	mx.m_elements_MaxNum;///5;
		m_elements = new Element[m_elements_MaxNum];
		CHeaders = new Element*[m_c];
		Rheaders = new Element*[m_r];
		for(int i=0; i<m_r; i++)
		{
			Rheaders[i]=NULL;
		}
		for(i=0; i<m_c; i++)
		{
			CHeaders[i]=NULL;
		}
		Element *t;
		for(i=0; i<m_r; i++)
		{
			t = mx.Rheaders[i];
			while(t)
			{
				InsertElement(t->i, t->j, t->value);
				t=t->right;
			}
			
		}
	}
	return *this;
}


CrossList CrossList::Transpose() const
{
	CrossList result(m_c, m_r);
	Element *t;
	for(int i=0; i<m_r; i++)
	{
		t = Rheaders[i];
		while(t)
		{
			result.InsertElement(t->j, t->i, t->value);
			t=t->right;
		}
	}
	return result;
}

double CrossList::QtbyAbyQ(double* q)
{
	if(m_c!=m_r)
	{
		return -1;
	}
	double *temp = new double[this->m_c];
	double result = 0;
	Element *e;
	memset(temp,  0, sizeof(double)*(this->m_c));
	for(int i=0; i<this->m_c; i++)
	{
		e = CHeaders[i];
		while(e)
		{
			temp[i]+=(e->value) * (q[e->i]);
			e=e->down;
		}
	}
	for(i=0; i<m_r; i++)
	{
		result +=q[i]*temp[i];
	}
	delete temp;
	return result;
}

CrossList CrossList::SubAlphak(double alphak)
{
	CrossList result(m_r,m_c);
	Element *e;
	for(int i=0; i<result.m_r; i++)
	{
		e = Rheaders[i];
		while(e)
		{
			if(e->j == i)
			{
				result.InsertElement(i,i,e->value-alphak);
				
			}
			else
			{
				result.InsertElement(i,e->j,e->value);
			}
			e=e->right;
		}
	}
	return result;
}


void CrossList::Arraydivdouble(double* Array, double n, double* outArray, int sz)
{
	for (int i=0; i<sz; i++)
	{
		outArray[i] = Array[i]/n;
	}
}

void CrossList::RightMultiplyArray(double* inArray, double* out_Array)
{
	int sz = m_r;
	Element *e;
	memset(out_Array,0,sizeof(double)*sz);
	for (int i=0; i<m_r; i++)
	{
		e = Rheaders[i];
		while(e)
		{
			out_Array[i]+=e->value * inArray[e->j];
			e=e->right;
		}
	}
}

//void CrossList::Lanczos(CString path)
//{
//	ofstream ofile(path, ios::out | ios::ate);
//	double *a,*b;
//	int dim = 0;
//	if(m_r!=m_c)
//		return;
////	int n=this->m_r;
//	int n =5;
//	a=new double[n+1];
//	b=new double[n];
//	int k=0,i=0,p=0;
//	Element *e,*e0;
//
//	double* q = new double[n];
//	double* q1= new double[n];
//	double* q2= new double[n];
//	double* r = new double[n];
//	double* t1 = new double[n];
//	double* t2 = new double[n];
//	double* t3 = new double[n];
//
//	memset(q,0,sizeof(double)*n);
//	memset(q1,0,sizeof(double)*n);
//	memset(q2,0,sizeof(double)*n);
//	memset(r,0,sizeof(double)*n);
//	CrossList temp;
//	temp= *this;
//	b[0]=1;
//	a[0]=0;
//
//	r[0]=1;
//	q2[0]=1;
//
//
///////////////////////////////////////////////////////////////////////   begin
//	while(k<n)//floor(sqrt(n)))
//	{
//		
//		Arraydivdouble(r,b[k],q2,n);
//		k++;
//		a[k]=(QtbyAbyQ(q2));    //??????
//		for(i=0; i<m_r; i++)
//		{
//			e0=Rheaders[i];
//			e=temp.Rheaders[i];
//			while(e0 && e)
//			{
//				if(e->j == e->i)
//				{
//					e->value = e0->value;
//					break;
//				}
//				e = e->right;
//				e0= e0->right;
//			}
//		}
//		temp.SubAlphak1(a[k]);          ////????????
//		temp.RightMultiplyArray(q2,t1);    //t1=(*this-I)*q2   ??????
//		
//		for(p=0; p<n;p++)
//		{
//			t2[p]=q1[p]*b[k-1];
//		}
//		ArraySubArray(t1,t2,t3,n);
//		memcpy(r,t3,sizeof(double)*n);
//		b[k]=(ArrayNorm(r,n));
//		memcpy(q1,q2,sizeof(double)*n);
////		Q->iFace_SetColVector(k-1,n,q1);
//
//	}
//	ofile<<"alpha:"<<endl;
//	for (i=0; i<k+1; i++)
//	{
//		ofile<<"   "<<a[i];
//	}
//	ofile<<endl<<endl;
//	ofile<<"beta:"<<endl;
//	for (i=0; i<k; i++)
//	{
//		ofile<<"   "<<b[i];
//	}
//	ofile<<endl;
////////////////////////////////////////////////////////////////////////////    en
//	delete q;
//	delete q1;
//	delete q2;
//	delete r;
//	delete t1;
//	delete t2;
//	delete t3;
//}
void CrossList::ArraySubArray(double* inArray1, double* inArray2, double* outArray, int sz)
{
//	memset(outArray,0,sizeof(double)*sz);
	for (int i=0; i<sz; i++)
	{
		outArray[i] = inArray1[i] - inArray2[i];
	}
}

double CrossList::ArrayNorm(double* inArray, int n)
{
	double result=0;
	for (int i=0; i<n; i++)
	{
		result+=inArray[i]*inArray[i];
	}
	return sqrt(result);
}



//void CrossList::Lanczos(iFace_CMatrix *eigenvecs, int dim)
//{
//	double *a,*b;
//	if(m_r!=m_c)
//		return;
//	int n=this->m_r;
//	a=new double[n+1];
//	b=new double[n];
//	int k=0,i=0,p=0;
//	Element *e,*e0;
//	iFace_CMatrix* Q = new iFace_CMatrix(n,dim);
//	iFace_CMatrix* U = new iFace_CMatrix(dim,dim);
//	iFace_CMatrix* S = new iFace_CMatrix(dim,dim);
//	iFace_CMatrix* V = new iFace_CMatrix(dim,dim);
//	
//	double* q = new double[n];
//	double* q1= new double[n];
//	double* q2= new double[n];
//	double* r = new double[n];
//	double* t1 = new double[n];
//	double* t2 = new double[n];
//	double* t3 = new double[n];
//
//	memset(q,0,sizeof(double)*n);
//	memset(q1,0,sizeof(double)*n);
//	memset(q2,0,sizeof(double)*n);
//	memset(r,0,sizeof(double)*n);
//	CrossList temp;
//	CrossList subSparse(dim,dim);
//	temp= *this;
//	b[0]=1;
//	a[0]=0;
//
//	r[0]=1;
//	q2[0]=1;
//
///////////////////////////////////////////////////////////////////////   begin
//	while(k<dim)//floor(sqrt(n)))
//	{
//		
//		Arraydivdouble(r,b[k],q2,n);
//		k++;
//		a[k]=(QtbyAbyQ(q2));    //??????
//		for(i=0; i<m_r; i++)
//		{
//			e0=Rheaders[i];
//			e=temp.Rheaders[i];
//			while(e0 && e)
//			{
//				if(e->j == e->i)
//				{
//					e->value = e0->value;
//					break;
//				}
//				e = e->right;
//				e0= e0->right;
//			}
//		}
//		temp.SubAlphak1(a[k]);          ////????????
//		temp.RightMultiplyArray(q2,t1);    //t1=(*this-I)*q2   ??????
//		
//		for(p=0; p<n;p++)
//		{
//			t2[p]=q1[p]*b[k-1];
//		}
//		ArraySubArray(t1,t2,t3,n);
//		memcpy(r,t3,sizeof(double)*n);
//		b[k]=(ArrayNorm(r,n));
//		memcpy(q1,q2,sizeof(double)*n);
//		Q->iFace_SetColVector(k-1,n,q1);
//
//	}
//	for (i=1; i<=dim; i++)
//	{
//		if(i<dim)
//		{
//			subSparse.InsertElement(i-1,i-1,a[i]);
//			subSparse.InsertElement(i,i-1,b[i]);
//			subSparse.InsertElement(i-1,i,b[i]);
//		}
//		else
//		{
//			subSparse.InsertElement(i-1,i-1,a[i]);
//		}
//	}
//	subSparse.ToMatrix(S);
//	S->iFace_SplitUV(*U,*V);
//	*eigenvecs=(*Q)*(*U);
////	eigenvecs->iFace_SaveToFileDouble("C:\\evecs0.txt");
////////////////////////////////////////////////////////////////////////////    en
//
//
//	delete q;
//	delete q1;
//	delete q2;
//	delete r;
//	delete t1;
//	delete t2;
//	delete t3;
//
//
//	delete a;
//	delete b;
//	delete Q;
//	delete U;
//	delete S;
//	delete V;
//}

//void CrossList::ToMatrix(iFace_CMatrix* mx)
//{
//	int r = mx->iFace_GetNumRows();
//	int c = mx->iFace_GetNumColumns();
//	mx->iFace_Init(r,c);
//	Element *e;
//	for(int i=0; i<r; i++)
//	{
//		e=Rheaders[i];
//		while(e)
//		{
//			mx->iFace_SetElement(e->i, e->j, e->value);
//			e=e->right;
//		}
//	}
//}

void CrossList::SubAlphak1(double alphak)
{
	Element *e;
	for(int i=0; i<m_r; i++)
	{
		e = Rheaders[i];
		while(e)
		{
			if(e->j == i)
			{
				e->value-=alphak;
				
			}
			e=e->right;
		}
	}
}

void CrossList::CrossListDivDouble(double n)
{
	Element *e;
	for(int i=0; i<m_r; i++)
	{
		e = Rheaders[i];
		while(e)
		{
			e->value = e->value/n;
			e=e->right;
		}
	}
}

double CrossList::MaxNonnegtiveValue()
{
	double t = 0;
	Element *e;
	for(int i=0; i<m_r; i++)
	{
		e = Rheaders[i];
		while(e)
		{
			if(e->value > t)
			{
				t = e->value;
			}
			e=e->right;
		}
	}
	return t;
}