tsmatrix.cpp

程序简单实现了对稀疏矩阵进行incomplete Cholesky分解的功能

// TSMatrix.cpp: implementation of the TSMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "multiply.h"
#include "TSMatrix.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
//
//TSMatrix::TSMatrix()
//{
//
//}
//
//TSMatrix::~TSMatrix()
//{
//
// }

void TSMatrix::mulmatrix(TSMatrix &N, TSMatrix &Q)
{
    int mrow,nrow,tp,ccol;
	int* mnum = new int[mu+1];
	int* mpot = new int[mu+1];
	int* nnum = new int[N.mu+1];
	int* npot = new int[N.mu+1];
	float* a  = new float[17];
	TSMatrix TQ,MQ;
	TQ.data = new Triple[17];
	MQ.data = new Triple[16*N.nu+1];

    if(nu!=N.mu)
	{
		//MessageBox("error");
		// cerr<<"error\n"; exit(0);
	}

	Q.mu = mu;
	Q.nu = N.nu;
	Q.tu = 0;
	
	for (int i=1;i<=16;i++)
	{
		a[i] = 0;
	}

    for(mrow=1;mrow<=mu;++mrow)  
		   mnum[mrow] = 0;
	for(int t=1;t<=tu;++t)  
		   ++mnum[data[t].ii];
	mpot[1] = 1;
       //*(cpot+1)=1; //M中第一行第一个非零元在M中的序号为1
    for(mrow=2;mrow<=mu;++mrow)
		   mpot[mrow] = mpot[mrow-1] + mnum[mrow-1];

	for(nrow=1;nrow<=N.mu;++nrow)  
		   nnum[nrow] = 0;
	for(t=1;t<=N.tu;++t)  
		   ++nnum[N.data[t].ii];
	npot[1] = 1;                        //M中第一行第一个非零元在N中的序号为1
    for(nrow=2;nrow<=N.mu;++nrow)
		   npot[nrow] = npot[nrow-1] + nnum[nrow-1];
	
    
	if(tu!=0 && N.tu!=0)
	{
     for(mrow=1;mrow<=mu;++mrow)
	 {
		if(mrow<mu)
	      tp=mpot[mrow+1];
	    else
	      tp=tu+1;
		int tep = 0;
		TQ.tu = 0;
	    for(int p=mpot[mrow];p<tp;++p)
		{
			nrow=data[p].jj;
	        if(nrow<N.mu)
		        t=npot[nrow+1];
	        else
		        t=N.tu+1;
	        for(int q=npot[nrow];q<t;++q)
	           {
//				   ccol=N.data[q].jj;
//	               a[(mrow-1)*N.nu+ccol]+=data[p].e*N.data[q].e;
				   // Q.tu++;
				   TQ.tu++;
				   tep++;
				   a[tep] = N.data[q].jj;
				   TQ.data[TQ.tu].e  = data[p].e*N.data[q].e;
				   TQ.data[TQ.tu].ii = mrow;
				   TQ.data[TQ.tu].jj = a[tep];
	           }
		}

		int r,s,pt;
		Triple ft;
		for (r=1;r<tep;r++)
			for (s=1;s<=tep-r;s++)
			{
                if (a[s]>a[s+1])
                {
					ft = TQ.data[s];
					TQ.data[s] = TQ.data[s+1];
					TQ.data[s+1] = ft;
					pt = a[s];
					a[s] = a[s+1];
					a[s+1] = pt;

                }
				else 
					if (a[s]==a[s+1])
					{
						TQ.data[s].e+= TQ.data[s+1].e;
						for (int fc=s+1;fc<tep;fc++)
						{
							a[fc] = a[fc+1];
							TQ.data[fc] = TQ.data[fc+1];
						}
						tep--;
						TQ.tu--;
						s--;
					}
			}
			for (int tc=1;tc<=TQ.tu;tc++)
			{
				MQ.data[Q.tu+tc] = TQ.data[tc];
			}
			Q.tu+=TQ.tu;

//		for (ccol=1;ccol<=Q.nu;ccol++)
//		{
//			if (a[(mrow-1)*N.nu+ccol])
//			{
//				++Q.tu;
//
//				Q.data[Q.tu].e  = a[(mrow-1)*N.nu+ccol];
//				Q.data[Q.tu].ii = mrow;
//				Q.data[Q.tu].jj = ccol;
// 
//			}
//		}
	 }
	 Q.data = new Triple[Q.tu+1];
	 for (i=1;i<=Q.tu;i++)
	 {
         Q.data[i] = MQ.data[i];
	 }
	 delete MQ.data;
	}
}

