• Main Page
  • Related Pages
  • Data Structures
  • Files
  • File List
  • Globals

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 "clapack_lite.h"
00045 #include "matrix.h"
00046 #include "err.h"
00047 #include "ckd_alloc.h"
00048 
00049 #ifndef WITH_LAPACK
00050 float64
00051 determinant(float32 **a, int32 n)
00052 {
00053     E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n");
00054     return 0.0;
00055 }
00056 int32
00057 invert(float32 **ainv, float32 **a, int32 n)
00058 {
00059     E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n");
00060     return 0;
00061 }
00062 int32
00063 solve(float32 **a, float32 *b, float32 *out_x, int32   n)
00064 {
00065     E_FATAL("No LAPACK library available, cannot solve linear equations (FIXME)\n");
00066     return 0;
00067 }
00068 
00069 void
00070 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
00071 {
00072     int32 i, j, k;
00073 
00074     memset(c[0], 0, n*n*sizeof(float32));
00075     for (i = 0; i < n; ++i) {
00076         for (j = 0; j < n; ++j) {
00077             for (k = 0; k < n; ++k) {
00078                 c[i][k] += a[i][j] * b[j][k];
00079             }
00080         }
00081     }
00082 }
00083 #else /* WITH_LAPACK */
00084 /* Find determinant through LU decomposition. */
00085 float64
00086 determinant(float32 ** a, int32 n)
00087 {
00088     float32 **tmp_a;
00089     float64 det;
00090     char uplo;
00091     int32 info, i;
00092 
00093     /* a is assumed to be symmetric, so we don't need to switch the
00094      * ordering of the data.  But we do need to copy it since it is
00095      * overwritten by LAPACK. */
00096     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00097     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00098 
00099     uplo = 'L';
00100     spotrf_(&uplo, &n, tmp_a[0], &n, &info);
00101     det = tmp_a[0][0];
00102     /* det = prod(diag(l))^2 */
00103     for (i = 1; i < n; ++i)
00104         det *= tmp_a[i][i];
00105     ckd_free_2d((void **)tmp_a);
00106     if (info > 0)
00107         return -1.0; /* Generic "not positive-definite" answer */
00108     else
00109         return det * det;
00110 }
00111 
00112 int32
00113 solve(float32 **a, /*Input : an n*n matrix A */
00114       float32 *b,  /*Input : a n dimesion vector b */
00115       float32 *out_x,  /*Output : a n dimesion vector x */
00116       int32   n)
00117 {
00118     char uplo;
00119     float32 **tmp_a;
00120     int32 info, nrhs;
00121 
00122     /* a is assumed to be symmetric, so we don't need to switch the
00123      * ordering of the data.  But we do need to copy it since it is
00124      * overwritten by LAPACK. */
00125     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00126     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00127     memcpy(out_x, b, n*sizeof(float32));
00128     uplo = 'L';
00129     nrhs = 1;
00130     sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, out_x, &n, &info);
00131     ckd_free_2d((void **)tmp_a);
00132 
00133     if (info != 0)
00134         return -1;
00135     else
00136         return info;
00137 }
00138 
00139 /* Find inverse by solving AX=I. */
00140 int32
00141 invert(float32 ** ainv, float32 ** a, int32 n)
00142 {
00143     char uplo;
00144     float32 **tmp_a;
00145     int32 info, nrhs, i;
00146 
00147     /* Construct an identity matrix. */
00148     memset(ainv[0], 0, sizeof(float32) * n * n);
00149     for (i = 0; i < n; i++)
00150         ainv[i][i] = 1.0;
00151     /* a is assumed to be symmetric, so we don't need to switch the
00152      * ordering of the data.  But we do need to copy it since it is
00153      * overwritten by LAPACK. */
00154     tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00155     memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00156     uplo = 'L';
00157     nrhs = n;
00158     sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, ainv[0], &n, &info);
00159     ckd_free_2d((void **)tmp_a);
00160 
00161     if (info != 0)
00162         return -1;
00163     else
00164         return info;
00165 }
00166 
00167 void
00168 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
00169 {
00170     char side, uplo;
00171     float32 alpha;
00172 
00173     side = 'L';
00174     uplo = 'L';
00175     alpha = 1.0;
00176     ssymm_(&side, &uplo, &n, &n, &alpha, a[0], &n, b[0], &n, &alpha, c[0], &n);
00177 }
00178 
00179 #endif /* WITH_LAPACK */
00180 
00181 void
00182 outerproduct(float32 ** a, float32 * x, float32 * y, int32 len)
00183 {
00184     int32 i, j;
00185 
00186     for (i = 0; i < len; ++i) {
00187         a[i][i] = x[i] * y[i];
00188         for (j = i + 1; j < len; ++j) {
00189             a[i][j] = x[i] * y[j];
00190             a[j][i] = x[j] * y[i];
00191         }
00192     }
00193 }
00194 
00195 void
00196 scalarmultiply(float32 ** a, float32 x, int32 n)
00197 {
00198     int32 i, j;
00199 
00200     for (i = 0; i < n; ++i) {
00201         a[i][i] *= x;
00202         for (j = i+1; j < n; ++j) {
00203             a[i][j] *= x;
00204             a[j][i] *= x;
00205         }
00206     }
00207 }
00208 
00209 void
00210 matrixadd(float32 ** a, float32 ** b, int32 n)
00211 {
00212     int32 i, j;
00213 
00214     for (i = 0; i < n; ++i)
00215         for (j = 0; j < n; ++j)
00216             a[i][j] += b[i][j];
00217 }
00218 
00219 
00220 /*
00221  * Log record.  Maintained by RCS.
00222  *
00223  * $Log$
00224  * Revision 1.4  2004/07/21  18:05:40  egouvea
00225  * Changed the license terms to make it the same as sphinx2 and sphinx3.
00226  * 
00227  * Revision 1.3  2001/04/05 20:02:30  awb
00228  * *** empty log message ***
00229  *
00230  * Revision 1.2  2000/09/29 22:35:13  awb
00231  * *** empty log message ***
00232  *
00233  * Revision 1.1  2000/09/24 21:38:31  awb
00234  * *** empty log message ***
00235  *
00236  * Revision 1.1  97/07/16  11:36:22  eht
00237  * Initial revision
00238  * 
00239  *
00240  */

Generated on Fri Jan 14 2011 for SphinxBase by  doxygen 1.7.1