path: root/src/libs/lprof/cmsmatn.cpp

/* */
/*  Little cms - profiler construction set */
/*  Copyright (C) 1998-2001 Marti Maria <marti@littlecms.com> */
/* */
/* THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, */
/* EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY */
/* WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. */
/* */
/* IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL, */
/* INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
/* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, */
/* WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF */
/* LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE */
/* OF THIS SOFTWARE. */
/* */
/* This file is free software; you can redistribute it and/or modify it */
/* under the terms of the GNU General Public License as published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, but */
/* WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
/* */
/* As a special exception to the GNU General Public License, if you */
/* distribute this file as part of a program that contains a */
/* configuration script generated by Autoconf, you may include it under */
/* the same distribution terms that you use for the rest of that program. */
/* */
/* Version 1.09a */


#include "lcmsprf.h"


LPMATN      cdecl MATNalloc(int Rows, int Cols);
void        cdecl MATNfree (LPMATN mat);
LPMATN      cdecl MATNmult(LPMATN a1, LPMATN a2);
double      cdecl MATNcross(LPMATN a);
void        cdecl MATNscalar (LPMATN a, double scl, LPMATN b);
LPMATN      cdecl MATNtranspose (LPMATN a);
BOOL        cdecl MATNsolve(LPMATN a, LPMATN b);


/* ------------------------------------------------------------ Implementation */

/* Free matrix */

void MATNfree(LPMATN mat)
{
    int i;

    if (mat == NULL) return;

    for (i = 0; i < mat->Rows; i++)
    {
        if (mat -> Values[i] != NULL)
            free (mat->Values[i]);
    }

    free(mat->Values);
    free(mat);
}


/* Allocate (and Zero) a new matrix */

LPMATN MATNalloc(int Rows, int Cols)
{
    int i;

    LPMATN mat = (LPMATN) malloc (sizeof (MATN));
    if (mat == NULL) return mat;

    ZeroMemory(mat, sizeof(MATN));

    mat->Rows = Rows;
    mat->Cols = Cols;
    mat->Values = (double**) malloc(Rows * sizeof (double*));

    if (mat->Values == NULL) {
        free(mat);
        return NULL;
    }

    ZeroMemory(mat -> Values, Rows * sizeof (double*));

    for (i = 0; i < Rows; i++)
    {
        mat-> Values [i] = (double*) malloc(Cols * sizeof (double));
        if (mat -> Values[i] == NULL) {
            MATNfree(mat);
            return NULL;
        }

    }

    return mat;
}

#define DO_SWAP(a, b, tmp)      { tmp = (a); (a) = (b); (b) = tmp; }

/* Gauss-Jordan elimination. There is also a more */
/* exahustive non-singular matrix checking part. */

