Subversion Repositories Scribus

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2454 craig 1
//
2
//  Little cms - profiler construction set
3
//  Copyright (C) 1998-2001 Marti Maria
4
//
5
// THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND,
6
// EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
7
// WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
8
//
9
// IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL,
10
// INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
11
// OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
12
// WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF
13
// LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
14
// OF THIS SOFTWARE.
15
//
16
// This file is free software; you can redistribute it and/or modify it
17
// under the terms of the GNU General Public License as published by
18
// the Free Software Foundation; either version 2 of the License, or
19
// (at your option) any later version.
20
//
21
// This program is distributed in the hope that it will be useful, but
22
// WITHOUT ANY WARRANTY; without even the implied warranty of
23
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
24
// General Public License for more details.
25
//
26
// You should have received a copy of the GNU General Public License
27
// along with this program; if not, write to the Free Software
28
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
29
//
30
// As a special exception to the GNU General Public License, if you
31
// distribute this file as part of a program that contains a
32
// configuration script generated by Autoconf, you may include it under
33
// the same distribution terms that you use for the rest of that program.
34
//
35
// Version 1.09a
36
 
37
 
38
#include "lcmsprf.h"
39
 
40
 
41
LPMATN      cdecl MATNalloc(int Rows, int Cols);
42
void        cdecl MATNfree (LPMATN mat);
43
LPMATN      cdecl MATNmult(LPMATN a1, LPMATN a2);
44
double      cdecl MATNcross(LPMATN a);
45
void        cdecl MATNscalar (LPMATN a, double scl, LPMATN b);
46
LPMATN      cdecl MATNtranspose (LPMATN a);
47
BOOL        cdecl MATNsolve(LPMATN a, LPMATN b);
48
 
49
 
50
// ------------------------------------------------------------ Implementation
51
 
52
// Free matrix
53
 
54
void MATNfree(LPMATN mat)
55
{
56
    int i;
57
 
58
    if (mat == NULL) return;
59
 
60
    for (i = 0; i < mat->Rows; i++)
61
    {
62
        if (mat -> Values[i] != NULL)
63
            free (mat->Values[i]);
64
    }
65
 
66
    free(mat->Values);
67
    free(mat);
68
}
69
 
70
 
71
// Allocate (and Zero) a new matrix
72
 
73
LPMATN MATNalloc(int Rows, int Cols)
74
{
75
    int i;
76
 
77
    LPMATN mat = (LPMATN) malloc (sizeof (MATN));
78
    if (mat == NULL) return mat;
79
 
80
    ZeroMemory(mat, sizeof(MATN));
81
 
82
    mat->Rows = Rows;
83
    mat->Cols = Cols;
84
    mat->Values = (double**) malloc(Rows * sizeof (double*));
85
 
86
    if (mat->Values == NULL) {
87
        free(mat);
88
        return NULL;
89
    }
90
 
91
    ZeroMemory(mat -> Values, Rows * sizeof (double*));
92
 
93
    for (i = 0; i < Rows; i++)
94
    {
95
        mat-> Values [i] = (double*) malloc(Cols * sizeof (double));
96
        if (mat -> Values[i] == NULL) {
97
            MATNfree(mat);
98
            return NULL;
99
        }
100
 
101
    }
102
 
103
    return mat;
104
}
105
 
106
#define DO_SWAP(a, b, tmp)      { tmp = (a); (a) = (b); (b) = tmp; }
107
 
108
// Gauss-Jordan elimination. There is also a more
109
// exahustive non-singular matrix checking part.
110
 
