SphinxBase 5prealpha
matrix.c
1/* -*- c-basic-offset: 4 -*- */
2/* ====================================================================
3 * Copyright (c) 1997-2006 Carnegie Mellon University. All rights
4 * reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 *
13 * 2. Redistributions in binary form must reproduce the above copyright
14 * notice, this list of conditions and the following disclaimer in
15 * the documentation and/or other materials provided with the
16 * distribution.
17 *
18 * This work was supported in part by funding from the Defense Advanced
19 * Research Projects Agency and the National Science Foundation of the
20 * United States of America, and the CMU Sphinx Speech Consortium.
21 *
22 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
23 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
24 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
25 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
26 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
27 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
28 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
29 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
30 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
31 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
32 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 *
34 * ====================================================================
35 *
36 */
37#include <string.h>
38#include <stdlib.h>
39
40#ifdef HAVE_CONFIG_H
41#include "config.h"
42#endif
43
44#include "sphinxbase/clapack_lite.h"
45#include "sphinxbase/matrix.h"
46#include "sphinxbase/err.h"
48
49void
50norm_3d(float32 ***arr,
51 uint32 d1,
52 uint32 d2,
53 uint32 d3)
54{
55 uint32 i, j, k;
56 float64 s;
57
58 for (i = 0; i < d1; i++) {
59 for (j = 0; j < d2; j++) {
60
61 /* compute sum (i, j) as over all k */
62 for (k = 0, s = 0; k < d3; k++) {
63 s += arr[i][j][k];
64 }
65
66 /* do 1 floating point divide */
67 s = 1.0 / s;
68
69 /* divide all k by sum over k */
70 for (k = 0; k < d3; k++) {
71 arr[i][j][k] *= s;
72 }
73 }
74 }
75}
76
77void
78accum_3d(float32 ***out,
79 float32 ***in,
80 uint32 d1,
81 uint32 d2,
82 uint32 d3)
83{
84 uint32 i, j, k;
85
86 for (i = 0; i < d1; i++) {
87 for (j = 0; j < d2; j++) {
88 for (k = 0; k < d3; k++) {
89 out[i][j][k] += in[i][j][k];
90 }
91 }
92 }
93}
94
95void
96floor_nz_3d(float32 ***m,
97 uint32 d1,
98 uint32 d2,
99 uint32 d3,
100 float32 floor)
101{
102 uint32 i, j, k;
103
104 for (i = 0; i < d1; i++) {
105 for (j = 0; j < d2; j++) {
106 for (k = 0; k < d3; k++) {
107 if ((m[i][j][k] != 0) && (m[i][j][k] < floor))
108 m[i][j][k] = floor;
109 }
110 }
111 }
112}
113void
114floor_nz_1d(float32 *v,
115 uint32 d1,
116 float32 floor)
117{
118 uint32 i;
119
120 for (i = 0; i < d1; i++) {
121 if ((v[i] != 0) && (v[i] < floor))
122 v[i] = floor;
123 }
124}
125
126void
127band_nz_1d(float32 *v,
128 uint32 d1,
129 float32 band)
130{
131 uint32 i;
132
133 for (i = 0; i < d1; i++) {
134 if (v[i] != 0) {
135 if ((v[i] > 0) && (v[i] < band)) {
136 v[i] = band;
137 }
138 else if ((v[i] < 0) && (v[i] > -band)) {
139 v[i] = -band;
140 }
141 }
142 }
143}
144
145#ifndef WITH_LAPACK
146float64
147determinant(float32 **a, int32 n)
148{
149 E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n");
150 return 0.0;
151}
152int32
153invert(float32 **ainv, float32 **a, int32 n)
154{
155 E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n");
156 return 0;
157}
158int32
159solve(float32 **a, float32 *b, float32 *out_x, int32 n)
160{
161 E_FATAL("No LAPACK library available, cannot solve linear equations (FIXME)\n");
162 return 0;
163}
164
165void
166matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
167{
168 int32 i, j, k;
169
170 memset(c[0], 0, n*n*sizeof(float32));
171 for (i = 0; i < n; ++i) {
172 for (j = 0; j < n; ++j) {
173 for (k = 0; k < n; ++k) {
174 c[i][k] += a[i][j] * b[j][k];
175 }
176 }
177 }
178}
179#else /* WITH_LAPACK */
180/* Find determinant through LU decomposition. */
181float64
182determinant(float32 ** a, int32 n)
183{
184 float32 **tmp_a;
185 float64 det;
186 char uplo;
187 int32 info, i;
188
189 /* a is assumed to be symmetric, so we don't need to switch the
190 * ordering of the data. But we do need to copy it since it is
191 * overwritten by LAPACK. */
192 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
193 memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
194
195 uplo = 'L';
196 spotrf_(&uplo, &n, tmp_a[0], &n, &info);
197 det = tmp_a[0][0];
198 /* det = prod(diag(l))^2 */
199 for (i = 1; i < n; ++i)
200 det *= tmp_a[i][i];
201 ckd_free_2d((void **)tmp_a);
202 if (info > 0)
203 return -1.0; /* Generic "not positive-definite" answer */
204 else
205 return det * det;
206}
207
208int32
209solve(float32 **a, /*Input : an n*n matrix A */
210 float32 *b, /*Input : a n dimesion vector b */
211 float32 *out_x, /*Output : a n dimesion vector x */
212 int32 n)
213{
214 char uplo;
215 float32 **tmp_a;
216 int32 info, nrhs;
217
218 /* a is assumed to be symmetric, so we don't need to switch the
219 * ordering of the data. But we do need to copy it since it is
220 * overwritten by LAPACK. */
221 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
222 memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
223 memcpy(out_x, b, n*sizeof(float32));
224 uplo = 'L';
225 nrhs = 1;
226 sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, out_x, &n, &info);
227 ckd_free_2d((void **)tmp_a);
228
229 if (info != 0)
230 return -1;
231 else
232 return info;
233}
234
235/* Find inverse by solving AX=I. */
236int32
237invert(float32 ** ainv, float32 ** a, int32 n)
238{
239 char uplo;
240 float32 **tmp_a;
241 int32 info, nrhs, i;
242
243 /* a is assumed to be symmetric, so we don't need to switch the
244 * ordering of the data. But we do need to copy it since it is
245 * overwritten by LAPACK. */
246 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
247 memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
248
249 /* Construct an identity matrix. */
250 memset(ainv[0], 0, sizeof(float32) * n * n);
251 for (i = 0; i < n; i++)
252 ainv[i][i] = 1.0;
253
254 uplo = 'L';
255 nrhs = n;
256 sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, ainv[0], &n, &info);
257
258 ckd_free_2d((void **)tmp_a);
259
260 if (info != 0)
261 return -1;
262 else
263 return info;
264}
265
266void
267matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
268{
269 char side, uplo;
270 float32 alpha;
271
272 side = 'L';
273 uplo = 'L';
274 alpha = 1.0;
275 ssymm_(&side, &uplo, &n, &n, &alpha, a[0], &n, b[0], &n, &alpha, c[0], &n);
276}
277
278#endif /* WITH_LAPACK */
279
280void
281outerproduct(float32 ** a, float32 * x, float32 * y, int32 len)
282{
283 int32 i, j;
284
285 for (i = 0; i < len; ++i) {
286 a[i][i] = x[i] * y[i];
287 for (j = i + 1; j < len; ++j) {
288 a[i][j] = x[i] * y[j];
289 a[j][i] = x[j] * y[i];
290 }
291 }
292}
293
294void
295scalarmultiply(float32 ** a, float32 x, int32 n)
296{
297 int32 i, j;
298
299 for (i = 0; i < n; ++i) {
300 a[i][i] *= x;
301 for (j = i+1; j < n; ++j) {
302 a[i][j] *= x;
303 a[j][i] *= x;
304 }
305 }
306}
307
308void
309matrixadd(float32 ** a, float32 ** b, int32 n)
310{
311 int32 i, j;
312
313 for (i = 0; i < n; ++i)
314 for (j = 0; j < n; ++j)
315 a[i][j] += b[i][j];
316}
Sphinx's memory allocation/deallocation routines.
SPHINXBASE_EXPORT void ckd_free_2d(void *ptr)
Free a 2-D array (ptr) previously allocated by ckd_calloc_2d.
Definition: ckd_alloc.c:255
#define ckd_calloc_2d(d1, d2, sz)
Macro for ckd_calloc_2d
Definition: ckd_alloc.h:270
Implementation of logging routines.
#define E_FATAL(...)
Exit with non-zero status after error message.
Definition: err.h:81
Matrix and linear algebra functions.
SPHINXBASE_EXPORT float64 determinant(float32 **a, int32 len)
Calculate the determinant of a positive definite matrix.
Definition: matrix.c:147
SPHINXBASE_EXPORT void floor_nz_1d(float32 *v, uint32 d1, float32 floor)
Floor 1-d array.
Definition: matrix.c:114
SPHINXBASE_EXPORT int32 solve(float32 **a, float32 *b, float32 *out_x, int32 n)
Solve (if possible) a positive-definite system of linear equations AX=B for X.
Definition: matrix.c:159
SPHINXBASE_EXPORT void matrixadd(float32 **inout_a, float32 **b, int32 n)
Add A += B.
Definition: matrix.c:309
SPHINXBASE_EXPORT void floor_nz_3d(float32 ***m, uint32 d1, uint32 d2, uint32 d3, float32 floor)
Floor 3-d array.
Definition: matrix.c:96
SPHINXBASE_EXPORT void band_nz_1d(float32 *v, uint32 d1, float32 band)
Ensures that non-zero values x such that -band < x < band, band > 0 are set to -band if x < 0 and ban...
Definition: matrix.c:127
SPHINXBASE_EXPORT int32 invert(float32 **out_ainv, float32 **a, int32 len)
Invert (if possible) a positive definite matrix with QR algorithm.
Definition: matrix.c:153
SPHINXBASE_EXPORT void outerproduct(float32 **out_a, float32 *x, float32 *y, int32 len)
Calculate the outer product of two vectors.
Definition: matrix.c:281
SPHINXBASE_EXPORT void accum_3d(float32 ***out, float32 ***in, uint32 d1, uint32 d2, uint32 d3)
Floor 3-d array.
Definition: matrix.c:78
SPHINXBASE_EXPORT void matrixmultiply(float32 **out_c, float32 **a, float32 **b, int32 n)
Multiply C=AB where A and B are symmetric matrices.
Definition: matrix.c:166
SPHINXBASE_EXPORT void scalarmultiply(float32 **inout_a, float32 x, int32 n)
Multiply a symmetric matrix by a constant in-place.
Definition: matrix.c:295
SPHINXBASE_EXPORT void norm_3d(float32 ***arr, uint32 d1, uint32 d2, uint32 d3)
Norm an array.
Definition: matrix.c:50