14a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project/*---------------------------------------------------------------------------* 24a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * matrix_i.c * 34a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * * 44a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * Copyright 2007, 2008 Nuance Communciations, Inc. * 54a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * * 64a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * Licensed under the Apache License, Version 2.0 (the 'License'); * 74a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * you may not use this file except in compliance with the License. * 84a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * * 94a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * You may obtain a copy of the License at * 104a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * http://www.apache.org/licenses/LICENSE-2.0 * 114a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * * 124a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * Unless required by applicable law or agreed to in writing, software * 134a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * distributed under the License is distributed on an 'AS IS' BASIS, * 144a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * 154a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * See the License for the specific language governing permissions and * 164a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * limitations under the License. * 174a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project * * 184a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project *---------------------------------------------------------------------------*/ 194a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 204a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 214a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include <stdlib.h> 224a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#ifndef _RTT 234a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include <stdio.h> 244a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#endif 254a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include <math.h> 264a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include <string.h> 274a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include <assert.h> 284a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 294a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include "hmmlib.h" 304a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include "prelib.h" 314a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#include "portable.h" 324a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 334a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#define PIVOT 1 344a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#define DEBUG 0 354a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 364a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#define TINY 1.0e-20 374a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#define SIGNIFICANT 0 /* 1.0e-20 */ 384a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 394a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Projectstatic const char matrix_i[] = "$Id: matrix_i.c,v 1.2.10.3 2007/10/15 18:06:24 dahan Exp $"; 404a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 414a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Projectvoid lubksb(double **mat, int dim, int *index, double *b); 424a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Projectint ludcmp(double **mat, int dim, int *index); 434a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 444a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Projectint invert_matrix(covdata **mat, covdata **inv, int dim) 454a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project{ 464a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project double *col, **input; 474a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project int ii, jj, *index, err_code; 484a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#if DEBUG 494a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project double sum; 504a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project int kk; 514a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#endif 524a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 534a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project ASSERT(mat); 544a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project ASSERT(inv); 554a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project index = (int *) CALLOC(dim, sizeof(int), "clib.index_imatrix"); 564a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project col = (double *) CALLOC(dim, sizeof(double), "clib.col"); 574a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project input = (double **) CALLOC(dim, sizeof(double *), "clib.input_imatrix"); 584a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 594a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 604a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project input[ii] = (double *) CALLOC(dim, sizeof(double), "clib.input_imatrix[]"); 614a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = 0; jj < dim; jj++) 624a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project input[ii][jj] = mat[ii][jj]; 634a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 644a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 654a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if ((err_code = ludcmp(input, dim, index)) > 0) return(err_code); 664a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = 0; jj < dim; jj++) 674a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 684a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 694a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project col[ii] = 0; 704a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project col[jj] = 1; 714a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project lubksb(input, dim, index, col); 724a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 734a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project inv[ii][jj] = col[ii]; 744a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 754a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 764a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project FREE((char *)input[ii]); 774a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project FREE((char *)input); 784a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project FREE((char *)col); 794a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project FREE((char *)index); 804a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 814a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#if DEBUG 824a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project printf("Testing the inverse:\n"); 834a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 844a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 854a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = 0; jj < dim; jj++) 864a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 874a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum = 0; 884a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (kk = 0; kk < dim; kk++) 894a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum += mat[ii][kk] * inv[kk][jj]; 904a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project printf("%.2f ", sum); 914a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 924a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project printf("\n"); 934a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 944a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#endif 954a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 964a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project return (0); 974a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project} 984a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 994a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Projectint ludcmp(double **mat, int dim, int *index) 1004a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project/* 1014a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project** This routine is straight out of the numerical recipes in C 1024a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project*/ 1034a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project{ 1044a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project int ii, imax = 0, jj, kk; 1054a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project double big, dumm, sum, temp; 1064a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project double *vv; 1074a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1084a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project vv = (double *) CALLOC(dim + 5, sizeof(double), "clib.ludcmp.vv"); 1094a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#if PIVOT 1104a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 1114a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1124a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project big = 0; 1134a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = 0; jj < dim; jj++) 1144a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if ((temp = (double) fabs(mat[ii][jj])) > big) big = temp; 1154a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if (big <= SIGNIFICANT) 1164a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1174a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project log_report("Singular matrix in routine ludcmp\n"); 1184a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project return (SINGULAR_MATRIX); 1194a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1204a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project vv[ii] = 1 / big; 1214a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1224a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#endif 1234a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1244a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = 0; jj < dim; jj++) 1254a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1264a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < jj; ii++) 1274a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1284a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum = mat[ii][jj]; 1294a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (kk = 0; kk < ii; kk++) 1304a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum -= mat[ii][kk] * mat[kk][jj]; 1314a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project mat[ii][jj] = sum; 1324a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1334a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project big = 0; 1344a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = jj; ii < dim; ii++) 1354a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1364a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum = mat[ii][jj]; 1374a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (kk = 0; kk < jj; kk++) 1384a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum -= mat[ii][kk] * mat[kk][jj]; 1394a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project mat[ii][jj] = sum; 1404a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if ((dumm = (double)(vv[ii] * fabs(sum))) >= big) 1414a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1424a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project big = dumm; 1434a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project imax = ii; 1444a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1454a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1464a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1474a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#if PIVOT 1484a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if (jj != imax) 1494a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1504a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (kk = 0; kk < dim; kk++) 1514a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1524a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project dumm = mat[imax][kk]; 1534a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project mat[imax][kk] = mat[jj][kk]; 1544a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project mat[jj][kk] = dumm; 1554a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1564a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project vv[imax] = vv[jj]; 1574a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1584a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project index[jj] = imax; 1594a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#else 1604a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project index[jj] = jj; 1614a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#endif 1624a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1634a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if (fabs(mat[jj][jj]) <= SIGNIFICANT) mat[jj][jj] = (double)TINY; 1644a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if (jj < (dim - 1)) 1654a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1664a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project dumm = 1 / mat[jj][jj]; 1674a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = jj + 1; ii < dim; ii++) 1684a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project mat[ii][jj] *= dumm; 1694a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1704a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1714a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1724a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project FREE((char *)vv); 1734a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project return (0); 1744a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project} 1754a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1764a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Projectvoid lubksb(double **mat, int dim, int *index, double *b) 1774a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project{ 1784a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project int ii, ix = -1, ip, jj; 1794a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project double sum; 1804a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1814a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = 0; ii < dim; ii++) 1824a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1834a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#if PIVOT 1844a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project ip = index[ii]; 1854a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum = b[ip]; 1864a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project b[ip] = b[ii]; 1874a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#else 1884a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum = b[ii]; 1894a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project#endif 1904a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project if (ix >= 0) 1914a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = ix; jj < ii; jj++) 1924a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum -= mat[ii][jj] * b[jj]; 1934a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project else if (sum) ix = ii; 1944a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project b[ii] = sum; 1954a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 1964a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project 1974a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (ii = dim - 1; ii >= 0; ii--) 1984a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project { 1994a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum = b[ii]; 2004a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project for (jj = ii + 1; jj < dim; jj++) 2014a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project sum -= mat[ii][jj] * b[jj]; 2024a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project b[ii] = sum / mat[ii][ii]; 2034a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project } 2044a68b3365c8c50aa93505e99ead2565ab73dcdb0The Android Open Source Project} 205