/* 
 * PSC98 HONSEN input data generation
 */

#include 
#include 
#include 
#include 

#define MATMAX	2000000
#define NZMAX	9

/* input data generatoin */

static	int	key = -1;

static int honsen1();
static int honsen2();
static int honsen3();
static int honsen4();
static int honsen5();
static int sample();

static void genmat_common(i, ja, a, b)
    int i, *ja;
    double *a, *b;
{
	char	*buff;

	if ( key == -1 ) {
		buff = getenv( "PSC98" );
		if ( buff == NULL ) buff = "0";
		key = atoi(buff);
	}

	switch (key) {
		case 1: honsen1(i, ja, a, b); break;
		case 2: honsen2(i, ja, a, b); break;
		case 3: honsen3(i, ja, a, b); break;
		case 4: honsen4(i, ja, a, b); break;
		case 5: honsen5(i, ja, a, b); break;
		default: key=0; sample(i, ja, a, b); break;
	}
}

/*
 *  Fortran77 Interface Routine
 */
/*
 *  call genmat(i, ja, a, b)
 */
void
genmat_(i, ja, a, b)
    int    *i, ja[];
    double a[], *b;
{
#if IDX_BEGINS_WITH_ZERO
    genmat_common(*i, ja, a, b);
#else
    int i_prime;
    if (*i == -1) i_prime = -1;
    else i_prime = *i - 1;
    genmat_common(i_prime, ja, a, b);
#endif
}

/*
 *  C Interface Routine
 */

void genmat(i, ja, a, b)
    int i, *ja;
    double *a, *b;
{
#if IDX_BEGINS_WITH_ZERO
    genmat_common(i, ja, a, b);
#else
    int i_prime;
    if (i == -1) i_prime = -1;
    else i_prime = i - 1;
    genmat_common(i_prime, ja, a, b);
#endif
}


void
chkval(s, n, x)
    FILE * s;			/* output stream */
    int n;
    double *x;
{
    double a[NZMAX];
    int ja[NZMAX], i, mat_size;
    double b, resmax, max_r = 0.0;

    genmat(-1, ja, a, &b);
    mat_size = ja[0];
    resmax = b;

    assert(n == mat_size);

    for (i = 0; i < mat_size; ++i) {
	double r;
	int j = 0;

#if IDX_BEGINS_WITH_ZERO
	genmat(i, ja, a, &b);
#else
	genmat(i + 1, ja, a, &b);
#endif
	r = b;

	for (j = 0; ja[j] != -1 && j < NZMAX; ++j) {
#if IDX_BEGINS_WITH_ZERO
	    r -= a[j] * x[ja[j]];
#else
	    r -= a[j] * x[ja[j]-1];
#endif
	}
	r = fabs(r);

	if (!(max_r >= r)) {
	    if (max_r >= 0) {
		max_r = r;
	    }
	}
    }

    fprintf(s, "Problem NO : %d\n", key);
    fprintf(s, "|b - Ax|_inf = %g", max_r);

    if (max_r < resmax) {
	fprintf (s, "  (OK)\n");
    }
    else {
	fprintf (s, "  (NG)\n");
    }
}

/*
 *  Simple sort
 */

static void
sort_PSC98(n, idx, perm)
    int n, *idx, *perm;
{
    int i, j;

    for (i = 0; i < n; ++i)
	perm[i] = i;

    for (i = 0; i < n; ++i) {
	int min = idx[perm[i]];
	int min_idx = i;
	for (j = i + 1; j < n; ++j) {
	    if (min > idx[perm[j]]) {
		min = idx[perm[j]];
		min_idx = j;
	    }
	}

	{
	    int t = perm[i];
	    perm[i] = perm[min_idx];
	    perm[min_idx] = t;
	}
    }
}

/*
 *  Right-hand term
 */

static double
trand1(i)
    int i;
{
    unsigned long t = 2 * i + 1;
    t *= 48828125;
    t *= 48828125;
    t *= 48828125;
    return .5 * (double)t / (double)(unsigned long)0x80000000;
}

static double
trand2(i)
    int i;
{
    unsigned long t = 2 * i + 1;
    t *= 1812433253;
    t *= 1812433253;
    t *= 1812433253;
    return .5 * (double)t / (double)(unsigned long)0x80000000;
}

