ldumat.c

Go to the documentation of this file.

/*  ldumat.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 ldumat(double *a, double *u, int m, int n)
{
    double *p0, *q0, *p, *q, *w;
 
    int i, j, k, mm;
 
    double s, h;
 
    w = (double *)calloc(m, sizeof(double));
    for (i = 0, mm = m * m, q = u; i < mm; ++i)
        *q++ = 0.;
    p0 = a + n * n - 1;
    q0 = u + m * m - 1;
    mm = m - n;
    i = n - 1;
    for (j = 0; j < mm; ++j, q0 -= m + 1)
        *q0 = 1.;
    if (mm == 0) {
        p0 -= n + 1;
        *q0 = 1.;
        q0 -= m + 1;
        --i;
        ++mm;
    }
    for (; i >= 0; --i, ++mm, p0 -= n + 1, q0 -= m + 1) {
        if (*p0 != 0.) {
            for (j = 0, p = p0 + n, h = 1.; j < mm; p += n)
                w[j++] = *p;
            h = *p0;
            *q0 = 1. - h;
            for (j = 0, q = q0 + m; j < mm; q += m)
                *q = -h * w[j++];
            for (k = i + 1, q = q0 + 1; k < m; ++k) {
                for (j = 0, p = q + m, s = 0.; j < mm; p += m)
                    s += w[j++] * *p;
                s *= h;
                for (j = 0, p = q + m; j < mm; p += m)
                    *p -= s * w[j++];
                *q++ = -s;
            }
        }
        else {
            *q0 = 1.;
            for (j = 0, p = q0 + 1, q = q0 + m; j < mm; ++j, q += m)
                *q = *p++ = 0.;
        }
    }
    free(w);
}