src/libsphinxbase/util/matrix.c

   1 /* -*- c-basic-offset: 4 -*- */
   2 /* ====================================================================
   3  * Copyright (c) 1997-2006 Carnegie Mellon University.  All rights
   4  * reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions
   8  * are met:
   9  *
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  *
  13  * 2. Redistributions in binary form must reproduce the above copyright
  14  *    notice, this list of conditions and the following disclaimer in
  15  *    the documentation and/or other materials provided with the
  16  *    distribution.
  17  *
  18  * This work was supported in part by funding from the Defense Advanced
  19  * Research Projects Agency and the National Science Foundation of the
  20  * United States of America, and the CMU Sphinx Speech Consortium.
  21  *
  22  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
  23  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
  24  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  25  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
  26  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  27  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  28  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  29  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  30  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  31  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  32  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  33  *
  34  * ====================================================================
  35  *
  36  */
  37 #include <string.h>
  38 #include <stdlib.h>
  39
  40 #ifdef HAVE_CONFIG_H
  41 #include "config.h"
  42 #endif
  43
  44 #include "sphinxbase/clapack_lite.h"
  45 #include "sphinxbase/matrix.h"
  46 #include "sphinxbase/err.h"
  47 #include "sphinxbase/ckd_alloc.h"
  48
  49 void
  50 norm_3d(float32 ***arr,
  51         uint32 d1,
  52         uint32 d2,
  53         uint32 d3)
  54 {
  55     uint32 i, j, k;
  56     float64 s;
  57
  58     for (i = 0; i < d1; i++) {
  59         for (j = 0; j < d2; j++) {
  60
  61             /* compute sum (i, j) as over all k */
  62             for (k = 0, s = 0; k < d3; k++) {
  63                 s += arr[i][j][k];
  64             }
  65
  66             /* do 1 floating point divide */
  67             s = 1.0 / s;
  68
  69             /* divide all k by sum over k */
  70             for (k = 0; k < d3; k++) {
  71                 arr[i][j][k] *= s;
  72             }
  73         }
  74     }
  75 }
  76 \f
  77 void
  78 accum_3d(float32 ***out,
  79          float32 ***in,
  80          uint32 d1,
  81          uint32 d2,
  82          uint32 d3)
  83 {
  84     uint32 i, j, k;
  85
  86     for (i = 0; i < d1; i++) {
  87         for (j = 0; j < d2; j++) {
  88             for (k = 0; k < d3; k++) {
  89                 out[i][j][k] += in[i][j][k];
  90             }
  91         }
  92     }
  93 }
  94
  95 void
  96 floor_nz_3d(float32 ***m,
  97             uint32 d1,
  98             uint32 d2,
  99             uint32 d3,
 100             float32 floor)
 101 {
 102     uint32 i, j, k;
 103
 104     for (i = 0; i < d1; i++) {
 105         for (j = 0; j < d2; j++) {
 106             for (k = 0; k < d3; k++) {
 107                 if ((m[i][j][k] != 0) && (m[i][j][k] < floor))
 108                     m[i][j][k] = floor;
 109             }
 110         }
 111     }
 112 }
 113 void
 114 floor_nz_1d(float32 *v,
 115             uint32 d1,
 116             float32 floor)
 117 {
 118     uint32 i;
 119
 120     for (i = 0; i < d1; i++) {
 121         if ((v[i] != 0) && (v[i] < floor))
 122             v[i] = floor;
 123     }
 124 }
 125
 126 void
 127 band_nz_1d(float32 *v,
 128            uint32 d1,
 129            float32 band)
 130 {
 131     uint32 i;
 132
 133     for (i = 0; i < d1; i++) {
 134         if (v[i] != 0) {
 135             if ((v[i] > 0) && (v[i] < band)) {
 136                 v[i] = band;
 137             }
 138             else if ((v[i] < 0) && (v[i] > -band)) {
 139                 v[i] = -band;
 140             }
 141         }
 142     }
 143 }
 144
 145 #ifndef WITH_LAPACK
 146 float64
 147 determinant(float32 **a, int32 n)
 148 {
 149     E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n");
 150     return 0.0;
 151 }
 152 int32
 153 invert(float32 **ainv, float32 **a, int32 n)
 154 {
 155     E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n");
 156     return 0;
 157 }
 158 int32
 159 solve(float32 **a, float32 *b, float32 *out_x, int32   n)
 160 {
 161     E_FATAL("No LAPACK library available, cannot solve linear equations (FIXME)\n");
 162     return 0;
 163 }
 164
 165 void
 166 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
 167 {
 168     int32 i, j, k;
 169
 170     memset(c[0], 0, n*n*sizeof(float32));
 171     for (i = 0; i < n; ++i) {
 172         for (j = 0; j < n; ++j) {
 173             for (k = 0; k < n; ++k) {
 174                 c[i][k] += a[i][j] * b[j][k];
 175             }
 176         }
 177     }
 178 }
 179 #else /* WITH_LAPACK */
 180 /* Find determinant through LU decomposition. */
 181 float64
 182 determinant(float32 ** a, int32 n)
 183 {
 184     float32 **tmp_a;
 185     float64 det;
 186     char uplo;
 187     int32 info, i;
 188
 189     /* a is assumed to be symmetric, so we don't need to switch the
 190      * ordering of the data.  But we do need to copy it since it is
 191      * overwritten by LAPACK. */
 192     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
 193     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
 194
 195     uplo = 'L';
 196     spotrf_(&uplo, &n, tmp_a[0], &n, &info);
 197     det = tmp_a[0][0];
 198     /* det = prod(diag(l))^2 */
 199     for (i = 1; i < n; ++i)
 200         det *= tmp_a[i][i];
 201     ckd_free_2d((void **)tmp_a);
 202     if (info > 0)
 203         return -1.0; /* Generic "not positive-definite" answer */
 204     else
 205         return det * det;
 206 }
 207
 208 int32
 209 solve(float32 **a, /*Input : an n*n matrix A */
 210       float32 *b,  /*Input : a n dimesion vector b */
 211       float32 *out_x,  /*Output : a n dimesion vector x */
 212       int32   n)
 213 {
 214     char uplo;
 215     float32 **tmp_a;
 216     int32 info, nrhs;
 217
 218     /* a is assumed to be symmetric, so we don't need to switch the
 219      * ordering of the data.  But we do need to copy it since it is
 220      * overwritten by LAPACK. */
 221     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
 222     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
 223     memcpy(out_x, b, n*sizeof(float32));
 224     uplo = 'L';
 225     nrhs = 1;
 226     sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, out_x, &n, &info);
 227     ckd_free_2d((void **)tmp_a);
 228
 229     if (info != 0)
 230         return -1;
 231     else
 232         return info;
 233 }
 234
 235 /* Find inverse by solving AX=I. */
 236 int32
 237 invert(float32 ** ainv, float32 ** a, int32 n)
 238 {
 239     char uplo;
 240     float32 **tmp_a;
 241     int32 info, nrhs, i;
 242
 243     /* Construct an identity matrix. */
 244     memset(ainv[0], 0, sizeof(float32) * n * n);
 245     for (i = 0; i < n; i++)
 246         ainv[i][i] = 1.0;
 247     /* a is assumed to be symmetric, so we don't need to switch the
 248      * ordering of the data.  But we do need to copy it since it is
 249      * overwritten by LAPACK. */
 250     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
 251     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
 252     uplo = 'L';
 253     nrhs = n;
 254     sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, ainv[0], &n, &info);
 255     ckd_free_2d((void **)tmp_a);
 256
 257     if (info != 0)
 258         return -1;
 259     else
 260         return info;
 261 }
 262
 263 void
 264 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
 265 {
 266     char side, uplo;
 267     float32 alpha;
 268
 269     side = 'L';
 270     uplo = 'L';
 271     alpha = 1.0;
 272     ssymm_(&side, &uplo, &n, &n, &alpha, a[0], &n, b[0], &n, &alpha, c[0], &n);
 273 }
 274
 275 #endif /* WITH_LAPACK */
 276
 277 void
 278 outerproduct(float32 ** a, float32 * x, float32 * y, int32 len)
 279 {
 280     int32 i, j;
 281
 282     for (i = 0; i < len; ++i) {
 283         a[i][i] = x[i] * y[i];
 284         for (j = i + 1; j < len; ++j) {
 285             a[i][j] = x[i] * y[j];
 286             a[j][i] = x[j] * y[i];
 287         }
 288     }
 289 }
 290
 291 void
 292 scalarmultiply(float32 ** a, float32 x, int32 n)
 293 {
 294     int32 i, j;
 295
 296     for (i = 0; i < n; ++i) {
 297         a[i][i] *= x;
 298         for (j = i+1; j < n; ++j) {
 299             a[i][j] *= x;
 300             a[j][i] *= x;
 301         }
 302     }
 303 }
 304
 305 void
 306 matrixadd(float32 ** a, float32 ** b, int32 n)
 307 {
 308     int32 i, j;
 309
 310     for (i = 0; i < n; ++i)
 311         for (j = 0; j < n; ++j)
 312             a[i][j] += b[i][j];
 313 }
 314
 315
 316 /*
 317  * Log record.  Maintained by RCS.
 318  *
 319  * $Log$
 320  * Revision 1.4  2004/07/21  18:05:40  egouvea
 321  * Changed the license terms to make it the same as sphinx2 and sphinx3.
 322  *
 323  * Revision 1.3  2001/04/05 20:02:30  awb
 324  * *** empty log message ***
 325  *
 326  * Revision 1.2  2000/09/29 22:35:13  awb
 327  * *** empty log message ***
 328  *
 329  * Revision 1.1  2000/09/24 21:38:31  awb
 330  * *** empty log message ***
 331  *
 332  * Revision 1.1  97/07/16  11:36:22  eht
 333  * Initial revision
 334  *
 335  *
 336  */