/*  atou1.c    CCMATH mathematics library source code.
 *
 *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
 *  This code may be redistributed under the terms of the GNU library
 *  public license (LGPL). ( See the lgpl.license file for details.)
 * ------------------------------------------------------------------------
 */
#include <stdlib.h>
void atou1(double *a, int m, int n)
{
    double *p0, *p, *q, *w;

    int i, j, k, mm;

    double s, h;

    w = (double *)calloc(m, sizeof(double));
    p0 = a + n * n - 1;
    i = n - 1;
    mm = m - n;
    if (mm == 0) {
	*p0 = 1.;
	p0 -= n + 1;
	--i;
	++mm;
    }
    for (; i >= 0; --i, ++mm, p0 -= n + 1) {
	if (*p0 != 0.) {
	    for (j = 0, p = p0 + n; j < mm; p += n)
		w[j++] = *p;
	    h = *p0;
	    *p0 = 1. - h;
	    for (j = 0, p = p0 + n; j < mm; p += n)
		*p = -h * w[j++];
	    for (k = i + 1, q = p0 + 1; k < n; ++k) {
		for (j = 0, p = q + n, s = 0.; j < mm; p += n)
		    s += w[j++] * *p;
		s *= h;
		for (j = 0, p = q + n; j < mm; p += n)
		    *p -= s * w[j++];
		*q++ = -s;
	    }
	}
	else {
	    *p0 = 1.;
	    for (j = 0, p = p0 + n, q = p0 + 1; j < mm; ++j, p += n)
		*p = *q++ = 0.;
	}
    }
    free(w);
}