test_tools_les.c

/*****************************************************************************
 *
 * MODULE:       Grass Gmath Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 *                 soerengebbert <at> gmx <dot> de
 *
 * PURPOSE:      functions to manage linear equation systems
 *                 part of the gmath library
 *
 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 *****************************************************************************/

#include "test_gmath_lib.h"
#include <stdlib.h>
#include <math.h>

/*!
 * \brief Allocate memory for a (not) quadratic linear equation system which
 * includes the Matrix A, vector x and vector b
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param cols int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_nquad_les(int rows, int cols, int type)
{
    return G_math_alloc_les_param(rows, cols, type, 2);
}

/*!
 * \brief Allocate memory for a (not) quadratic linear equation system which
 * includes the Matrix A and vector x
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param cols int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_nquad_les_Ax(int rows, int cols, int type)
{
    return G_math_alloc_les_param(rows, cols, type, 1);
}

/*!
 * \brief Allocate memory for a (not) quadratic linear equation system which
 * includes the Matrix A
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param cols int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_nquad_les_A(int rows, int cols, int type)
{
    return G_math_alloc_les_param(rows, cols, type, 0);
}

/*!
 * \brief Allocate memory for a (not) quadratic linear equation system which
 * includes the Matrix A, vector x and vector b
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param cols int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_nquad_les_Ax_b(int rows, int cols, int type)
{
    return G_math_alloc_les_param(rows, cols, type, 2);
}

/*!
 * \brief Allocate memory for a quadratic linear equation system which includes
 * the Matrix A, vector x and vector b
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_les(int rows, int type)
{
    return G_math_alloc_les_param(rows, rows, type, 2);
}

/*!
 * \brief Allocate memory for a quadratic linear equation system which includes
 * the Matrix A and vector x
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_les_Ax(int rows, int type)
{
    return G_math_alloc_les_param(rows, rows, type, 1);
}

/*!
 * \brief Allocate memory for a quadratic linear equation system which includes
 * the Matrix A
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_les_A(int rows, int type)
{
    return G_math_alloc_les_param(rows, rows, type, 0);
}

/*!
 * \brief Allocate memory for a quadratic linear equation system which includes
 * the Matrix A, vector x and vector b
 *
 * This function calls #G_math_alloc_les_param
 *
 * \param rows int
 * \param type int
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_les_Ax_b(int rows, int type)
{
    return G_math_alloc_les_param(rows, rows, type, 2);
}

/*!
 * \brief Allocate memory for a quadratic or not quadratic linear equation
 * system
 *
 * The type of the linear equation system must be G_MATH_NORMAL_LES for
 * a regular quadratic matrix or G_MATH_SPARSE_LES for a sparse matrix
 *
 * <p>
 * In case of G_MATH_NORMAL_LES
 *
 * A quadratic matrix of size rows*rows*sizeof(double) will allocated
 *
 * <p>
 * In case of G_MATH_SPARSE_LES
 *
 * a vector of size row will be allocated, ready to hold additional allocated
 * sparse vectors. each sparse vector may have a different size.
 *
 * Parameter parts defines which parts of the les should be allocated.
 * The number of columns and rows defines if the matrix is quadratic.
 *
 * \param rows int
 * \param cols int
 * \param type int
 * \param parts int -- 2 = A, x and b; 1 = A and x; 0 = A allocated
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_les_param(int rows, int cols, int type, int parts)
{
    G_math_les *les;

    if (type == G_MATH_SPARSE_LES)
        G_debug(2,
                "Allocate memory for a sparse linear equation system with %i "
                "rows\n",
                rows);
    else
        G_debug(2,
                "Allocate memory for a regular linear equation system with %i "
                "rows and %i cols\n",
                rows, cols);

    les = (G_math_les *)G_calloc(1, sizeof(G_math_les));
    les->x = NULL;
    les->b = NULL;

    if (parts > 0) {
        les->x = (double *)G_calloc(cols, sizeof(double));
    }

    if (parts > 1) {
        les->b = (double *)G_calloc(cols, sizeof(double));
    }

    les->A = NULL;
    les->data = NULL;
    les->Asp = NULL;
    les->rows = rows;
    les->cols = cols;
    les->symm = 0;
    les->bandwidth = cols;
    if (rows == cols)
        les->quad = 1;
    else
        les->quad = 0;

    if (type == G_MATH_SPARSE_LES) {
        les->Asp =
            (G_math_spvector **)G_calloc(rows, sizeof(G_math_spvector *));
        les->type = G_MATH_SPARSE_LES;
    }
    else {
        les->A = G_alloc_matrix(rows, cols);
        /*save the start pointer of the matrix */
        les->data = les->A[0];
        les->type = G_MATH_NORMAL_LES;
    }

    return les;
}

/***************** Floating point version ************************/

G_math_f_les *G_math_alloc_f_les(int rows, int type)
{
    return G_math_alloc_f_les_param(rows, rows, type, 2);
}

G_math_f_les *G_math_alloc_f_nquad_les_A(int rows, int cols, int type)
{
    return G_math_alloc_f_les_param(rows, cols, type, 0);
}