static double
trand3(i)
    int i;
{
    unsigned long t = 2 * i + 1;
    t *= 1566083941;
    t *= 1566083941;
    t *= 1566083941;
    return .5 * (double)t / (double)(unsigned long)0x80000000;
}

static double
trand4(i)
    int i;
{
    unsigned long t = 2 * i + 1;
    t *= 69069;
    t *= 69069;
    t *= 69069;
    t *= 69069;
    t *= 69069;
    return .5 * (double)t / (double)(unsigned long)0x80000000;
}

static double
trand5(i)
    int i;
{
    unsigned long t = 2 * i + 1;
    t *= 1664525;
    t *= 1664525;
    t *= 1664525;
    return .5 * (double)t / (double)(unsigned long)0x80000000;
}


/***************  sample  ***************/

#define SIZE_X		127
#define SIZE_Y		127
#define PROB_SIZE	(SIZE_X * SIZE_Y)
#define MAX_NZERO	5
#define RESMAX		1e-5

static int
sample(i, ja, a, b)
    int i, *ja;
    double *a, *b;
{
    double inv_h = (SIZE_X + 1) * (SIZE_X + 1);
    double d0 = 4.0 * inv_h;
    double d1 = -1.0 * inv_h;

    if (i == -1) {
	ja[0] = PROB_SIZE;
	ja[1] = PROB_SIZE
	    + 2 * (PROB_SIZE - SIZE_Y)
	    + 2 * (PROB_SIZE - SIZE_X);
	ja[2] = MAX_NZERO;
	*b = RESMAX;

	return 0;
    }
    else if (i >= 0 && i < PROB_SIZE) {
	int k = 0;

	if (i - SIZE_X >= 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - SIZE_X;
#else
	    ja[k] = i + 1 - SIZE_X;
#endif
	    a[k] = d1;
	    ++k;
	}
	if (i % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - 1;
#else
	    ja[k] = i;
#endif
	    a[k] = d1;
	    ++k;
	}
#if IDX_BEGINS_WITH_ZERO
	ja[k] = i;
#else
	ja[k] = i + 1;
#endif
	a[k] = d0;
	++k;
	if ((i + 1) % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + 1;
#else
	    ja[k] = i + 2;
#endif
	    a[k] = d1;
	    ++k;
	}
	if (i + SIZE_X < PROB_SIZE) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + SIZE_X;
#else
	    ja[k] = i + 1 + SIZE_X;
#endif
	    a[k] = d1;
	    ++k;
	}
	for (; k < NZMAX; ++k) {
	    ja[k] = -1;
	    a[k] = 0.0;
	}

	/**  right-hand side  **/

	*b = trand3(i);

	return 0;
    }
    else {
	perror("sample");

	return 1;
    }
}

#undef SIZE_X
#undef SIZE_Y
#undef PROB_SIZE
#undef MAX_NZERO
#undef RESMAX


/***************  honsen1  ***************/

/*
 *  2-dimensional Poisson Equation with Dirichlet Condition
 *
 *  Spiral Ordering
 *
 *  Assume SIZE_X == SIZE_Y, and SIZE_X is odd.
 */

#define SIZE_X		1413
#define SIZE_Y		(SIZE_X)
#define MAX_NZERO	5
#define RESMAX		1e-5

static void
r_to_ij(n, r, i, j)
    int n, r;
    int *i, *j;
{
    int c, d, e, f;
    double d_n = (double)n;

    c = (int)floor(.5 * (d_n - sqrt(d_n * d_n - r)));
    d = r - 4 * c * (n - c);
    e = n - 2 * c - 1;
    if (e > 0)
	f = d / e;
    else
	f = 0;

    switch (f) {
	char errmsg[255];
    case 0:  /**  upper edge  **/
	*i = c;
	*j = c + d;
	break;
    case 1:  /**  right edge  **/
	*i = c + d - e;
	*j = n - 1 - c;
	break;
    case 2:  /**  lower edge  **/
	*i = n - 1 - c;
	*j = n - 1 - c - (d - 2 * e);
	break;
    case 3:  /**  left edge  **/
	*i = n - 1 - c - (d - 3 * e);
	*j = c;
	break;
    default:
	sprintf(errmsg, "r_to_ij: r = %d, f = %d", r, f);
	perror(errmsg);
	exit(1);
    }
}

