solve.c

an introduction to boundary element methods一书源码

#include "cbox5.h"


/*[ function  Solve : Solution of the linear systems of equations by the Gauss elimination method providing for interchanging rows when encountering a zero diagonal coeficient.

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[101][101],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 non zero 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++)  /*[Interchange rows to get Nonzero ]*/
	
	        {
            	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 \n",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 \n",k);
	   	(*D) = 0.0;
   }          

  return;	
}