ilut.h

这是用于线性方程组求解的ILUT预处理算法的实现。在VC++编译通过。矩阵采用压缩稀疏行格式存储（CSR） 很容易移植到自己定义数值计算软件包中。经本人测试计算效率比Fortran写的高很多（比如与S

#include"malloc.h"
#include"math.h"
struct Distmatrix
{
double **ma;
int **ja;
int dimension;
int *nnzrow;
};
typedef Distmatrix FacLU;
struct Preconmat
{
	FacLU *L;
	FacLU *U;
};
struct Prepar
{
	int lfil,ufil;
	double droptol;
};
int quick_search(double *w,int*jw,int len_s,int len);
int ilut(struct Distmatrix *dm,struct Preconmat *precon,struct Prepar *ipar)
{
	int n,lfil,ufil,*rowj,*jw,lenl,lenu,len;
	int i,ii,j,jj,jpos,ierr,k,jcol,jrow,nnz;
	double *w,*rowm,tnorm,fact,s,t,droptol;
	n=dm->dimension;
	lfil=ipar->lfil;
	ufil=ipar->ufil;
	droptol=ipar->droptol;
	w=(double *)malloc(2*n*sizeof(double));
	jw=(int *)malloc(2*n*sizeof(int));
	for(j=0;j<n;j++)jw[n+j]=-1;
	ierr=0;
	precon->L=(FacLU *)malloc(sizeof(FacLU));
	precon->U=(FacLU *)malloc(sizeof(FacLU));
	precon->L->dimension=n;
	precon->U->dimension=n;
	precon->L->ja=(int**)malloc(n*sizeof(int *));
	precon->L->ma=(double **)malloc(n*sizeof(double *));
	precon->L->nnzrow=(int *)malloc(n*sizeof(int));
	precon->U->ja=(int**)malloc(n*sizeof(int *));
	precon->U->ma=(double **)malloc(n*sizeof(double *));
	precon->U->nnzrow=(int *)malloc(n*sizeof(int));
	for(ii=0;ii<n;ii++)
	{
		tnorm=0.0;
		for(k=0;k<dm->nnzrow[ii];k++)tnorm+=fabs(dm->ma[ii][k]);
		tnorm/=double(dm->nnzrow[ii]);
		lenu=1;
		lenl=0;
		jw[ii]=ii;
		w[ii]=0.0;
		jw[n+ii]=ii;
		for(j=0;j<dm->nnzrow[ii];j++)
		{
			jcol=dm->ja[ii][j];
			t=dm->ma[ii][j];
			if(jcol<ii)
			{
				jw[lenl]=jcol;
				jw[n+jcol]=lenl;
				w[lenl]=t;
				++lenl;
			}
			else if(jcol==ii)w[ii]=t;
			else
			{
				jpos=ii+lenu;
				jw[jpos]=jcol;
				w[jpos]=t;
				jw[n+jcol]=jpos;
				++lenu;
			}
		}
		len=0;
		/*eliminate previous rows*/
		for(jj=0;jj<lenl;jj++)
		{
			/*in order to eliminate the row ii in the correct order
			we must select the smallest column index among jw[k],
			k=jj+1,...,lenl.*/
			jrow=jw[jj];
			k=jj;
			for(j=jj+1;j<lenl;j++)
			{
				if(jw[j]<jrow)
				{
					jrow=jw[j];
					k=j;
				}
			}
			if(k!=jj)
			{
				/*exchange in jw*/
				j=jw[jj];
				jw[jj]=jw[k];
				jw[k]=j;
				/*exchange in jw[n+position]*/
				jw[n+jrow]=jj;
				jw[n+j]=k;
				/*exchange in w*/
				s=w[jj];
				w[jj]=w[k];
				w[k]=s;
			}
			jw[n+jrow]=-1;
			rowm=precon->U->ma[jrow];
			fact=w[jj]*rowm[0];
			if(fabs(fact)<droptol)continue;
			rowj=precon->U->ja[jrow];
			nnz=precon->U->nnzrow[jrow];
			/*combine the current row with the row jrow*/
			for(k=1;k<nnz;k++)
			{
				s=fact*rowm[k];
				jcol=rowj[k];
				jpos=jw[n+jcol];
				if(jcol>=ii)
				{
					if(jpos==-1)
					{
						/*this is a fill-in element*/
						i=ii+lenu;
						jw[i]=jcol;
						w[i]=-s;
						jw[n+jcol]=i;
						lenu++;
						if(lenu>n)
						{
							ierr=-1;
							break;
						}
					}
					else
					{
						w[jpos]-=s;
					}
				}
				else
				{
					if(jpos==-1)
					{
						/*this is a fill-in element*/
						jw[lenl]=jcol;
						w[lenl]=-s;
						jw[n+jcol]=lenl;
						++lenl;
						if(lenl>n)
						{
							ierr=-1;
							break;
						}
					}
					else w[jpos]-=s;
				}
			}
			w[len]=fact;
			jw[len]=jrow;
			len++;
		}
		/*reset jw[n+pointer] to -1---only U-part necessarily*/
		for(k=0;k<lenu;k++) jw[n+jw[ii+k]]=-1;
		/*update the row ii of the Lpart*/
		lenl=len>lfil?lfil:len;
		precon->L->nnzrow[ii]=lenl;
		/*weight the elements before sort*/
		for(k=0;k<len;k++)w[k]*=(droptol+w[n+jw[k]]);/*questioning?*/
		if(len>lenl)quick_search(w,jw,lenl,len);
		if(lenl>0)
		{
			precon->L->ja[ii]=(int *)malloc(lenl*sizeof(int));
			precon->L->ma[ii]=(double *)malloc(lenl*sizeof(double));
		}
		rowj=precon->L->ja[ii];
		rowm=precon->L->ma[ii];
		for(k=0;k<lenl;k++)rowj[k]=jw[k];
		for(k=0;k<lenl;k++)rowm[k]=w[k]/(droptol+w[n+jw[k]]);
		/*update the row ii of the U-part*/
		len=0;
		for(k=1;k<lenu;k++)
		{
			if(fabs(w[ii+k])>droptol*fabs(w[ii]))
			{
				++len;
				w[ii+len]=w[ii+k];
				jw[ii+len]=jw[ii+k];
			}
		}
		lenu=len+1>ufil?ufil:len+1;
		jpos=lenu-1;
		if(len>jpos)quick_search(&w[ii+1],&jw[ii+1],jpos,len);
		precon->U->nnzrow[ii]=lenu;
		precon->U->ja[ii]=(int *)malloc(lenu*sizeof(int));
		precon->U->ma[ii]=(double *)malloc(lenu*sizeof(double));
		rowj=precon->U->ja[ii];
		rowm=precon->U->ma[ii];
		for(k=0;k<lenu;k++)
		{
			rowj[k]=jw[ii+k];
			rowm[k]=w[ii+k];
		}
		t=fabs(w[ii]);
		for(k=1;k<lenu;k++)t+=fabs(w[ii+k]);
		w[n+ii]=t/(double)(lenu+1);
		if(w[ii]==0.0)w[ii]=(0.0001+droptol)*tnorm;
		rowm[0]=1.0/w[ii];
		rowj[0]=ii;		
	}
	return 0;
	free(w);
	free(jw);
	return ierr;
}
int quick_search(double *w,int*jw,int len_s,int len)
{
	double temp,abskey;
	int first,last,mid;
	int j,itemp;
	first=0;
	last=len-1;
	if(len_s<first||len_s>=last)return 0;
label1:
	mid=first;
	abskey=fabs(w[mid]);
	for(j=first+1;j<=last;j++)
	{
		if(fabs(w[j])>abskey)
		{
			++mid;
			temp=w[mid];
			w[mid]=w[j];
			w[j]=w[mid];
			itemp=jw[mid];
			jw[mid]=jw[j];
			jw[j]=itemp;
		}
	}
	/*interchange w[first] with w[mid]*/
	temp=w[first];
	w[first]=w[mid];
	w[mid]=temp;
	itemp=jw[first];
	jw[first]=jw[mid];
	jw[mid]=itemp;
	if(mid==len_s) return 0;
	if(mid>len_s)last=mid-1;
	else
		first=mid+1;
	goto label1;
}