1/*---------------------------------------------------------------------------*
2 *  jacobi.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#ifndef _RTT
22#include <stdio.h>
23#endif
24#include <stdlib.h>
25#include <math.h>
26#include <string.h>
27#include <assert.h>
28
29#include "hmmlib.h"
30#include "portable.h"
31
32#undef EPSILON
33#define EPSILON         0.00000001
34#define MAX_ITER      100
35#define DEFAULT_MAX_PC_DIFF     0.0001
36#define DEFAULT_MAX_OFF_DIAG 0.0001
37
38static void   Rotate(double **a, int dim, int i, int j, int k, int l, double s,
39                     double tau);
40
41static double SumOffDiag(double **mat, int dim);
42
43
44/* ------------------------------------------------------------------------ */
45
46void Jacobi(double **matrix, int dim, double *egval, double **egvec)
47/*
48//    Compute all eigenvalues and eigenvectors of the real symmetric matrix
49//    <mat>  of size  <dim>x<dim>.  Fills in <egval> with the eigenvalues
50//    and  <egvec[i]>  with the i-th normalized eignevector of  <mat>.
51*/
52{
53  int    i, j, k;
54  int    nRotations,  /* number of Jacobi rotations that are performed */
55  iter;
56  double g, thresh, sum, c, s, t, tau, h;
57  double *b, *d, *z;
58  double **v, **a;
59
60  ASSERT(matrix);
61  ASSERT(egval);
62  ASSERT(egvec);
63
64  b = (double *) CALLOC(dim, sizeof(double), "clib.jacobi.b");
65  d = (double *) CALLOC(dim, sizeof(double), "clib.jacobi.d");
66  z = (double *) CALLOC(dim, sizeof(double), "clib.jacobi.z");
67  a = (double **) CALLOC(dim, sizeof(double *), "clib.jacobi.input_jacobi");
68  v = (double **) CALLOC(dim, sizeof(double *), "clib.jacobi.input_jacobi");
69  for (i = 0; i < dim; i++)
70  {
71    a[i] = (double *) CALLOC(dim, sizeof(double), "clib.jacobi.input_jacobi[]");
72    v[i] = (double *) CALLOC(dim, sizeof(double), "clib.jacobi.input_jacobi[]");
73    for (j = 0; j < dim; j++)
74      a[i][j] = (float) matrix[i][j];
75  }
76
77  /* initialize v to identity matrix, d and b to the diagonal of mat */
78  for (i = 0; i < dim; i++)
79  {
80    v[i][i] = 1.0;
81    b[i] = d[i] = a[i][i];
82  }
83
84  nRotations = 0;
85  iter = 0;
86
87  while (1)
88  {
89    sum = SumOffDiag(a, dim);
90
91    if (sum < EPSILON)    /* normal convergence */
92    {
93      log_report("\nConverged after %u iterations", iter);
94      break;
95    }
96    if (iter >= MAX_ITER)
97    {
98      log_report("\nMax number %u of iterations reached",
99                 MAX_ITER);
100      break;
101    }
102
103    if (iter < 3)
104      thresh = 20.0 * sum / (dim * dim);    /* .. first 3 iterations only */
105    else
106      thresh = 0.0;                    /* .. thereafter */
107
108    for (i = 0; i < dim - 1; i++)
109    {
110      for (j = i + 1; j < dim; j++)
111      {
112        g = 100.0 * fabs(a[i][i]);
113
114        /* after 4 iter's, skip rotation if off-diag elmt is small */
115        if ((iter >= 4) && (g < EPSILON*fabs(d[i]))
116            && (g < EPSILON*fabs(d[j])))
117        {
118          a[i][i] = 0.0;
119        }
120        else if (g > thresh)
121        {
122          h = d[j] - d[i];
123
124          if (g < EPSILON*fabs(h))
125            t = a[i][j] / h;
126          else
127          {
128            double theta = 0.5 * h / a[i][j];
129            t = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta));
130            if (theta < 0.0)  t = -t;
131          }
132
133          c = 1.0 / sqrt(1 + t * t);
134          s = t * c;
135          tau = s / (1.0 + c);
136          h = t * a[i][j];
137
138          z[i] -= h;
139          z[j] += h;
140          d[i] -= h;
141          d[j] += h;
142          a[i][j] = 0.0;
143
144          for (k = 0  ; k < i;   k++)  Rotate(a, dim, k, i, k, j, s, tau);
145          for (k = i + 1; k < j;   k++)  Rotate(a, dim, i, k, k, j, s, tau);
146          for (k = j + 1; k < dim; k++)  Rotate(a, dim, i, k, j, k, s, tau);
147
148          for (k = 0;   k < dim; k++)  Rotate(v, dim, k, i, k, j, s, tau);
149
150          nRotations++;
151        }
152      }
153    }
154
155    for (i = 0; i < dim; i++)    /* update d[] and re-initialize z[] */
156    {
157      b[i] += z[i];
158      d[i] = b[i];
159      z[i] = 0.0;
160    }
161
162    iter++;
163  }
164
165  /* save eigenvectors and eigenvalues */
166  for (i = 0; i < dim; i++)
167  {
168    /* the i-th column of v is the i-th eigenvector */
169    for (j = 0; j < dim; j++)  egvec[i][j] = v[j][i]; /* TODO: should this be egvec[j][i] */
170
171    /* the i-th entry of d is the i-th eigenvalue */
172    egval[i] = d[i];
173  }
174
175  log_report("\nDiagonalization required %u Jacobi rotations",
176             nRotations);
177
178  FREE((char *)b);
179  FREE((char *)d);
180  FREE((char *)z);
181  for (i = 0; i < dim; i++)
182  {
183    FREE((char *)a[i]);
184    FREE((char *)v[i]);
185  }
186  FREE((char *)a);
187  FREE((char *)v);
188  return;
189}
190
191/* ------------------------------------------------------------------------ */
192
193static void Rotate(double **a, int dim, int i, int j, int k, int l, double s,
194                   double tau)
195{
196  double g = a[i][j],
197             h = a[k][l];
198
199  a[i][j] = g - s * (h + g * tau);
200  a[k][l] = h + s * (g - h * tau);
201  return;
202}
203
204static double SumOffDiag(double **mat, int dim)
205/*
206//    Compute the sum of the absolute values of the off-diagonal elements
207**/
208{
209  int    i, j;
210  double sum = 0.0;
211
212  for (i = 0; i < dim - 1; i++)
213    for (j = i + 1; j < dim; j++)
214      sum += fabs(mat[i][j]);
215
216  return (sum);
217}
218
219