#define MAX(i, j)  ((i) > (j) ? (i) : (j))

static int
ij_to_r(n, i, j)
    int n, i, j;
{
    int i1, j1, c;
    int n2 = n / 2;

    i1 = i - n2;
    j1 = j - n2;
    c = n2 - MAX(abs(i1), abs(j1));

    if (i - j <= 0) {
	return 4 * c * (n - c) + i + j - 2 * c;
    }
    else {
	return 4 * (c + 1)* (n - c - 1) - (i + j - 2 * c);
    }
}

#undef MAX

static int
honsen1(i, ja, a, rhs)
    int i, *ja;
    double *a, *rhs;
{
    double inv_h = (SIZE_X + 1) * (SIZE_X + 1);
    double d0 = 4.0 * inv_h;
    double d1 = -1.0 * inv_h;
    int N = SIZE_X * SIZE_Y;

    if (i == -1) {
	ja[0] = N;
	ja[1] = N + 2 * (N - SIZE_Y)
	    + 2 * (N - SIZE_X);
	ja[2] = MAX_NZERO;
	*rhs  = RESMAX;

	assert(N <= MATMAX);

	return 0;
    }
    else if (i >= 0 && i < N) {
	double b[MAX_NZERO];
	int jb[MAX_NZERO];
	int perm[MAX_NZERO];

	int k = 0;
	int s, t;
	int j;

	r_to_ij(SIZE_X, i, &s, &t);

	/**  Diagonal elements  **/

	b[k] = d0;
#if IDX_BEGINS_WITH_ZERO
	jb[k] = i;
#else
	jb[k] = i + 1;
#endif
	++k;

	/**  Subdiagonal elements  **/

	if (s - 1 >= 0) {
	    b[k] = d1;
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s - 1, t);
#else
	    jb[k] = ij_to_r(SIZE_X, s - 1, t) + 1;
#endif
	    ++k;
	}
	if (t - 1 >= 0) {
	    b[k] = d1;
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s, t - 1);
#else
	    jb[k] = ij_to_r(SIZE_X, s, t - 1) + 1;
#endif
	    ++k;
	}
	if (s + 1 < SIZE_X) {
	    b[k] = d1;
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s + 1, t);
#else
	    jb[k] = ij_to_r(SIZE_X, s + 1, t) + 1;
#endif
	    ++k;
	}
	if (t + 1 < SIZE_Y) {
	    b[k] = d1;
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s, t + 1);
#else
	    jb[k] = ij_to_r(SIZE_X, s, t + 1) + 1;
#endif
	    ++k;
	}

	sort_PSC98(k, jb, perm);

	for (j = 0; j < k; ++j) {
	    ja[j] = jb[perm[j]];
	    a[j] = b[perm[j]];
	}

	for (; k < NZMAX; ++k) {
	    ja[k] = -1;
	    a[k] = 0.0;
	}

	/**  right-hand side  **/

	*rhs = trand1(i);

	return 0;
    }
    else {
	perror("genmat");

	return 1;
    }
}

#undef SIZE_X
#undef SIZE_Y
#undef MAX_NZERO
#undef RESMAX


/***************  honsen2  ***************/

/*
 *  2-dimensional Poisson Equation with Dirichlet Condition
 *
 *  Spiral Ordering
 *
 *  Assume SIZE_X == SIZE_Y, and SIZE_X is odd.
 */

#define SIZE_X		511
#define SIZE_Y		(SIZE_X)
#define MAX_NZERO	5
#define RESMAX		1e-5

/*
 *  +----------+
 *  | k1       |
 *  |    +--+  |
 *  |    |k2|  |
 *  |    +--+  |
 *  |          |
 *  +----------+
 */

static void
get_coeff2(i, c)
    int i;
    double *c;
{
    double k1 = 1;
    double k2 = 100;
    int ix, iy;
    int cx1, cy1, cx2, cy2;

    iy = i / SIZE_X;
    ix = i % SIZE_X;

    cx1 = SIZE_X / 2;
    cy1 = SIZE_Y / 2;
    cx2 = 3 * SIZE_X / 4;
    cy2 = 3 * SIZE_Y / 4;

    if ((ix > cx1 && iy > cy1)
	&& (ix <= cx2 && iy <= cy2))
	c[0] = k2;
    else
	c[0] = k1;
	
    if ((ix >= cx1 && iy > cy1)
	&& (ix < cx2 && iy <= cy2))
	c[1] = k2;
    else
	c[1] = k1;

    if ((ix > cx1 && iy >= cy1)
	&& (ix <= cx2 && iy < cy2))
	c[2] = k2;
    else
	c[2] = k1;

    if ((ix >= cx1 && iy >= cy1)
	&& (ix < cx2 && iy < cy2))
	c[3] = k2;
    else
	c[3] = k1;
}

