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"
47 #include "sphinxbase/ckd_alloc.h"
48 
49 void
50 norm_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 
77 void
78 accum_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 
95 void
96 floor_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 }
113 void
114 floor_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 
126 void
127 band_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
146 float64
147 determinant(float32 **a, int32 n)
148 {
149  E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n");
150  return 0.0;
151 }
152 int32
153 invert(float32 **ainv, float32 **a, int32 n)
154 {
155  E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n");
156  return 0;
157 }
158 int32
159 solve(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 
165 void
166 matrixmultiply(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. */
181 float64
182 determinant(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 
208 int32
209 solve(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. */
236 int32
237 invert(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 
266 void
267 matrixmultiply(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 
280 void
281 outerproduct(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 
294 void
295 scalarmultiply(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 
308 void
309 matrixadd(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 }
#define ckd_calloc_2d(d1, d2, sz)
Macro for ckd_calloc_2d
Definition: ckd_alloc.h:270
Sphinx&#39;s memory allocation/deallocation routines.
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 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 int32 invert(float32 **out_ainv, float32 **a, int32 len)
Invert (if possible) a positive definite matrix with QR algorithm.
Definition: matrix.c:153
Implementation of logging routines.
Matrix and linear algebra functions.
SPHINXBASE_EXPORT void matrixadd(float32 **inout_a, float32 **b, int32 n)
Add A += B.
Definition: matrix.c:309
SPHINXBASE_EXPORT void band_nz_1d(float32 *v, uint32 d1, float32 band)
Ensures that non-zero values x such that -band &lt; x &lt; band, band &gt; 0 are set to -band if x &lt; 0 and ban...
Definition: matrix.c:127
SPHINXBASE_EXPORT float64 determinant(float32 **a, int32 len)
Calculate the determinant of a positive definite matrix.
Definition: matrix.c:147
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
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
#define E_FATAL(...)
Exit with non-zero status after error message.
Definition: err.h:81
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 floor_nz_1d(float32 *v, uint32 d1, float32 floor)
Floor 1-d array.
Definition: matrix.c:114
SPHINXBASE_EXPORT void norm_3d(float32 ***arr, uint32 d1, uint32 d2, uint32 d3)
Norm an array.
Definition: matrix.c:50
SPHINXBASE_EXPORT void floor_nz_3d(float32 ***m, uint32 d1, uint32 d2, uint32 d3, float32 floor)
Floor 3-d array.
Definition: matrix.c:96