[alt.sources] lin: solves systems of linear equations.

joe@dayton.UUCP (Joseph P. Larson) (01/04/90)

This is a really simple program for anyone about to take Linear Algebra
or perhaps who just wants to solve systems of linear equations.  To save
it, cut everything above the comment line plus the signature at the end,
save it to lin.c and type "make lin".

Caveat: if you're about to study linear algebra, I would *not* recommend
you use lin in lieu of learning the material.  Use it to help in your studies,
not replace them.  (Of course, you could check your problem sets over...)

Watch here in the future for more advance programs. -Joe

/*
 * lin: solve systems of linear equations.
 *
 * Description:
 *
 *		lin solves systems of linear equations.  Actually, it's
 *		not currently that smart -- all it does is solve very simple
 *		systems -- the kind with only one solution.  I whipped this
 *		while bored one day and will probably eventually make it a
 *		very smart little program.  But I thought, with a new term
 *		starting, there might be a few students out there about to
 *		take linear algebra who could use a tiny bit of help.  All
 *		it will do is try to mutilate the initial set of equations
 *		so that you can directly read off the values for x1, x2, etc.
 *
 *		A future version will solve systems with an infinite number
 *		of solutions.
 *
 *		See the code or the "usage" line for any help you need running it.
 *
 * Author:
 *
 *		Joe Larson
 *		(currently joe@dayton.dhdsc.mn.org)
 *
 * I hereby place this program in the public domain, as if it's big enough
 * for anyone to care, anyways. -Joe
 */

#include <stdio.h>
#include <ctype.h>

#define		MAXVAR	20
#define		MAXEQN	20

double strtod();
short	vcount, ecount;
double	matrix[MAXEQN][MAXVAR+1];
char	showstep = 0;

main(argc,argv)
int argc;
char **argv;
{
	char	*c,
			line[132],
			x;
	short	i,j;

	for (i = 1; i < argc; i++)
	{
		if (!strcmp(argv[i], "-s")) showstep = 1;
		else usage(argv[0]);
	}
	printf("This program solves systems of linear equations involving\n");
	printf("n variables and m unknowns.\n");
	while(1)
	{
		while (1)
		{
			printf("Number of variables: "); fflush(stdout);
			gets(line);
			vcount = atoi(line);
			if (vcount > MAXVAR)
			{
				printf("Max variables = %d.  ", MAXVAR);
				continue;
			}
			if (vcount) break;
			exit(0);
		}
		while (1)
		{
			printf("Number of equations: "); fflush(stdout);
			gets(line);
			ecount = atoi(line);
			if (ecount > MAXEQN)
			{
				printf("Max equations = %d.  ", MAXEQN);
				continue;
			}
			if (ecount) break;
			exit(0);
		}

		printf("Equations are of the form: \n");
		for (i = 0; i < vcount; i++)
		{
			if (i) printf("+ ");
			printf(" c%d * x%d ", i, i);
		}
		printf("= result\n");
		printf("For each equation, enter the %d constants plus the result,\n",
				vcount);
		printf("separated by spaces.\n");

		for (i = 0; i < ecount; i++)
		{
			printf("\nEquation %d\n", i+1);
			gets(line);
			for (c = line, j = 0; *c && j < ecount; j++)
				matrix[i][j] = strtod(c, &c);
			matrix[i][ecount] = strtod(c, &c);
		}

		printf("The system: \n");
		show();

	/*
	 * Now solve the equation in the following manner:
	 *
	 *		for each variable:
	 *				find an equation with a 1 for constant.
	 *					(if non, find one with a non-zero for constant
	 *					and multiply to make it a 1.)
	 *				add multiples of this equiation to others as needed
	 *					to zero-out their constants.
	 */

		for (i = 0; i < ecount; i++)
		{
			/* Find an equation to place here */
			x = 0;
			for (j = i; j < ecount; j++)
				if (matrix[j][i] == 1.0)
				{
					swap(i, j);
					x = 1;
					break;
				}
			if (!x)
				for (j = i; j < ecount; j++)
					if (matrix[j][i] != 0.0)
					{
						swap(i,j);
						x = 1;
						mult(i, 1.0 / matrix[i][i]);
						break;
					}
			if (!x)
			{
				/* This variable doesn't matter. */
				continue;
			}

			/*
			 * Now -- zero out everyone else.
			 */
			for (j = 0; j < ecount; j++)
			{
				if (j == i) continue;
				if (matrix[j][i] == 0.0) continue;
				add(j, i, -matrix[j][i]);
			}
		}

		printf("The system: \n");
		show();
	}
}

swap(i,j)
short i,j;
/*
 * Swap two rows of the matrix.
 */
{
	double	temp;
	short	x;

	if (i == j) return;
	for (x = 0; x <= vcount; x++)
	{
		temp = matrix[i][x];
		matrix[i][x] = matrix[j][x];
		matrix[j][x] = temp;
	}
	if (showstep)
	{
		printf("Swap row %d and %d leaving:\n", i, j);
		show();
	}
}

mult(i, f)
short i;
double f;
/*
 * Multiply row i by constant f.
 */
{
	short x;

	for (x = 0; x <= vcount; x++)
		matrix[i][x] *= f;
	if (showstep)
	{
		printf("Multiply row %d by %f leaving:\n", i, f);
		show();
	}
}

add(i,j,f)
short i,j;
double f;
/*
 * Add to row i row (j*f).
 */
{
	short x;

	for (x = 0; x <= vcount; x++)
		matrix[i][x] += matrix[j][x] * f;
	if (showstep)
	{
		printf("Add (%f times row %d) to row %d:\n", f, j, i);
		show();
	}
}

show()
/*
 * Show the equation
 */
{
	short i,j;

	for (i = 0; i < ecount; i++)
	{
		for (j = 0; j < vcount; j++)
			printf(" %8.4f", matrix[i][j]);
		printf("   %8.4f\n", matrix[i][j]);
	}
}

usage(arg)
char *arg;
{
	fprintf(stderr, "\
Usage: %d [-s]\n\
where -s selects display of steps used to solve the set of equations.\n", arg);
	
	exit(0);
}
-- 
UUCP: rutgers!dayton!joe   (Picts 1-16 are   DHDSC - Joe Larson/MIS 1060
ATT : (612) 375-3537       now ready.)       700 on the Mall, Mpls, Mn. 55402