/*
 *  Spiral problem
 */

static int
honsen2(i, ja, a, rhs)
    int i, *ja;
    double *a, *rhs;
{
    double inv_h = (SIZE_X + 1) * (SIZE_X + 1);
    int N = SIZE_X * SIZE_Y;

    if (i == -1) {
	ja[0] = N;
	ja[1] = N + 2 * (N - SIZE_Y)
	    + 2 * (N - SIZE_X);
	ja[2] = MAX_NZERO;
	*rhs  = RESMAX;

	assert(N <= MATMAX);

	return 0;
    }
    else if (i >= 0 && i < N) {
	double b[MAX_NZERO];
	int jb[MAX_NZERO];
	int perm[MAX_NZERO];

	int k = 0;
	int s, t;
	int j;
	double c[4];
	double c0, c1, c2, c3;

	r_to_ij(SIZE_X, i, &s, &t);

	/*
	 *  c[2] | c[3]
	 *  -----+-----
	 *  c[0] | c[1]
	 */

	get_coeff2(t * SIZE_X + s, c);

	c0 = c[0] * inv_h;
	c1 = c[1] * inv_h;
	c2 = c[2] * inv_h;
	c3 = c[3] * inv_h;

	/**  Diagonal elements  **/

	b[k] = c0 + c1 + c2 + c3;
#if IDX_BEGINS_WITH_ZERO
	jb[k] = i;
#else
	jb[k] = i + 1;
#endif
	++k;

	/**  Subdiagonal elements  **/

	if (s - 1 >= 0) {
	    b[k] = -.5 * (c0 + c2);
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s - 1, t);
#else
	    jb[k] = ij_to_r(SIZE_X, s - 1, t) + 1;
#endif
	    ++k;
	}
	if (t - 1 >= 0) {
	    b[k] = -.5 * (c0 + c1);
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s, t - 1);
#else
	    jb[k] = ij_to_r(SIZE_X, s, t - 1) + 1;
#endif
	    ++k;
	}
	if (s + 1 < SIZE_X) {
	    b[k] = -.5 * (c1 + c3);
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s + 1, t);
#else
	    jb[k] = ij_to_r(SIZE_X, s + 1, t) + 1;
#endif
	    ++k;
	}
	if (t + 1 < SIZE_Y) {
	    b[k] = -.5 * (c2 + c3);
#if IDX_BEGINS_WITH_ZERO
	    jb[k] = ij_to_r(SIZE_X, s, t + 1);
#else
	    jb[k] = ij_to_r(SIZE_X, s, t + 1) + 1;
#endif
	    ++k;
	}

	sort_PSC98(k, jb, perm);

	for (j = 0; j < k; ++j) {
	    ja[j] = jb[perm[j]];
	    a[j] = b[perm[j]];
	}

	for (; k < NZMAX; ++k) {
	    ja[k] = -1;
	    a[k] = 0.0;
	}

	/**  right-hand side  **/

	*rhs = trand2(i);

	return 0;
    }
    else {
	perror("genmat");

	return 1;
    }
}

#undef SIZE_X
#undef SIZE_Y
#undef MAX_NZERO
#undef RESMAX


/***************  honsen3  ***************/

/*
 *  2-dimensional Poisson Equation with Dirichlet Condition
 *
 *  9-point stencil
 */

#define SIZE_X		1413
#define SIZE_Y		1413
#define PROB_SIZE	(SIZE_X * SIZE_Y)
#define MAX_NZERO	9
#define RESMAX		1e-5

