SphinxBase  0.6
src/libsphinxbase/util/matrix.c
00001 /* -*- c-basic-offset: 4 -*- */
00002 /* ====================================================================
00003  * Copyright (c) 1997-2006 Carnegie Mellon University.  All rights 
00004  * reserved.
00005  *
00006  * Redistribution and use in source and binary forms, with or without
00007  * modification, are permitted provided that the following conditions
00008  * are met:
00009  *
00010  * 1. Redistributions of source code must retain the above copyright
00011  *    notice, this list of conditions and the following disclaimer. 
00012  *
00013  * 2. Redistributions in binary form must reproduce the above copyright
00014  *    notice, this list of conditions and the following disclaimer in
00015  *    the documentation and/or other materials provided with the
00016  *    distribution.
00017  *
00018  * This work was supported in part by funding from the Defense Advanced 
00019  * Research Projects Agency and the National Science Foundation of the 
00020  * United States of America, and the CMU Sphinx Speech Consortium.
00021  *
00022  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
00023  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
00024  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00025  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
00026  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00027  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
00028  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
00029  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
00030  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00031  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
00032  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00033  *
00034  * ====================================================================
00035  *
00036  */
00037 #include <string.h>
00038 #include <stdlib.h>
00039 
00040 #ifdef HAVE_CONFIG_H
00041 #include "config.h"
00042 #endif
00043 
00044 #include "sphinxbase/clapack_lite.h"
00045 #include "sphinxbase/matrix.h"
00046 #include "sphinxbase/err.h"
00047 #include "sphinxbase/ckd_alloc.h"
00048 
00049 void
00050 norm_3d(float32 ***arr,
00051         uint32 d1,
00052         uint32 d2,
00053         uint32 d3)
00054 {
00055     uint32 i, j, k;
00056     float64 s;
00057 
00058     for (i = 0; i < d1; i++) {
00059         for (j = 0; j < d2; j++) {
00060 
00061             /* compute sum (i, j) as over all k */
00062             for (k = 0, s = 0; k < d3; k++) {
00063                 s += arr[i][j][k];
00064             }
00065 
00066             /* do 1 floating point divide */
00067             s = 1.0 / s;
00068 
00069             /* divide all k by sum over k */
00070             for (k = 0; k < d3; k++) {
00071                 arr[i][j][k] *= s;
00072             }
00073         }
00074     }
00075 }
00076 
00077 void
00078 accum_3d(float32 ***out,
00079          float32 ***in,
00080          uint32 d1,
00081          uint32 d2,
00082          uint32 d3)
00083 {
00084     uint32 i, j, k;
00085 
00086     for (i = 0; i < d1; i++) {
00087         for (j = 0; j < d2; j++) {
00088             for (k = 0; k < d3; k++) {
00089                 out[i][j][k] += in[i][j][k];
00090             }
00091         }
00092     }
00093 }
00094 
00095 void
00096 floor_nz_3d(float32 ***m,
00097             uint32 d1,
00098             uint32 d2,
00099             uint32 d3,
00100             float32 floor)
00101 {
00102     uint32 i, j, k;
00103 
00104     for (i = 0; i < d1; i++) {
00105         for (j = 0; j < d2; j++) {
00106             for (k = 0; k < d3; k++) {
00107                 if ((m[i][j][k] != 0) && (m[i][j][k] < floor))
00108                     m[i][j][k] = floor;
00109             }
00110         }
00111     }
00112 }
00113 void
00114 floor_nz_1d(float32 *v,
00115             uint32 d1,
00116             float32 floor)
00117 {
00118     uint32 i;
00119 
00120     for (i = 0; i < d1; i++) {
00121         if ((v[i] != 0) && (v[i] < floor))
00122             v[i] = floor;
00123     }
00124 }
00125 
00126 void
00127 band_nz_1d(float32 *v,
00128            uint32 d1,
00129            float32 band)
00130 {
00131     uint32 i;
00132 
00133     for (i = 0; i < d1; i++) {
00134         if (v[i] != 0) {
00135             if ((v[i] > 0) && (v[i] < band)) {
00136                 v[i] = band;
00137             }
00138             else if ((v[i] < 0) && (v[i] > -band)) {
00139                 v[i] = -band;
00140             }
00141         }
00142     }
00143 }
00144 
00145 #ifndef WITH_LAPACK
00146 float64
00147 determinant(float32 **a, int32 n)
00148 {
00149     E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n");
00150     return 0.0;
00151 }
00152 int32
00153 invert(float32 **ainv, float32 **a, int32 n)
00154 {
00155     E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n");
00156     return 0;
00157 }
00158 int32
00159 solve(float32 **a, float32 *b, float32 *out_x, int32   n)
00160 {
00161     E_FATAL("No LAPACK library available, cannot solve linear equations (FIXME)\n");
00162     return 0;
00163 }
00164 
00165 void
00166 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
00167 {
00168     int32 i, j, k;
00169 
00170     memset(c[0], 0, n*n*sizeof(float32));
00171     for (i = 0; i < n; ++i) {
00172         for (j = 0; j < n; ++j) {
00173             for (k = 0; k < n; ++k) {
00174                 c[i][k] += a[i][j] * b[j][k];
00175             }
00176         }
00177     }
00178 }
00179 #else /* WITH_LAPACK */
00180 /* Find determinant through LU decomposition. */
00181 float64
00182 determinant(float32 ** a, int32 n)
00183 {
00184     float32 **tmp_a;
00185     float64 det;
00186     char uplo;
00187     int32 info, i;
00188 
00189     /* a is assumed to be symmetric, so we don't need to switch the
00190      * ordering of the data.  But we do need to copy it since it is
00191      * overwritten by LAPACK. */
00192     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00193     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00194 
00195     uplo = 'L';
00196     spotrf_(&uplo, &n, tmp_a[0], &n, &info);
00197     det = tmp_a[0][0];
00198     /* det = prod(diag(l))^2 */
00199     for (i = 1; i < n; ++i)
00200         det *= tmp_a[i][i];
00201     ckd_free_2d((void **)tmp_a);
00202     if (info > 0)
00203         return -1.0; /* Generic "not positive-definite" answer */
00204     else
00205         return det * det;
00206 }
00207 
00208 int32
00209 solve(float32 **a, /*Input : an n*n matrix A */
00210       float32 *b,  /*Input : a n dimesion vector b */
00211       float32 *out_x,  /*Output : a n dimesion vector x */
00212       int32   n)
00213 {
00214     char uplo;
00215     float32 **tmp_a;
00216     int32 info, nrhs;
00217 
00218     /* a is assumed to be symmetric, so we don't need to switch the
00219      * ordering of the data.  But we do need to copy it since it is
00220      * overwritten by LAPACK. */
00221     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00222     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00223     memcpy(out_x, b, n*sizeof(float32));
00224     uplo = 'L';
00225     nrhs = 1;
00226     sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, out_x, &n, &info);
00227     ckd_free_2d((void **)tmp_a);
00228 
00229     if (info != 0)
00230         return -1;
00231     else
00232         return info;
00233 }
00234 
00235 /* Find inverse by solving AX=I. */
00236 int32
00237 invert(float32 ** ainv, float32 ** a, int32 n)
00238 {
00239     char uplo;
00240     float32 **tmp_a;
00241     int32 info, nrhs, i;
00242 
00243     /* Construct an identity matrix. */
00244     memset(ainv[0], 0, sizeof(float32) * n * n);
00245     for (i = 0; i < n; i++)
00246         ainv[i][i] = 1.0;
00247     /* a is assumed to be symmetric, so we don't need to switch the
00248      * ordering of the data.  But we do need to copy it since it is
00249      * overwritten by LAPACK. */
00250     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00251     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00252     uplo = 'L';
00253     nrhs = n;
00254     sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, ainv[0], &n, &info);
00255     ckd_free_2d((void **)tmp_a);
00256 
00257     if (info != 0)
00258         return -1;
00259     else
00260         return info;
00261 }
00262 
00263 void
00264 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
00265 {
00266     char side, uplo;
00267     float32 alpha;
00268 
00269     side = 'L';
00270     uplo = 'L';
00271     alpha = 1.0;
00272     ssymm_(&side, &uplo, &n, &n, &alpha, a[0], &n, b[0], &n, &alpha, c[0], &n);
00273 }
00274 
00275 #endif /* WITH_LAPACK */
00276 
00277 void
00278 outerproduct(float32 ** a, float32 * x, float32 * y, int32 len)
00279 {
00280     int32 i, j;
00281 
00282     for (i = 0; i < len; ++i) {
00283         a[i][i] = x[i] * y[i];
00284         for (j = i + 1; j < len; ++j) {
00285             a[i][j] = x[i] * y[j];
00286             a[j][i] = x[j] * y[i];
00287         }
00288     }
00289 }
00290 
00291 void
00292 scalarmultiply(float32 ** a, float32 x, int32 n)
00293 {
00294     int32 i, j;
00295 
00296     for (i = 0; i < n; ++i) {
00297         a[i][i] *= x;
00298         for (j = i+1; j < n; ++j) {
00299             a[i][j] *= x;
00300             a[j][i] *= x;
00301         }
00302     }
00303 }
00304 
00305 void
00306 matrixadd(float32 ** a, float32 ** b, int32 n)
00307 {
00308     int32 i, j;
00309 
00310     for (i = 0; i < n; ++i)
00311         for (j = 0; j < n; ++j)
00312             a[i][j] += b[i][j];
00313 }
00314 
00315 
00316 /*
00317  * Log record.  Maintained by RCS.
00318  *
00319  * $Log$
00320  * Revision 1.4  2004/07/21  18:05:40  egouvea
00321  * Changed the license terms to make it the same as sphinx2 and sphinx3.
00322  * 
00323  * Revision 1.3  2001/04/05 20:02:30  awb
00324  * *** empty log message ***
00325  *
00326  * Revision 1.2  2000/09/29 22:35:13  awb
00327  * *** empty log message ***
00328  *
00329  * Revision 1.1  2000/09/24 21:38:31  awb
00330  * *** empty log message ***
00331  *
00332  * Revision 1.1  97/07/16  11:36:22  eht
00333  * Initial revision
00334  * 
00335  *
00336  */