BOOL MATNsolve(LPMATN a, LPMATN b)
{
    BOOL      status;
    int       n = a->Rows;
    int       i, iCol=0, iRow=0, j, k;
    double    fMax, fAbs, fSave, fInf, temp;
    int*      aiColIndex;
    int*      aiRowIndex=0;
    int*      aiPivoted=0;


    if (a->Rows != a->Cols) return false;

    status = false;
    if((aiColIndex = (int*) malloc(n * sizeof(int))) == NULL)
        goto GotError;

    if((aiRowIndex = (int*) malloc(n * sizeof(int))) == NULL)
        goto GotError;

    if((aiPivoted = (int*) malloc(n * sizeof(int))) == NULL)
        goto GotError;

    ZeroMemory(aiPivoted, n * sizeof(int));


    for(i = 0; i < n; i++) {

        /* search matrix (excluding pivoted rows) for maximum absolute entry */

        fMax = 0.0;
        for (j = 0; j < n; j++)
            if (aiPivoted[j] != 1)
                for (k = 0; k < n; k++)
                {
                   fAbs = fabs(a->Values[j][k]);
                   if (fAbs >= fMax) {

                        fMax = fAbs;
                        iRow = j;
                        iCol = k;
                    }
                    else
                    if (aiPivoted[k] > 1) {

                        status = false;
                        goto GotError;
                    }
                }

        aiPivoted[iCol]++;

        /* swap rows so that A[iCol][iCol] contains the pivot entry */

        if (iRow != iCol) {

            for(j = 0; j < n; j++)
                DO_SWAP(a->Values[iRow][j], a->Values[iCol][j], temp)

            DO_SWAP(b->Values[iRow][0], b->Values[iCol][0], temp)
        }

        /* keep track of the permutations of the rows */

        aiRowIndex[i] = iRow;
        aiColIndex[i] = iCol;

        if (a->Values[iCol][iCol] == 0.0)
        {
            status = false;
            goto GotError;
        }

       /* scale the row so that the pivot entry is 1 */

        fInf = 1.0 / a->Values[iCol][iCol];
        a->Values[iCol][iCol] = 1.0;

                for(j = 0; j < n; j++)
                        a->Values[iCol][j] *= fInf;

                b->Values[iCol][0] *= fInf;

                /* zero out the pivot column locations in the other rows */

                for(j = 0; j < n; j++)
                        if (j != iCol) {

                                fSave = a->Values[j][iCol];
                                a->Values[j][iCol] = 0.0;

                                for(k = 0; k < n; k++)
                                        a->Values[j][k] -= a->Values[iCol][k] * fSave;

                                b->Values[j][0] -= b->Values[iCol][0] * fSave;
                        }
        }

        /* reorder rows so that A[][] stores the inverse of the original matrix */

        for(i = n - 1; i >= 0; i--) {

                if(aiRowIndex[i] != aiColIndex[i])
                        for(j = 0; j < n; j++)
                                DO_SWAP(a->Values[j][aiRowIndex[i]], a->Values[j][aiColIndex[i]], temp)
        }

        status = true;

GotError:
        if(aiColIndex) free(aiColIndex);
        if(aiRowIndex) free(aiRowIndex);
        if(aiPivoted)  free(aiPivoted);
        return status;

}

#undef DO_SWAP


LPMATN MATNmult(LPMATN a1, LPMATN a2)
{
    int     i, j, k;
    LPMATN b;

    if (a1->Cols != a2->Rows)
        return NULL;

        b = MATNalloc (a1->Rows, a2->Cols);
        if (b == NULL)
                return NULL;

        for (i = 0; i < b->Rows; i++) {

                for (j = 0; j < b->Cols; j++) {

                        b->Values[i][j] = 0.0;

                        for (k = 0; k < a1->Cols; k++) {

                                b->Values[i][j] += a1->Values[i][k] * a2->Values[k][j];
                        }
                }
        }

        return b;
}


double MATNcross(LPMATN a)
{
        int i;
        double prod = 0.0;

        for (i = 0; i < a->Rows; i++) {

        prod += a->Values[i][0]*a->Values[i][0];
        }
        return prod;
}


void MATNscalar(LPMATN a, double scl, LPMATN b)
{
    int i, j;

    if (a->Rows != b->Rows || a->Cols != b->Cols)
        return;

    for (i = 0; i < a->Rows; i++) {

        for (j = 0; j < a->Cols; j++)
            b->Values[i][j] = a->Values[i][j] * scl;
    }
}


LPMATN MATNtranspose(LPMATN a)
{
    LPMATN b = MATNalloc(a->Cols, a->Rows);
    if (b != NULL) {

        int i, j;

        for (i = 0; i < a->Rows; i++)
        {
            for (j = 0; j < a->Cols; j++)
                b->Values[j][i] = a->Values [i][j];
        }
    }
    return b;
}



/* Used for debug purposes */
#ifdef DEBUG
void MATNprintf(char* name, LPMATN mat)
{
    int i, j;

    printf ("%s:\n", name);
    for (i= 0; i < mat->Rows; i++) {

        printf ("%3d", i);
        for (j = 0; j < mat->Cols; j++)
            printf ("  %.5f", mat->Values[i][j]);
        printf ("\n");
    }
}
#endif