static double
rhs3(i)
    int i;
{
    double q1, q2, q3, q4, q5;
    if (i - SIZE_X >= 0)
	q1 = trand3(i - SIZE_X);
    else
	q1 = 0.0;

    if (i % SIZE_X > 0)
	q2 = trand3(i - 1);
    else
	q2 = 0.0;

    q3 = trand3(i);

    if ((i + 1) % SIZE_X > 0)
	q4 = trand3(i + 1);
    else
	q4 = 0.0;

    if (i + SIZE_X < PROB_SIZE)
	q5 = trand3(i + SIZE_X);
    else
	q5 = 0.0;

    return q1 + q2 + 8.0 * q3 + q4 + q5;
}

static int
honsen3(i, ja, a, b)
    int i, *ja;
    double *a, *b;
{
    double inv_h = (SIZE_X + 1) * (SIZE_X + 1);
    double d0 = 40.0 * inv_h;
    double d1 = -8.0 * inv_h;
    double d2 = -2.0 * inv_h;

    if (i == -1) {
	ja[0] = PROB_SIZE;
	ja[1] = PROB_SIZE
	    + 2 * (PROB_SIZE - SIZE_Y)
	    + 2 * (PROB_SIZE - SIZE_X)
	    + 4 * (PROB_SIZE - SIZE_X - SIZE_Y + 1);
	ja[2] = MAX_NZERO;
	*b    = RESMAX;

	assert(PROB_SIZE <= MATMAX);

	return 0;
    }
    else if (i >= 0 && i < PROB_SIZE) {
	int k = 0;

	if (i - SIZE_X >= 0 && i % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - SIZE_X - 1;
#else
	    ja[k] = i - SIZE_X - 1 + 1;
#endif
	    a[k] = d2;
	    ++k;
	}
	if (i - SIZE_X >= 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - SIZE_X;
#else
	    ja[k] = i - SIZE_X + 1;
#endif
	    a[k] = d1;
	    ++k;
	}
	if (i - SIZE_X >= 0 && (i + 1) % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - SIZE_X + 1;
#else
	    ja[k] = i - SIZE_X + 1 + 1;
#endif
	    a[k] = d2;
	    ++k;
	}
	if (i % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - 1;
#else
	    ja[k] = i - 1 + 1;
#endif
	    a[k] = d1;
	    ++k;
	}
#if IDX_BEGINS_WITH_ZERO
	ja[k] = i;
#else
	ja[k] = i + 1;
#endif
	a[k] = d0;
	++k;
	if ((i + 1) % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + 1;
#else
	    ja[k] = i + 1 + 1;
#endif
	    a[k] = d1;
	    ++k;
	}
	if (i + SIZE_X < PROB_SIZE && i % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + SIZE_X - 1;
#else
	    ja[k] = i + SIZE_X - 1 + 1;
#endif
	    a[k] = d2;
	    ++k;
	}
	if (i + SIZE_X < PROB_SIZE) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + SIZE_X;
#else
	    ja[k] = i + SIZE_X + 1;
#endif
	    a[k] = d1;
	    ++k;
	}
	if (i + SIZE_X < PROB_SIZE && (i + 1) % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + SIZE_X + 1;
#else
	    ja[k] = i + SIZE_X + 1 + 1;
#endif
	    a[k] = d2;
	    ++k;
	}
	for (; k < NZMAX; ++k) {
	    ja[k] = -1;
	    a[k] = 0.0;
	}

	/**  right-hand side  **/

	*b = rhs3(i);

	return 0;
    }
    else {
	perror("genmat");

	return 1;
    }
}

#undef SIZE_X
#undef SIZE_Y
#undef PROB_SIZE
#undef MAX_NZERO
#undef RESMAX


/***************  honsen4  ***************/

/*
 *  3-dimensional Poisson Equation with Dirichlet Condition
 *
 *  Non-uniform diffusion constant
 */

#define SIZE_X		101
#define SIZE_Y		101
#define SIZE_Z		101
#define PROB_SIZE	(SIZE_X * SIZE_Y * SIZE_Z)
#define SIZE_XY		(SIZE_X * SIZE_Y)
#define MAX_NZERO	7
#define	RESMAX		1e-5

/*
 *  +----------+
 *  | k1       |
 *  |    +--+  |
 *  |    |k2|  |
 *  |    +--+  |
 *  |          |
 *  +----------+
 */

