/* luband3.c */

#include <math.h>
#include "luband3.h"

/* 帯行列の形をした連立方程式をGaussの消去法で解く。 */

/* LU分解 */
void bandlu(int m, int n, double *am[])
{
    int i, j, k, m1, mm;
    double aa;

    m1 = m - 1;
    for (i = 0; i < m1; i++) {
	mm = i + n + 1;
	if (mm > m)
	    mm = m;
	for (j = i + 1; j < mm; j++) {
	    aa = am[j][i - j + n] / am[i][n];	/* aa = a[j][i]/a[i][i]; */
	    for (k = i; k < mm; k++)
		am[j][k - j + n] -= aa * am[i][k - i + n];	/* a[j][k]-= aa * a[i][k]; */
	    am[j][i - j + n] = aa;
	}
    }
}

/* LU分解したものを使って方程式を解く */
void bandsol(int m, int n, double *am[], double f[])
{
    int i, j, jj, k, m1, mm;

    m1 = m - 1;
    for (i = 1; i < m; i++) {
	jj = i - n;
	if (jj < 0)
	    jj = 0;
	for (j = jj; j < i; j++)
	    f[i] -= f[j] * am[i][j - i + n];
    }

    f[m1] /= am[m1][n];
    for (i = m1 - 1; i >= 0; i--) {
	mm = i + n + 1;
	if (mm > m)
	    mm = m;
	for (j = i + 1; j < mm; j++)
	    f[i] -= f[j] * am[i][j - i + n];
	f[i] /= am[i][n];
    }
}