連立１次方程式ライブラリ 使用例

/*		test13.c				*/
#include <stdio.h>
#include "sslib.h"

main()
{
	int error, i, iter, j, k, l, m, n;

	static double a1[12] = { 5, 1, 1,10, 1,-7, 2,-7, 1, 1, 6,21};
	static double a2[15] = { 5, 1, 1,10,18, 1,-7, 2,-7,-9, 1, 1, 6,21,11};
	static double ar[12] = { 2, 1, 3,10, 1,-2, 5,-3,-4, 5, 1, 2};
	static double ai[12] = {-1, 2, 1, 2, 1,-4, 3,-1,-1, 1, 4, 4};
	Complex a[12], true[3], cerr;
	double ax[12], aw[15], x[3], err[3], eps;
	double al[6], am[6], au[6], b[6], xx[6];

	printf("*** Gauss-Seidel method ***\n");
	m = 3;
	l = m + 1;
	iter = 20;
	eps = 1.e-6;
	for(i = k = 0; i < m; i++)	for(j = 0; j < l; j++, k++)	ax[k] = a1[k];
	for(i = k = 0; i < m; i++)
	{
		printf("    [ ");
		for(j = 0; j < m; j++)	printf("%5.1f  ", ax[k++]);
		printf("] [ x%1d ] = [ %5.1f ]\n", i + 1, ax[k++]);
	}
	gausei(ax, l, m, iter, eps, x);
	printf("\n* test result from gausei()\n");
	err[0] = x[0] - 1.0;
	err[1] = x[1] - 2.0;
	err[2] = x[2] - 3.0;
	for(i = 0; i < m; i++)	printf(" x%1d = %22.15e  (err=%e)\n", i + 1, x[i], err[i]);

	printf("\n\n*** Gauss method ***\n");
	m = 3;
	l = n = 5;
	eps = 1.e-6;
	for(i = k = 0; i < m; i++)	for(j = 0; j < n; j++, k++)	aw[k] = a2[k];
	for(i = k = 0; i < m; i++)
	{
		printf("   [ ");
		for(j = 0; j < m; j++)	printf("%5.1f ", aw[k++]);
		printf("] [ x%1d ] = ", i + 1);
		for(j = m; j < l; j++)	printf("[ %5.1f ] ", aw[k++]);
		putchar('\n');
	}
	gauss(aw, l, m, n, eps);
	printf("\n* test result from gauss()\n");
	for(i = k = 0; i < m; i++)
	{
		printf(" x%1d = ", i + 1);
		for(j = m, k += m; j < l; j++)	printf("%22.15e  ", aw[k++]);
		putchar('\n');
	}

	printf("\n\n*** Gauss-Jordan method ***\n");
	m = 3;
	l = n = 5;
	eps = 1.e-6;
	for(i = k = 0; i < m; i++)	for (j = 0; j < l; j++, k++)	aw[k] = a2[k];
	for(i = k = 0; i < m; i++)
	{
		printf("   [ ");
		for(j = 0; j < m; j++)	printf("%5.1f ", aw[k++]);
		printf("] [ x%1d ] = ", i + 1);
		for(j = m; j < l; j++)	printf("[ %5.1f ] ", aw[k++]);
		putchar('\n');
	}
	gaujor(aw, l, m, n, eps);
	printf("\n* test result from gaujor()\n");
	for(i = k = 0; i < m; i++)
	{
		printf(" x%1d = ", i + 1);
		for(j = m, k += m; j < l; j++)	printf("%22.15e  ", aw[k++]);
		
		putchar('\n');
	}

	printf("\n\n*** Solution linear equation of tridiagonal matrix\n");
	n = 5;
	for(i = 0; i <= n; i++)
	{
		al[i] = i + 1;
		au[i] = i + 2;
		am[i] = (al[i] + au[i]) * 2.0;
		b[i] = al[i] + am[i] + au[i];
	}
	am[0] += al[0];
	am[n] += au[n];
	printf("  Coefficient   A\n");
	for(i = 0; i <= n; i++)
	{
		if(i == 0)
			printf("                 am[%2d]=%4.1f   au[%2d]=%4.1f    b[%2d]=%4.1f\n",
				i, am[i], i, au[i], i, b[i]);
		if(i == n)
			printf("   al[%2d]=%4.1f   am[%2d]=%4.1f                  b[%2d]=%4.1f\n",
				i, al[i], i, am[i], i, b[i]);
		if((i != 0) & (i != n))
			printf("   al[%2d]=%4.1f   am[%2d]=%4.1f   au[%2d]=%4.1f    b[%2d]=%4.1f\n",
				i, al[i], i, am[i], i, au[i], i, b[i]);
	}
	error = trdiam(al, am, au, b, n, xx);
	for(i = 0; i <= n; i++)
	{
		if(i % 3 == 0)	printf("\n");
		printf("  x[%2d]= %7.5f  ", i, xx[i]);
	}
	printf("\n   Error  = %5d\n", error);

	printf("\n\n*** Gauss-Jordan method (Complex) ***\n");
	l = n = 4;
	m = 3;
	eps = 1.e-6;
	for(i = k = 0; i < m; i++)	for (j = 0; j < l; j++, k++)	a[k] = tocomplex(ar[k], ai[k]);
	for(i = k = 0; i < m; i++)
	{
		printf(" [ ");
		for(j = 0; j < m; j++, k++)	printf("(%5.1f,%5.1f)  ", a[k].r, a[k].r);
		printf("] [ x%1d ] = ", i + 1);
		for(j = m; j < l; j++, k++)	printf("[ (%5.1f,%5.1f) ] ", a[k].r, a[k].i);
		putchar('\n');
	}
	cgauj(a, l, m, n, eps);
	true[0] = tocomplex(1876238. / 853450., -1020591. / 853450);
	true[1] = tocomplex(2112162. / 853450., -914459. / 853450);
	true[2] = tocomplex(60762. / 65650., 39091. / 65650.);
	printf("\n* test result from cgauj()\n");
	for(i = k = 0; i < m; i++)
	{
		printf(" x%1d = ", i + 1);
		for(j = m, k += m; j < l; j++, k++)
		{
			printf("(%10.7f,%10.7f)    ", a[k].r, a[k].i);
			cerr = csub(a[k], true[i]);
			printf("err=(%e,%e)   ", cerr.r, cerr.i);
		}
		putchar('\n');
	}

	return 1;
}