SphinxBase
0.6
|
00001 /* -*- c-basic-offset: 4 -*- */ 00002 /* ==================================================================== 00003 * Copyright (c) 1997-2006 Carnegie Mellon University. All rights 00004 * reserved. 00005 * 00006 * Redistribution and use in source and binary forms, with or without 00007 * modification, are permitted provided that the following conditions 00008 * are met: 00009 * 00010 * 1. Redistributions of source code must retain the above copyright 00011 * notice, this list of conditions and the following disclaimer. 00012 * 00013 * 2. Redistributions in binary form must reproduce the above copyright 00014 * notice, this list of conditions and the following disclaimer in 00015 * the documentation and/or other materials provided with the 00016 * distribution. 00017 * 00018 * This work was supported in part by funding from the Defense Advanced 00019 * Research Projects Agency and the National Science Foundation of the 00020 * United States of America, and the CMU Sphinx Speech Consortium. 00021 * 00022 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 00023 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 00024 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00025 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY 00026 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00027 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 00028 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00029 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00030 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00031 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 00032 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00033 * 00034 * ==================================================================== 00035 * 00036 */ 00037 #include <string.h> 00038 #include <stdlib.h> 00039 00040 #ifdef HAVE_CONFIG_H 00041 #include "config.h" 00042 #endif 00043 00044 #include "sphinxbase/clapack_lite.h" 00045 #include "sphinxbase/matrix.h" 00046 #include "sphinxbase/err.h" 00047 #include "sphinxbase/ckd_alloc.h" 00048 00049 void 00050 norm_3d(float32 ***arr, 00051 uint32 d1, 00052 uint32 d2, 00053 uint32 d3) 00054 { 00055 uint32 i, j, k; 00056 float64 s; 00057 00058 for (i = 0; i < d1; i++) { 00059 for (j = 0; j < d2; j++) { 00060 00061 /* compute sum (i, j) as over all k */ 00062 for (k = 0, s = 0; k < d3; k++) { 00063 s += arr[i][j][k]; 00064 } 00065 00066 /* do 1 floating point divide */ 00067 s = 1.0 / s; 00068 00069 /* divide all k by sum over k */ 00070 for (k = 0; k < d3; k++) { 00071 arr[i][j][k] *= s; 00072 } 00073 } 00074 } 00075 } 00076 00077 void 00078 accum_3d(float32 ***out, 00079 float32 ***in, 00080 uint32 d1, 00081 uint32 d2, 00082 uint32 d3) 00083 { 00084 uint32 i, j, k; 00085 00086 for (i = 0; i < d1; i++) { 00087 for (j = 0; j < d2; j++) { 00088 for (k = 0; k < d3; k++) { 00089 out[i][j][k] += in[i][j][k]; 00090 } 00091 } 00092 } 00093 } 00094 00095 void 00096 floor_nz_3d(float32 ***m, 00097 uint32 d1, 00098 uint32 d2, 00099 uint32 d3, 00100 float32 floor) 00101 { 00102 uint32 i, j, k; 00103 00104 for (i = 0; i < d1; i++) { 00105 for (j = 0; j < d2; j++) { 00106 for (k = 0; k < d3; k++) { 00107 if ((m[i][j][k] != 0) && (m[i][j][k] < floor)) 00108 m[i][j][k] = floor; 00109 } 00110 } 00111 } 00112 } 00113 void 00114 floor_nz_1d(float32 *v, 00115 uint32 d1, 00116 float32 floor) 00117 { 00118 uint32 i; 00119 00120 for (i = 0; i < d1; i++) { 00121 if ((v[i] != 0) && (v[i] < floor)) 00122 v[i] = floor; 00123 } 00124 } 00125 00126 void 00127 band_nz_1d(float32 *v, 00128 uint32 d1, 00129 float32 band) 00130 { 00131 uint32 i; 00132 00133 for (i = 0; i < d1; i++) { 00134 if (v[i] != 0) { 00135 if ((v[i] > 0) && (v[i] < band)) { 00136 v[i] = band; 00137 } 00138 else if ((v[i] < 0) && (v[i] > -band)) { 00139 v[i] = -band; 00140 } 00141 } 00142 } 00143 } 00144 00145 #ifndef WITH_LAPACK 00146 float64 00147 determinant(float32 **a, int32 n) 00148 { 00149 E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n"); 00150 return 0.0; 00151 } 00152 int32 00153 invert(float32 **ainv, float32 **a, int32 n) 00154 { 00155 E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n"); 00156 return 0; 00157 } 00158 int32 00159 solve(float32 **a, float32 *b, float32 *out_x, int32 n) 00160 { 00161 E_FATAL("No LAPACK library available, cannot solve linear equations (FIXME)\n"); 00162 return 0; 00163 } 00164 00165 void 00166 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n) 00167 { 00168 int32 i, j, k; 00169 00170 memset(c[0], 0, n*n*sizeof(float32)); 00171 for (i = 0; i < n; ++i) { 00172 for (j = 0; j < n; ++j) { 00173 for (k = 0; k < n; ++k) { 00174 c[i][k] += a[i][j] * b[j][k]; 00175 } 00176 } 00177 } 00178 } 00179 #else /* WITH_LAPACK */ 00180 /* Find determinant through LU decomposition. */ 00181 float64 00182 determinant(float32 ** a, int32 n) 00183 { 00184 float32 **tmp_a; 00185 float64 det; 00186 char uplo; 00187 int32 info, i; 00188 00189 /* a is assumed to be symmetric, so we don't need to switch the 00190 * ordering of the data. But we do need to copy it since it is 00191 * overwritten by LAPACK. */ 00192 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32)); 00193 memcpy(tmp_a[0], a[0], n*n*sizeof(float32)); 00194 00195 uplo = 'L'; 00196 spotrf_(&uplo, &n, tmp_a[0], &n, &info); 00197 det = tmp_a[0][0]; 00198 /* det = prod(diag(l))^2 */ 00199 for (i = 1; i < n; ++i) 00200 det *= tmp_a[i][i]; 00201 ckd_free_2d((void **)tmp_a); 00202 if (info > 0) 00203 return -1.0; /* Generic "not positive-definite" answer */ 00204 else 00205 return det * det; 00206 } 00207 00208 int32 00209 solve(float32 **a, /*Input : an n*n matrix A */ 00210 float32 *b, /*Input : a n dimesion vector b */ 00211 float32 *out_x, /*Output : a n dimesion vector x */ 00212 int32 n) 00213 { 00214 char uplo; 00215 float32 **tmp_a; 00216 int32 info, nrhs; 00217 00218 /* a is assumed to be symmetric, so we don't need to switch the 00219 * ordering of the data. But we do need to copy it since it is 00220 * overwritten by LAPACK. */ 00221 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32)); 00222 memcpy(tmp_a[0], a[0], n*n*sizeof(float32)); 00223 memcpy(out_x, b, n*sizeof(float32)); 00224 uplo = 'L'; 00225 nrhs = 1; 00226 sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, out_x, &n, &info); 00227 ckd_free_2d((void **)tmp_a); 00228 00229 if (info != 0) 00230 return -1; 00231 else 00232 return info; 00233 } 00234 00235 /* Find inverse by solving AX=I. */ 00236 int32 00237 invert(float32 ** ainv, float32 ** a, int32 n) 00238 { 00239 char uplo; 00240 float32 **tmp_a; 00241 int32 info, nrhs, i; 00242 00243 /* Construct an identity matrix. */ 00244 memset(ainv[0], 0, sizeof(float32) * n * n); 00245 for (i = 0; i < n; i++) 00246 ainv[i][i] = 1.0; 00247 /* a is assumed to be symmetric, so we don't need to switch the 00248 * ordering of the data. But we do need to copy it since it is 00249 * overwritten by LAPACK. */ 00250 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32)); 00251 memcpy(tmp_a[0], a[0], n*n*sizeof(float32)); 00252 uplo = 'L'; 00253 nrhs = n; 00254 sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, ainv[0], &n, &info); 00255 ckd_free_2d((void **)tmp_a); 00256 00257 if (info != 0) 00258 return -1; 00259 else 00260 return info; 00261 } 00262 00263 void 00264 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n) 00265 { 00266 char side, uplo; 00267 float32 alpha; 00268 00269 side = 'L'; 00270 uplo = 'L'; 00271 alpha = 1.0; 00272 ssymm_(&side, &uplo, &n, &n, &alpha, a[0], &n, b[0], &n, &alpha, c[0], &n); 00273 } 00274 00275 #endif /* WITH_LAPACK */ 00276 00277 void 00278 outerproduct(float32 ** a, float32 * x, float32 * y, int32 len) 00279 { 00280 int32 i, j; 00281 00282 for (i = 0; i < len; ++i) { 00283 a[i][i] = x[i] * y[i]; 00284 for (j = i + 1; j < len; ++j) { 00285 a[i][j] = x[i] * y[j]; 00286 a[j][i] = x[j] * y[i]; 00287 } 00288 } 00289 } 00290 00291 void 00292 scalarmultiply(float32 ** a, float32 x, int32 n) 00293 { 00294 int32 i, j; 00295 00296 for (i = 0; i < n; ++i) { 00297 a[i][i] *= x; 00298 for (j = i+1; j < n; ++j) { 00299 a[i][j] *= x; 00300 a[j][i] *= x; 00301 } 00302 } 00303 } 00304 00305 void 00306 matrixadd(float32 ** a, float32 ** b, int32 n) 00307 { 00308 int32 i, j; 00309 00310 for (i = 0; i < n; ++i) 00311 for (j = 0; j < n; ++j) 00312 a[i][j] += b[i][j]; 00313 } 00314 00315 00316 /* 00317 * Log record. Maintained by RCS. 00318 * 00319 * $Log$ 00320 * Revision 1.4 2004/07/21 18:05:40 egouvea 00321 * Changed the license terms to make it the same as sphinx2 and sphinx3. 00322 * 00323 * Revision 1.3 2001/04/05 20:02:30 awb 00324 * *** empty log message *** 00325 * 00326 * Revision 1.2 2000/09/29 22:35:13 awb 00327 * *** empty log message *** 00328 * 00329 * Revision 1.1 2000/09/24 21:38:31 awb 00330 * *** empty log message *** 00331 * 00332 * Revision 1.1 97/07/16 11:36:22 eht 00333 * Initial revision 00334 * 00335 * 00336 */