static void
get_coeff4(i, c)
    int i;
    double *c;
{
    double k1 = 1;
    double k2 = 100;
    int ix, iy, iz;
    int cx1, cy1, cz1, cx2, cy2, cz2;

    iz = i / SIZE_XY;
    iy = (i % SIZE_XY) / SIZE_X;
    ix = i % SIZE_XY % SIZE_X;

    cx1 = SIZE_X / 2;
    cy1 = SIZE_Y / 2;
    cz1 = SIZE_Z / 2;
    cx2 = 3 * SIZE_X / 4;
    cy2 = 3 * SIZE_Y / 4;
    cz2 = 3 * SIZE_Z / 4;

    /*
     *  c[2] | c[3]
     *  -----+-----
     *  c[0] | c[1]
     *
     *     c[6] | c[7]
     *     -----+-----
     *     c[4] | c[5]
     */

    if ((ix > cx1 && iy > cy1 && iz > cz1)
	&& (ix <= cx2 && iy <= cy2 && iz <= cz2))
	c[0] = k2;
    else
	c[0] = k1;
	
    if ((ix >= cx1 && iy > cy1 && iz > cz1)
	&& (ix < cx2 && iy <= cy2 && iz <= cz2))
	c[1] = k2;
    else
	c[1] = k1;

    if ((ix > cx1 && iy >= cy1 && iz > cz1)
	&& (ix <= cx2 && iy < cy2 && iz <= cz2))
	c[2] = k2;
    else
	c[2] = k1;

    if ((ix >= cx1 && iy >= cy1 && iz > cz1)
	&& (ix < cx2 && iy < cy2 && iz <= cz2))
	c[3] = k2;
    else
	c[3] = k1;

    if ((ix > cx1 && iy > cy1 && iz >= cz1)
	&& (ix <= cx2 && iy <= cy2 && iz < cz2))
	c[4] = k2;
    else
	c[4] = k1;
	
    if ((ix >= cx1 && iy > cy1 && iz >= cz1)
	&& (ix < cx2 && iy <= cy2 && iz < cz2))
	c[5] = k2;
    else
	c[5] = k1;

    if ((ix > cx1 && iy >= cy1 && iz >= cz1)
	&& (ix <= cx2 && iy < cy2 && iz < cz2))
	c[6] = k2;
    else
	c[6] = k1;

    if ((ix >= cx1 && iy >= cy1 && iz >= cz1)
	&& (ix < cx2 && iy < cy2 && iz < cz2))
	c[7] = k2;
    else
	c[7] = k1;
}

static int
honsen4(i, ja, a, b)
    int i, *ja;
    double *a, *b;
{
    double inv_h = (SIZE_X + 1) * (SIZE_X + 1);

    if (i == -1) {
	ja[0] = PROB_SIZE;
	ja[1] = PROB_SIZE
	    + 2 * (PROB_SIZE - SIZE_XY)
	    + 2 * (PROB_SIZE - SIZE_Y * SIZE_Z)
	    + 2 * (PROB_SIZE - SIZE_Z * SIZE_X);
	ja[2] = MAX_NZERO;
	*b    = RESMAX;

	assert(PROB_SIZE <= MATMAX);

	return 0;
    }
    else if (i >= 0 && i < PROB_SIZE) {
	int k = 0;
	double c[8];
	double c0, c1, c2, c3;
	double c4, c5, c6, c7;

	/*
	 *  c[2] | c[3]
	 *  -----+-----
	 *  c[0] | c[1]
	 *
	 *     c[6] | c[7]
	 *     -----+-----
	 *     c[4] | c[5]
	 */

	get_coeff4(i, c);

	c0 = c[0] * inv_h;
	c1 = c[1] * inv_h;
	c2 = c[2] * inv_h;
	c3 = c[3] * inv_h;
	c4 = c[4] * inv_h;
	c5 = c[5] * inv_h;
	c6 = c[6] * inv_h;
	c7 = c[7] * inv_h;

	if (i >= SIZE_XY) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - SIZE_XY;
#else
	    ja[k] = i - SIZE_XY + 1;
#endif
	    a[k] = -.25 * (c0 + c1 + c2 + c3);
	    ++k;
	}
	if (i % SIZE_XY >= SIZE_X) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - SIZE_X;
#else
	    ja[k] = i - SIZE_X + 1;
#endif
	    a[k] = -.25 * (c0 + c1 + c4 + c5);
	    ++k;
	}
	if (i % SIZE_XY % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - 1;
#else
	    ja[k] = i - 1 + 1;
#endif
	    a[k] = -.25 * (c0 + c2 + c4 + c6);
	    ++k;
	}