111
BOOL MATNsolve(LPMATN a, LPMATN b)
112
{
113
    BOOL      status;
114
    int       n = a->Rows;
115
    int       i, iCol=0, iRow=0, j, k;
116
    double    fMax, fAbs, fSave, fInf, temp;
117
    int*      aiColIndex;
118
    int*      aiRowIndex=0;
119
    int*      aiPivoted=0;
120
 
121
 
122
    if (a->Rows != a->Cols) return FALSE;
123
 
124
    status = FALSE;
125
    if((aiColIndex = (int*) malloc(n * sizeof(int))) == NULL)
126
        goto GotError;
127
 
128
    if((aiRowIndex = (int*) malloc(n * sizeof(int))) == NULL)
129
        goto GotError;
130
 
131
    if((aiPivoted = (int*) malloc(n * sizeof(int))) == NULL)
132
        goto GotError;
133
 
134
    ZeroMemory(aiPivoted, n * sizeof(int));
135
 
136
 
137
    for(i = 0; i < n; i++) {
138
 
139
        // search matrix (excluding pivoted rows) for maximum absolute entry
140
 
141
        fMax = 0.0;
142
        for (j = 0; j < n; j++)
143
            if (aiPivoted[j] != 1)
144
                for (k = 0; k < n; k++)
145
                {
146
                   fAbs = fabs(a->Values[j][k]);
147
                   if (fAbs >= fMax) {
148
 
149
                        fMax = fAbs;
150
                        iRow = j;
151
                        iCol = k;
152
                    }
153
                    else
154
                    if (aiPivoted[k] > 1) {
155
 
156
                        status = FALSE;
157
                        goto GotError;
158
                    }
159
                }
160
 
161
        aiPivoted[iCol]++;
162
 
163
        // swap rows so that A[iCol][iCol] contains the pivot entry
164
 
165
        if (iRow != iCol) {
166
 
167
            for(j = 0; j < n; j++)
168
                DO_SWAP(a->Values[iRow][j], a->Values[iCol][j], temp)
169
 
170
            DO_SWAP(b->Values[iRow][0], b->Values[iCol][0], temp)
171
        }
172
 
173
        // keep track of the permutations of the rows
174
 
175
        aiRowIndex[i] = iRow;
176
        aiColIndex[i] = iCol;
177
 
178
        if (a->Values[iCol][iCol] == 0.0)
179
        {
180
            status = FALSE;
181
            goto GotError;
182
        }
183
 
184
       // scale the row so that the pivot entry is 1
185
 
186
        fInf = 1.0 / a->Values[iCol][iCol];
187
        a->Values[iCol][iCol] = 1.0;
188
 
189
                for(j = 0; j < n; j++)
190
                        a->Values[iCol][j] *= fInf;
191
 
192
                b->Values[iCol][0] *= fInf;
193
 
194
                // zero out the pivot column locations in the other rows
195
 
196
                for(j = 0; j < n; j++)
197
                        if (j != iCol) {
198
 
199
                                fSave = a->Values[j][iCol];
200
                                a->Values[j][iCol] = 0.0;
201
 
202
                                for(k = 0; k < n; k++)
203
                                        a->Values[j][k] -= a->Values[iCol][k] * fSave;
204
 
205
                                b->Values[j][0] -= b->Values[iCol][0] * fSave;
206
                        }
207
        }
208
 
209
        // reorder rows so that A[][] stores the inverse of the original matrix
210
 
211
        for(i = n - 1; i >= 0; i--) {
212
 
213
                if(aiRowIndex[i] != aiColIndex[i])
214
                        for(j = 0; j < n; j++)
215
                                DO_SWAP(a->Values[j][aiRowIndex[i]], a->Values[j][aiColIndex[i]], temp)
216
        }
217
 
218
        status = TRUE;
219
 
220
GotError:
221
        if(aiColIndex) free(aiColIndex);
222
        if(aiRowIndex) free(aiRowIndex);
223
        if(aiPivoted)  free(aiPivoted);
224
        return status;
225
 
226
}
227
 
228
#undef DO_SWAP
229
 
230
 
231
LPMATN MATNmult(LPMATN a1, LPMATN a2)
232
{
233
    int     i, j, k;
234
    LPMATN b;
235
 
236
    if (a1->Cols != a2->Rows)
237
        return NULL;
238
 
239
        b = MATNalloc (a1->Rows, a2->Cols);
240
        if (b == NULL)
241
                return NULL;
242
 
243
        for (i = 0; i < b->Rows; i++) {
244
 
245
                for (j = 0; j < b->Cols; j++) {
246
 
247
                        b->Values[i][j] = 0.0;
248
 
249
                        for (k = 0; k < a1->Cols; k++) {
250
 
251
                                b->Values[i][j] += a1->Values[i][k] * a2->Values[k][j];
252
                        }
253
                }
254
        }
255
 
256
        return b;
257
}
258
 
259
 
260
double MATNcross(LPMATN a)
261
{
262
        int i;
263
        double prod = 0.0;
264
 
265
        for (i = 0; i < a->Rows; i++) {
266
 
267
        prod += a->Values[i][0]*a->Values[i][0];
268
        }
269
        return prod;
270
}
271
 
272
 
273
void MATNscalar(LPMATN a, double scl, LPMATN b)
274
{
275
    int i, j;
276
 
277
    if (a->Rows != b->Rows || a->Cols != b->Cols)
278
        return;
279
 
280
    for (i = 0; i < a->Rows; i++) {
281
 
282
        for (j = 0; j < a->Cols; j++)
283
            b->Values[i][j] = a->Values[i][j] * scl;
284
    }
285
}
286
 
287
 
288
LPMATN MATNtranspose(LPMATN a)
289
{
290
    LPMATN b = MATNalloc(a->Cols, a->Rows);
291
    if (b != NULL) {
292
 
293
        int i, j;
294
 
295
        for (i = 0; i < a->Rows; i++)
296
        {
297
            for (j = 0; j < a->Cols; j++)
298
                b->Values[j][i] = a->Values [i][j];
299
        }
300
    }
301
    return b;
302
}
303
 
304
 
305
 
306
// Used for debug purposes
307
#ifdef DEBUG
308
void MATNprintf(char* name, LPMATN mat)
309
{
310
    int i, j;
311
 
312
    printf ("%s:\n", name);
313
    for (i= 0; i < mat->Rows; i++) {
314
 
315
        printf ("%3d", i);
316
        for (j = 0; j < mat->Cols; j++)
317
            printf ("  %.5f", mat->Values[i][j]);
318
        printf ("\n");
319
    }
320
}
321
#endif
322
 
323