G_math_f_les *G_math_alloc_f_les_param(int rows, int cols, int type, int parts)
{
    G_math_f_les *les;

    G_debug(2,
            "Allocate memory for a regular float linear equation system with "
            "%i rows\n",
            rows);

    les = (G_math_f_les *)G_calloc(1, sizeof(G_math_f_les));
    les->x = NULL;
    les->b = NULL;

    if (parts > 0) {
        les->x = (float *)G_calloc(cols, sizeof(float));
    }

    if (parts > 1) {
        les->b = (float *)G_calloc(cols, sizeof(float));
    }

    les->A = NULL;
    les->data = NULL;
    les->rows = rows;
    les->cols = cols;
    les->symm = 0;
    les->bandwidth = cols;
    if (rows == cols)
        les->quad = 1;
    else
        les->quad = 0;

    les->A = G_alloc_fmatrix(rows, cols);
    /*save the start pointer of the matrix */
    les->data = les->A[0];
    les->type = G_MATH_NORMAL_LES;

    return les;
}

/*!
 * \brief Adds a sparse vector to a sparse linear equation system at position
 * row
 *
 * Return 1 for success and -1 for failure
 *
 * \param les G_math_les *
 * \param spvector G_math_spvector *
 * \param row int
 * \return int 0 success, -1 failure
 *
 * */
int G_math_add_spvector_to_les(G_math_les *les, G_math_spvector *spvector,
                               int row)
{

    if (les != NULL) {
        if (les->type != G_MATH_SPARSE_LES)
            return -1;

        if (les->rows > row) {
            G_debug(5,
                    "Add sparse vector %p to the sparse linear equation system "
                    "at row %i\n",
                    spvector, row);
            les->Asp[row] = spvector;
        }
        else
            return -1;
    }
    else {
        return -1;
    }

    return 1;
}

/*!
 *
 * \brief prints the linear equation system to stdout
 *
 * <p>
 * Format:
 * A*x = b
 *
 * <p>
 * Example
 \verbatim

 2 1 1 1 * 2 = 0.1
 1 2 0 0 * 3 = 0.2
 1 0 2 0 * 3 = 0.2
 1 0 0 2 * 2 = 0.1

 \endverbatim
 *
 * \param les G_math_les *
 * \return void
 *
 * */
void G_math_print_les(G_math_les *les)
{
    int i, j, k, out;

    if (les->type == G_MATH_SPARSE_LES) {
        for (i = 0; i < les->rows; i++) {
            for (j = 0; j < les->cols; j++) {
                out = 0;
                for (k = 0; k < les->Asp[i]->cols; k++) {
                    if (les->Asp[i]->index[k] == j) {
                        fprintf(stdout, "%4.5f ", les->Asp[i]->values[k]);
                        out = 1;
                    }
                }
                if (!out)
                    fprintf(stdout, "%4.5f ", 0.0);
            }
            if (les->x)
                fprintf(stdout, "  *  %4.5f", les->x[i]);
            if (les->b)
                fprintf(stdout, " =  %4.5f ", les->b[i]);

            fprintf(stdout, "\n");
        }
    }
    else {

        for (i = 0; i < les->rows; i++) {
            for (j = 0; j < les->cols; j++) {
                fprintf(stdout, "%4.5f ", les->A[i][j]);
            }
            if (les->x)
                fprintf(stdout, "  *  %4.5f", les->x[i]);
            if (les->b)
                fprintf(stdout, " =  %4.5f ", les->b[i]);

            fprintf(stdout, "\n");
        }
    }
    return;
}

/*!
 * \brief Release the memory of the linear equation system
 *
 * \param les G_math_les *
 * \return void
 *
 * */

void G_math_free_les(G_math_les *les)
{
    int i;

    if (les->type == G_MATH_SPARSE_LES)
        G_debug(2, "Releasing memory of a sparse linear equation system\n");
    else
        G_debug(2, "Releasing memory of a regular linear equation system\n");

    if (les) {

        if (les->x)
            G_free(les->x);
        if (les->b)
            G_free(les->b);

        if (les->type == G_MATH_SPARSE_LES) {

            if (les->Asp) {
                for (i = 0; i < les->rows; i++)
                    if (les->Asp[i])
                        G_math_free_spvector(les->Asp[i]);

                G_free(les->Asp);
            }
        }
        else {
            /*We don't know if the rows have been changed by pivoting,
             * so we restore the data pointer*/
            les->A[0] = les->data;
            G_free_matrix(les->A);
        }

        free(les);
    }

    return;
}

/*!
 * \brief Release the memory of the float linear equation system
 *
 * \param les G_math_f_les *
 * \return void
 *
 * */

void G_math_free_f_les(G_math_f_les *les)
{
    G_debug(2, "Releasing memory of a regular float linear equation system\n");

    if (les) {

        if (les->x)
            G_free(les->x);
        if (les->b)
            G_free(les->b);

        /*We don't know if the rows have been changed by pivoting,
         * so we restore the data pointer*/
        les->A[0] = les->data;
        G_free_fmatrix(les->A);

        free(les);
    }

    return;
}