6.6 求解一般带型方程组 band.c

许士良常用算法程序集C语言,包括c++一些常用算法代码

#include "math.h"
#include "stdio.h"
int band(b,d,n,l,il,m)
int n,l,il,m;
double b[],d[];
{ 
	int ls,k,i,j,is,u,v;
    double p,t;
    if (il!=(2*l+1))
    { 
		printf("fail\n"); 
		return(-2);
	}
    ls=l;
    for (k=0;k<=n-2;k++)
    { 
		p=0.0;
        for (i=k;i<=ls;i++)
        { 
			t=fabs(b[i*il]);
            if (t>p) {p=t; is=i;}
        }
        if (p+1.0==1.0)
        { 
			printf("fail\n"); 
			return(0);
		}
        for (j=0;j<=m-1;j++)
        { 
			u=k*m+j; v=is*m+j;
            t=d[u]; d[u]=d[v]; d[v]=t;
        }
        for (j=0;j<=il-1;j++)
        { 
			u=k*il+j; v=is*il+j;
            t=b[u]; b[u]=b[v]; b[v]=t;
        }
        for (j=0;j<=m-1;j++)
        { 
			u=k*m+j; d[u]=d[u]/b[k*il];
		}
        for (j=1;j<=il-1;j++)
        { 
			u=k*il+j; b[u]=b[u]/b[k*il];
		}
        for (i=k+1;i<=ls;i++)
        { 
			t=b[i*il];
            for (j=0;j<=m-1;j++)
            { 
				u=i*m+j; v=k*m+j;
                d[u]=d[u]-t*d[v];
            }
            for (j=1;j<=il-1;j++)
            { 
				u=i*il+j; v=k*il+j;
                b[u-1]=b[u]-t*b[v];
            }
            u=i*il+il-1; b[u]=0.0;
        }
        if (ls!=(n-1)) ls=ls+1;
    }
    p=b[(n-1)*il];
    if (fabs(p)+1.0==1.0)
    { 
		printf("fail\n"); 
		return(0);
	}
    for (j=0;j<=m-1;j++)
    { 
		u=(n-1)*m+j; d[u]=d[u]/p;
	}
    ls=1;
    for (i=n-2;i>=0;i--)
    { 
		for (k=0;k<=m-1;k++)
        { 
			u=i*m+k;
            for (j=1;j<=ls;j++)
            { 
				v=i*il+j; is=(i+j)*m+k;
                d[u]=d[u]-b[v]*d[is];
            }
        }
        if (ls!=(il-1)) ls=ls+1;
    }
    return(2);
}