#if IDX_BEGINS_WITH_ZERO
	ja[k] = i;
#else
	ja[k] = i + 1;
#endif
	a[k] = c0 + c1 + c2 + c3 + c4 + c5 + c6 + c7;
	++k;
	if ((i + 1) % SIZE_XY % SIZE_X > 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + 1;
#else
	    ja[k] = i + 1 + 1;
#endif
	    a[k] = -.25 * (c1 + c3 + c5 + c7);
	    ++k;
	}
	if (i % SIZE_XY < SIZE_XY - SIZE_X) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + SIZE_X;
#else
	    ja[k] = i + SIZE_X + 1;
#endif
	    a[k] = -.25 * (c2 + c3 + c6 + c7);
	    ++k;
	}
	if (i < PROB_SIZE - SIZE_XY) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + SIZE_XY;
#else
	    ja[k] = i + SIZE_XY + 1;
#endif
	    a[k] = -.25 * (c4 + c5 + c6 + c7);
	    ++k;
	}
	for (; k < NZMAX; ++k) {
	    ja[k] = -1;
	    a[k] = 0.0;
	}

	/**  right-hand side  **/

	*b = trand4(i);

	return 0;
    }
    else {
	perror("genmat");

	return 1;
    }
}

#undef SIZE_X
#undef SIZE_Y
#undef SIZE_Z
#undef SIZE_XY
#undef PROB_SIZE
#undef MAX_NZERO
#undef RESMAX


/***************  honsen5  ***************/

/*
 *  1-dimensional Poisson Equation with Dirichlet Condition
 *
 *  Non-uniform diffusion constant
 */

#define SIZE_X	        24575
#define PROB_SIZE	(SIZE_X)
#define MAX_NZERO	3
#define	RESMAX		1e-3

/*
 *     k1  k2 k1
 *  +-----+--+--+
 */

static void
get_coeff5(i, c)
    int i;
    double *c;
{
    double k1 = 1;
    double k2 = 10;
    int cx1, cx2;

    cx1 = SIZE_X / 2;
    cx2 = 3 * SIZE_X / 4;

    if ((i > cx1) && (i <= cx2))
	c[0] = k2;
    else
	c[0] = k1;
	
    if ((i >= cx1) && (i < cx2))
	c[1] = k2;
    else
	c[1] = k1;
}

int
honsen5(i, ja, a, b)
    int i, *ja;
    double *a, *b;
{
    double inv_h = (SIZE_X + 1) * (SIZE_X + 1);

    if (i == -1) {
	ja[0] = PROB_SIZE;
	ja[1] = PROB_SIZE + 2 * (PROB_SIZE - 1);
	ja[2] = MAX_NZERO;
	*b    = RESMAX;

	assert(PROB_SIZE <= MATMAX);

	return 0;
    }
    else if (i >= 0 && i < PROB_SIZE) {
	int k = 0;
	double c[2];
	double c0, c1;

	/*
	 *  c[0] | c[1]
	 */

	get_coeff5(i, c);

	c0 = c[0] * inv_h;
	c1 = c[1] * inv_h;

	if (i - 1 >= 0) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i - 1;
#else
	    ja[k] = i - 1 + 1;
#endif
	    a[k] = -c0;
	    ++k;
	}
#if IDX_BEGINS_WITH_ZERO
	ja[k] = i;
#else
	ja[k] = i + 1;
#endif
	a[k] = c0 + c1;
	++k;
	if (i + 1 < PROB_SIZE) {
#if IDX_BEGINS_WITH_ZERO
	    ja[k] = i + 1;
#else
	    ja[k] = i + 1 + 1;
#endif
	    a[k] = -c1;
	    ++k;
	}
	for (; k < NZMAX; ++k) {
	    ja[k] = -1;
	    a[k] = 0.0;
	}

	/**  right-hand side  **/

	*b = trand5(i);

	return 0;
    }
    else {
	perror("genmat");

	return 1;
    }
}

#undef SIZE_X
#undef PROB_SIZE
#undef MAX_NZERO
#undef RESMAX