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 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 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 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