void TSMatrix::FastTrMatrix(TSMatrix &R)
{
     int col,p,q,i,t;
  
   R.mu=nu;   R.nu=mu;   R.tu=tu;
   int* num  = new int[nu+1];
   int* cpot = new int[nu+1];
   R.data = new Triple[tu+1];
   if(tu)
   {
	   for(col=1;col<=nu;++col)  
		   num[col] = 0;
		   //*(num+col) = 0;
       //先置M每列非0元个数均为0
       for(t=1;t<=tu;++t)  
		   ++num[data[t].jj];
		   //++(*(num+T.data[t].jj));
       //求M中每一列非0元个数
	   cpot[1] = 1;
       //*(cpot+1)=1; //M中第一列第一个非元在T中的序号为1
       for(col=2;col<=nu;++col)
		   cpot[col] = cpot[col-1] + num[col-1];
          //*(cpot+col)=*(cpot+col-1) + *(num+col-1);
    //求M中第col列中第一个非0元在T中的序号
       for(p=1;p<=tu;++p)
	   {
		   col=data[p].jj; //记下M中第p个元素的列号
           q=cpot[col];     //该列中第一个非0元在T中的序号
           R.data[q].ii = data[p].jj;
           R.data[q].jj = data[p].ii;
           R.data[q].e  = data[p].e;//行、列交换,并赋值
           ++cpot[col];
	   } //该列下一个非0元序号=上一个序号+1
   }
		int errorr = 0;
	for (i=1;i<=tu-1;i++)
	{
		if ((R.data[i].ii==R.data[i+1].ii)&&(R.data[i].jj>=R.data[i+1].jj))
		{
			errorr =1;
		}
	}   
   
  
}

void TSMatrix::vectMulSMat(float *rp, float *result)
{
    int mrow;
	int i,tp,t;
	int* mnum = new int[mu+1];
	int* mpot = new int[mu+1];
    //result = new float[nu+1];
    for(mrow=1;mrow<=mu;++mrow)  
		   mnum[mrow] = 0;
	for( t=1;t<=tu;++t)  
		   ++mnum[data[t].ii];
	mpot[1] = 1;
       //*(cpot+1)=1; //M中第一行第一个非零元在M中的序号为1
    for(mrow=2;mrow<=mu;++mrow)
		   mpot[mrow] = mpot[mrow-1] + mnum[mrow-1];
    
	for (i=1;i<=mu;i++)
	{
		result[i]=0;
        if(i<mu)
	      tp=mpot[i+1];
	    else
	      tp=tu+1;
	    for(int p=mpot[i];p<tp;++p)
		{
			t = data[p].jj;
			result[i]+=data[p].e*rp[t];
		}
	}
	
}

void TSMatrix::addSMat(TSMatrix &N, TSMatrix &Q)
{
    Triple *Mp,*Me,*Np,*Ne,*Qh,*Qe;
	if (mu==N.mu && nu==N.nu)
	{
//		Q.mu = mu;
//		Q.nu = nu;
//		Q.data = new Triple[tu+N.tu+1];
		Mp = &data[1];         //Mp的初值指向矩阵M的非零元素首地址
		Np = &N.data[1];       //Np的初值指向矩阵N的非零元素首地址
		Me = &data[tu];        //Me的初值指向矩阵M的非零元素尾地址
		Ne = &N.data[N.tu];      //Ne的初值指向矩阵N的非零元素尾地址
		Qh = Qe = &Q.data[0];      //Qh,Qe的初值指向矩阵Q的非零元素首地址的前一地址
		while (Mp<=Me&&Np<=Ne)
		{
			Qe++;
			switch(comp(Mp->ii,Np->ii))
			{
			    case 1: *Qe = *Mp;
				        Mp++;
					    break;
				case 0: switch(comp(Mp->jj,Np->jj)) //M,N矩阵当前非零元素的行相等，继续比较列
						{
				            case 1: *Qe = *Mp;
								    Mp++;
									break;
							case 0: *Qe = *Mp;
								     Qe->e+=Np->e;
									if (!Qe->e)
									{
										Qe--;
									}
									Mp++;
									Np++;
									break;
							case -1:*Qe = *Np;
								    Np++;
									break;
									
						}
					    break;
				case -1: *Qe = *Np;
					     Np++;
						 break;
			}
		}
    if(Mp>Me)          //矩阵M的元素全部处理完毕
		while (Np<=Ne)
		{
			Qe++;
			*Qe = *Np;
			Np++;
		}
	if(Np>Ne)         //矩阵N的元素全部处理完毕
		while (Mp<=Me)
		{
			Qe++;
			*Qe = *Mp;
			Mp++;
		}
	Q.tu = Qe - Qh; 
	}
}

int TSMatrix::comp(int c1, int c2)
{
     int i;
	if(c1<c2)
		i = 1;
	else if(c1==c2)
		i = 0;
	else
		i = -1;
	return i;
}