iluk.h

这是用于线性方程组求解的ILUK预处理算法的实现。在VC++编译通过。矩阵采用压缩稀疏行格式存储（CSR）

#include "malloc.h"
#define min(a, b)  a < b ? a : b
struct Distmatrix
{
double **ma;
int    **ja,dimension,*nnzrow;
};
typedef Distmatrix FacLU;
struct Preconmat
{
	FacLU *L;
	FacLU *U;
};
int iluk(struct Distmatrix *dm,struct Preconmat *precon,int fil_level)
{
	int *jw,jcol,jrow,lenl,lenu,*lev,nnz,*rowj;
	int **levs;
	int jlev,i,j,ii,jj,jpos,k,ierr=0;
	double *w,*rowm,t,s,fact;
	int n,n2;
	n=dm->dimension;
	if(fil_level<0)
	{
		ierr=-2;
		return ierr;
	}
	levs=(int **)malloc(n*sizeof(int*));
	w=(double *)malloc(n*sizeof(double));
	jw=(int *)malloc(3*n*sizeof(int));
	precon->L->ma=(double **)malloc(n*sizeof(double *));
	precon->L->ja=(int **)malloc(n*sizeof(int *));
	precon->L->nnzrow=(int *)malloc(n*sizeof(int));
	precon->U->ma=(double **)malloc(n*sizeof(double *));
	precon->U->ja=(int **)malloc(n*sizeof(int *));
	precon->U->nnzrow=(int *)malloc(n*sizeof(int));
	n2=n<<1;
	for(i=n;i<n2;i++)
	{
		jw[i]=-1;
		jw[n+i]=0;
	}
	for(ii=0;ii<n;ii++)
	{
		lenl=0;
		lenu=1;
		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(t==0.0)continue;
			if(jcol<ii)
			{
				jw[lenl]=jcol;
				w[lenl]=t;
				jw[n+jcol]=lenl;
				jw[n2+lenl]=0;
				++lenl;
			}
			else if(jcol==ii)
			{
				w[ii]=t;
				jw[n2+ii]=0;
			}
			else
			{
				jpos=ii+lenu;
				jw[jpos]=jcol;
				w[jpos]=t;
				jw[n+jcol]=jpos;
				jw[n2+jpos]=0;
				++lenu;
			}
		}
		/*taking the linear combination of the previous rows with row ii*/
		for(jj=0;jj<lenl;jj++)
		{
			/*in order to eliminate in the correct order, we must select
			the smallest column index among jw[k],k=jj,...,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+(points/nonzero indicator)]*/
				jw[n+jrow]=jj;
				jw[n+j]=k;
				/*exchange in jw[n2+(levels)]*/
				j=jw[n2+jj];
				jw[n2+jj]=jw[n2+k];
				jw[n2+k]=j;
				/*exchange in w*/
				s=w[jj];
				w[jj]=w[k];
				w[k]=s;
			}
			jw[n+jrow]=-1;
			rowj=precon->U->ja[jrow];
			rowm=precon->U->ma[jrow];
			nnz=precon->U->nnzrow[jrow];
			lev=levs[jrow];
			fact=w[jj]*rowm[0];
			jlev=jw[n2+jj];
			if(jlev>fil_level)continue;
			/*combinate the current row ii with 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;
						jw[n+jcol]=i;
						w[i]=-s;
						jw[n2+i]=jlev+lev[k]+1;
						++lenu;
						if(lenu>n)
						{
							ierr=-1;
							break;
						}
					}
					else
					{
						w[jpos]-=s;
						jw[n2+jpos]=min(jw[n2+jpos],lev[k]+jlev+1);
					}									
				}
				else
				{
					if(jpos==-1)
					{
						/*This is a fill-in element*/
						jw[lenl]=jcol;
						jw[n+lenl]=lenl;
						w[lenl]=-s;
						jw[n2+lenl]=jlev+lev[k]+1;
						++lenl;
						if(lenl>n)
						{
							ierr=-1;
							break;
						}
					}
					else
					{
						/*this is not a fill-in element*/
						w[jpos]-=s;
						jw[n2+jpos]=min(jw[n2+jpos],lev[k]+jlev+1);
					}
				}
			}
			w[jj]=fact;
			jw[jj]=jrow;
		}
		for(k=0;k<lenu;k++)jw[n+jw[ii+k]]=-1;
		j=0;
		for(k=0;k<lenl;k++)
		{
			if(jw[n2+k]<=fil_level)j++;
		}
		precon->L->nnzrow[ii]=j;
		if(j>0)
		{
			precon->L->ma[ii]=(double *)malloc(j*sizeof(double));
			precon->L->ja[ii]=(int *)malloc(j*sizeof(int));
		}
		rowj=precon->L->ja[ii];
		rowm=precon->L->ma[ii];
		j=0;
		for(k=0;k<lenl;k++)
		{
			if(jw[n2+k]<=fil_level)
			{
				rowj[j]=jw[k];
				rowm[j]=w[k];
				++j;
			}
		}
		if(w[ii]==0.0)
		{
			ierr=-3;
			break;
		}
		j=1;
		for(k=ii+1;k<ii+lenu;k++)
		{
			if(jw[n2+k]<=fil_level)j++;
		}
		precon->U->nnzrow[ii]=j;
		precon->U->ma[ii]=(double *)malloc(j*sizeof(double));
		precon->U->ja[ii]=(int *)malloc(j*sizeof(int));
		rowj=precon->U->ja[ii];
		rowm=precon->U->ma[ii];
		levs[ii]=(int *)malloc(j*sizeof(int));
		lev=levs[ii];
		rowm[0]=1.0/w[ii];
		rowj[0]=ii;
		j=1;
		for(k=ii+1;k<ii+lenu;k++)
		{
			if(jw[n2+k]<=fil_level)
			{
				rowj[j]=jw[k];
				rowm[j]=w[k];
				lev[j]=jw[n2+k];
				j++;
			}
		}		
	}
	for(i=0;i<ii;i++) free(levs[i]);
	free(levs);
	free(w);
	free(jw);
	return ierr;
}