band.cpp

采用压缩存储方法

#include "stdio.h"
#include "math.h"
#include "malloc.h"
#include "stdlib.h"
#include "string.h"

int band(int n,int h,int ih,int m,double b[],double d[])
{ 
  int ls,k,i,j,is,u,v;
  double p,t;
  if(ih!=(2*h+1))
  {
    printf("fail\n");
	return(-2);;   /*带宽与半带宽不匹配，错误返回*/
  }

  ls=h;
  for(k=0;k<=n-2;k++)
  {
    p=0.0;
	for(i=k;i<=ls;i++)     /*列选主元*/
	{
	  t=fabs(b[i*ih]);
	  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++)   /*m组常数向量行交换*/
	{
	  u=k*m+j;
	  v=is*m+j;
	  t=d[u];
	  d[u]=d[v];
	  d[v]=t;
	}

	for(j=0;j<=ih-1;j++)   /*系数矩阵行交换*/
	{
	  u=k*ih+j;
	  v=is*ih+j;
	  t=b[u];
	  b[u]=b[v];
	  b[v]=t;
	}

	for(j=0;j<=m-1;j++)   /*m组常数向量归一化*/
	{
	  u=k*m+j;
	  d[u]=d[u]/b[k*ih];
	}

	for(j=1;j<=ih-1;j++)  /*系数矩阵归一化*/
	{
	  u=k*ih+j;
	  b[u]=b[u]/b[k*ih];
	}

	for(i=k+1;i<=ls;i++)
	{
	  t=b[i*ih];
	  for(j=0;j<=m-1;j++)   /*m组常数向量消元*/
	  {
	    u=i*m+j;
		v=k*m+j;
		d[u]=d[u]-t*d[v];
	  }
	  for(j=1;j<=ih-1;j++)  /*系数矩阵消元，且左移一列*/
	  {
	    u=i*ih+j;
		v=k*ih+j;
		b[u-1]=b[u]-t*b[v];
	  }
	  u=i*ih+ih-1;
	  b[u]=0.0;
	}
	if(ls!=(n-1))ls=ls+1;
  }

  p=b[(n-1)*ih];
  if(fabs(p)+1.0==1.0)
  {
    printf("fail\n");  /*方程组系数奇异，错误返回*/
	return(0);  
  }

  for(j=0;j<=m-1;j++)  /*计算m组解向量中的最后一个分量*/
  {
    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*ih+j;
		is=(i+j)*m+k;
		d[u]=d[u]-b[v]*d[is];
	  }
	}
	if(ls!=(ih-1))ls=ls+1;
  }
  return(2);      /*正常返回*/
}



main()
{
  int i;
  static double b[8][5]={{3.0,-4.0,1.0,0.0,0.0},{-2.0,-5.0,6.0,1.0,0.0},{1.0,3.0,-1.0,2.0,-3.0},{2.0,5.0,-5.0,6.0,-1.0},{-3.0,1.0,-1.0,2.0,-5.0},{6.0,1.0,-3.0,2.0,-9.0},{-4.0,1.0,-1.0,2.0,0.0},{5.0,1.0,-7.0,0.0,0.0}};
  static double d[8][3]={{13.0,29.0,-13.0},{-6.0,17.0,-21.0},{-31.0,-6.0,4.0},{64.0,3.0,16.0},{-20.0,1.0,-5.0},{-22.0,-41.0,56.0},{-29.0,10.0,-21.0},{7.0,-24.0,20.0}};

  double *temp1 = &b[0][0];
  double *temp2 = &d[0][0];
  
  if(band(8,2,5,3,temp1,temp2)>0)
	  for(i=0;i<=7;i++)
		  printf("x(%d)=%15.7e,    %15.7e,     %15.7e\n",i,d[i][0],d[i][1],d[i][2]);
}