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
39 static const char matrix_i[] = "$Id: matrix_i.c,v 1.2.10.3 2007/10/15 18:06:24 dahan Exp $";
40
41 void lubksb(double **mat, int dim, int *index, double *b);
42 int ludcmp(double **mat, int dim, int *index);
43
invert_matrix(covdata ** mat,covdata ** inv,int dim)44 int 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
ludcmp(double ** mat,int dim,int * index)99 int 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
lubksb(double ** mat,int dim,int * index,double * b)176 void 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