solve.c

an introduction to boundary element methods一书源码

#include "cbox1.h"

/*[function Solve for the solution of a linear system of equations by the Gauss Elimination method providing for interchange of rows when a zero diagonal coeficient is encountered.
A : system matrix
B:  Originally it contains the Independent coefficients.  After solutions it
	contains the values of the system unknowns.                            				   				
N:  Actual number of unknowns.
Dim: Row and Column Dimension of A.  ]*/

void Solve(A,B,D,N,Dim)
	float A[100][100],B[101],*D;  /*[ D is assigned some value here ]*/
  	int N,Dim;
{
	int N1,i,j,k,i1,l,k1,found;
	float c;
/*[ found is a flag which is used to check if any nonzero coeff is found ]*/	
	N1 = N - 1;
	for(k=1;k<=N1;k++)
	{
	k1 = k + 1;
	c = A[k][k];
	if( (fabs(c) - tol) <=0 )
	{
	found = 0;
	for(j=k1;j<=N;j++)  /*[To interchange rows to get nonzero coeff  ]*/
 	{
         if( (fabs(A[j][k]) - tol) > 0.0 ) 
           {
		for(l=k;l<=N;l++)
                    {
                     	c = A[k][l];
                        A[k][l] = A[j][l];
                        A[j][l] = c;
                     }
		c = B[k];
		B[k] = B[j];
		B[j] = c;
		c = A[k][k];
		found = 1;		/*[ coeff is found ]*/
		break;
	      }
	 }
	}

		if(!found)
		{
		printf("Singularity in Row %d 1",k);
	   	(*D) = 0.0;          
	    	return;   /*[ If no coeff is found the control is transferred to main ]*/
		}
	/*[ Divide Row by Diagonal coefficient ]*/		
  		c = A[k][k];
		for(j = k1;j<=N;j++)
		   A[k][j] /= c;
		   B[k] /= c;  
/*[ Eliminate unknown X[k] from Row i ]*/
		for(i = k1;i<=N; i++)
		{
		  c = A[i][k];
		  for(j = k1;j<= N;j++)
		 	A[i][j] -= c*A[k][j];
		        B[i] -= c*B[k];
		}
	}
/*[ Compute the last unknown ]*/
	if((fabs(A[N][N]) - tol) > 0.0)
	{
	   B[N] /= A[N][N];
/*[Apply back substitution to compute the remaining unknowns  ]*/
	   for(l = 1;l<= N1;l++)
	   {
	      k = N - l;
	      k1 = k +1;
	      for(j = k1;j<=N;j++)
			B[k] -= A[k][j]*B[j];
	   }
/*[Compute the value of the determinent ]*/
	(*D) = 1.0;
	for(i=1;i<=N;i++)
	  (*D) *= A[i][i];

   }
  else
  {
  		printf("Singularity in Row %d 2",k);
	   	(*D) = 0.0;
   }          
  return;	
}