1/*---------------------------------------------------------------------------*
2 *  matrix_i.c  *
3 *                                                                           *
4 *  Copyright 2007, 2008 Nuance Communciations, Inc.                               *
5 *                                                                           *
6 *  Licensed under the Apache License, Version 2.0 (the 'License');          *
7 *  you may not use this file except in compliance with the License.         *
8 *                                                                           *
9 *  You may obtain a copy of the License at                                  *
10 *      http://www.apache.org/licenses/LICENSE-2.0                           *
11 *                                                                           *
12 *  Unless required by applicable law or agreed to in writing, software      *
13 *  distributed under the License is distributed on an 'AS IS' BASIS,        *
14 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
15 *  See the License for the specific language governing permissions and      *
16 *  limitations under the License.                                           *
17 *                                                                           *
18 *---------------------------------------------------------------------------*/
19
20
21#include <stdlib.h>
22#ifndef _RTT
23#include <stdio.h>
24#endif
25#include <math.h>
26#include <string.h>
27#include <assert.h>
28
29#include "hmmlib.h"
30#include "prelib.h"
31#include "portable.h"
32
33#define PIVOT 1
34#define DEBUG 0
35
36#define TINY  1.0e-20
37#define SIGNIFICANT 0 /* 1.0e-20 */
38
39static const char matrix_i[] = "$Id: matrix_i.c,v 1.2.10.3 2007/10/15 18:06:24 dahan Exp $";
40
41void lubksb(double **mat, int dim, int *index, double *b);
42int  ludcmp(double **mat, int dim, int *index);
43
44int invert_matrix(covdata **mat, covdata **inv, int dim)
45{
46  double *col, **input;
47  int ii, jj, *index, err_code;
48#if DEBUG
49  double  sum;
50  int     kk;
51#endif
52
53  ASSERT(mat);
54  ASSERT(inv);
55  index = (int *) CALLOC(dim, sizeof(int), "clib.index_imatrix");
56  col = (double *) CALLOC(dim, sizeof(double), "clib.col");
57  input = (double **) CALLOC(dim, sizeof(double *), "clib.input_imatrix");
58  for (ii = 0; ii < dim; ii++)
59  {
60    input[ii] = (double *) CALLOC(dim, sizeof(double), "clib.input_imatrix[]");
61    for (jj = 0; jj < dim; jj++)
62      input[ii][jj] = mat[ii][jj];
63  }
64
65  if ((err_code = ludcmp(input, dim, index)) > 0) return(err_code);
66  for (jj = 0; jj < dim; jj++)
67  {
68    for (ii = 0; ii < dim; ii++)
69      col[ii] = 0;
70    col[jj] = 1;
71    lubksb(input, dim, index, col);
72    for (ii = 0; ii < dim; ii++)
73      inv[ii][jj] = col[ii];
74  }
75  for (ii = 0; ii < dim; ii++)
76    FREE((char *)input[ii]);
77  FREE((char *)input);
78  FREE((char *)col);
79  FREE((char *)index);
80
81#if DEBUG
82  printf("Testing the inverse:\n");
83  for (ii = 0; ii < dim; ii++)
84  {
85    for (jj = 0; jj < dim; jj++)
86    {
87      sum = 0;
88      for (kk = 0; kk < dim; kk++)
89        sum += mat[ii][kk] * inv[kk][jj];
90      printf("%.2f ", sum);
91    }
92    printf("\n");
93  }
94#endif
95
96  return (0);
97}
98
99int ludcmp(double **mat, int dim, int *index)
100/*
101**  This routine is straight out of the numerical recipes in C
102*/
103{
104  int ii, imax = 0, jj, kk;
105  double big, dumm, sum, temp;
106  double *vv;
107
108  vv = (double *) CALLOC(dim + 5, sizeof(double), "clib.ludcmp.vv");
109#if PIVOT
110  for (ii = 0; ii < dim; ii++)
111  {
112    big = 0;
113    for (jj = 0; jj < dim; jj++)
114      if ((temp = (double) fabs(mat[ii][jj])) > big) big = temp;
115    if (big <= SIGNIFICANT)
116    {
117      log_report("Singular matrix in routine ludcmp\n");
118      return (SINGULAR_MATRIX);
119    }
120    vv[ii] = 1 / big;
121  }
122#endif
123
124  for (jj = 0; jj < dim; jj++)
125  {
126    for (ii = 0; ii < jj; ii++)
127    {
128      sum = mat[ii][jj];
129      for (kk = 0; kk < ii; kk++)
130        sum -= mat[ii][kk] * mat[kk][jj];
131      mat[ii][jj] = sum;
132    }
133    big = 0;
134    for (ii = jj; ii < dim; ii++)
135    {
136      sum = mat[ii][jj];
137      for (kk = 0; kk < jj; kk++)
138        sum -= mat[ii][kk] * mat[kk][jj];
139      mat[ii][jj] = sum;
140      if ((dumm = (double)(vv[ii] * fabs(sum))) >= big)
141      {
142        big = dumm;
143        imax = ii;
144      }
145    }
146
147#if PIVOT
148    if (jj != imax)
149    {
150      for (kk = 0; kk < dim; kk++)
151      {
152        dumm = mat[imax][kk];
153        mat[imax][kk] = mat[jj][kk];
154        mat[jj][kk] = dumm;
155      }
156      vv[imax] = vv[jj];
157    }
158    index[jj] = imax;
159#else
160    index[jj] = jj;
161#endif
162
163    if (fabs(mat[jj][jj]) <= SIGNIFICANT) mat[jj][jj] = (double)TINY;
164    if (jj < (dim - 1))
165    {
166      dumm = 1 / mat[jj][jj];
167      for (ii = jj + 1; ii < dim; ii++)
168        mat[ii][jj] *= dumm;
169    }
170  }
171
172  FREE((char *)vv);
173  return (0);
174}
175
176void lubksb(double **mat, int dim, int *index, double *b)
177{
178  int ii, ix = -1, ip, jj;
179  double sum;
180
181  for (ii = 0; ii < dim; ii++)
182  {
183#if PIVOT
184    ip = index[ii];
185    sum = b[ip];
186    b[ip] = b[ii];
187#else
188    sum = b[ii];
189#endif
190    if (ix >= 0)
191      for (jj = ix; jj < ii; jj++)
192        sum -= mat[ii][jj] * b[jj];
193    else if (sum) ix = ii;
194    b[ii] = sum;
195  }
196
197  for (ii = dim - 1; ii >= 0; ii--)
198  {
199    sum = b[ii];
200    for (jj = ii + 1; jj < dim; jj++)
201      sum -= mat[ii][jj] * b[jj];
202    b[ii] = sum / mat[ii][ii];
203  }
204}
205