matrix.c

一个通用的隐性马尔可夫C代码库 开发环境:C语言 简要说明:这是一个通用的隐性马尔可夫C代码库

/*******************************************************************************
  author       : Bernhard Knab
  filename     : ghmm/ghmm/matrix.c
  created      : TIME: ?            DATE: ?
  $Id: matrix.c,v 1.7 2004/01/05 19:23:12 schliep Exp $

Copyright (C) 1998-2001, ZAIK/ZPR, Universit鋞 zu K鰈n

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


*******************************************************************************/

#include <math.h>
#include <float.h>
#include "matrix.h"
#include "vector.h"
#include "model.h"
#include "rng.h"
#include "randvar.h"
#include "mes.h"

static void lrdecomp(int dim, double **a, double *p);
static void lyequalsb(double **a, double *b, double *p, int dim, double *y);
static void ltranspxequalsy(double **a, double *y, double *p, int dim, 
			    double *x);

/*============================================================================*/

int matrix_d_read(scanner_t *s, double **matrix, int max_zeile, int max_spalte){
#define CUR_PROC "matrix_d_read"
  int len =0, zeile = 0;
  scanner_consume(s, '{'); if (s->err) return(-1);
  while( !s->eof && !s->err && s->c - '}') {
    if (zeile >= max_zeile) {
      scanner_error(s, "too many rows in matrix");
      return(-1);
    } 
    /* Memory allocation takes place in scanner_get_double_earray */
    matrix[zeile] = scanner_get_double_earray(s, &len);
    if (len != max_spalte) {
      scanner_error(s, "wrong number of elements in matrix");
      return(-1);
    }
    scanner_consume(s, ';'); 
    if (s->err) {
      scanner_error(s, "missing ';' or wrong number of columns");
      return(-1);
    }
    zeile++;
  }
  scanner_consume(s, '}'); if (s->err) return(-1);
  if (zeile < max_zeile){
    scanner_error(s, "rows missing in matrix");
    return(-1);
  } 
  return(0);  
#undef CUR_PROC
} /* matrix_d_read */

/*============================================================================*/

int matrix_i_read(scanner_t *s, int **matrix, int max_zeile, int max_spalte) {
#define CUR_PROC "matrix_i_read" 
  int len =0, zeile = 0;
  scanner_consume(s, '{'); if (s->err) return(-1);
  while( !s->eof && !s->err && s->c - '}') {
    if (zeile >= max_zeile) {
      scanner_error(s, "too many rows in matrix");
      return(-1);
    } 
    /* Memory allocation takes place in scanner_get_int_array */
    matrix[zeile] = scanner_get_int_array(s, &len);
    if (len != max_spalte) {
      scanner_error(s, "wrong number of elements in matrix");
      return(-1);
    }
    scanner_consume(s, ';'); 
    if (s->err) {
      scanner_error(s, "missing ';' or wrong number of columns");
      return(-1);
    }
    zeile++;
  }
  scanner_consume(s, '}'); if (s->err) return(-1);
  if (zeile < max_zeile){
    scanner_error(s, "rows missing in matrix");
    return(-1);
  } 
  return(0);  
#undef CUR_PROC
} /* matrix_i_read */

/*============================================================================*/

/* allocation of matrices with fixed dimensions  */
double** stat_matrix_d_alloc(int n, int m) {
#define CUR_PROC "stat_matrix_d_alloc"
  int i, j;
  double **A;
  double *tmp;
  
  if (!m_calloc(A,  n * sizeof(double*) +  n * m * sizeof(double)) ){
	mes_proc(); 
	goto STOP;
  }
  
  tmp = (double*)(A + n);
  for (i = 0; i < n; i++) {
    A[i] = tmp;
    tmp += m;
  }
  return A;
STOP:
  stat_matrix_d_free(&A);
  return NULL;
#undef CUR_PROC
}


int stat_matrix_d_free(double ***matrix) {
#define CUR_PROC "stat_matrix_d_free"
  mes_check_ptr(matrix, return(-1));
  if ( !*matrix) return(0);
  free(*matrix);
  return 0;
#undef CUR_PROC
}


/*============================================================================*/


double** matrix_d_alloc(int zeilen, int spalten) {
#define CUR_PROC "matrix_d_alloc"
  double **matrix;
  int i;
  
  //printf("*** matrix_d_alloc %d zeilen, %d spalten:\n",zeilen, spalten);
  
  if (!m_calloc(matrix, zeilen)) {mes_proc(); goto STOP;}
  for (i = 0; i < zeilen; i++)
    if (!m_calloc(matrix[i], spalten)) {mes_proc(); goto STOP;}
  return matrix;
STOP:
  matrix_d_free(&matrix, zeilen);
  return NULL;
#undef CUR_PROC
} /* matrix_d_alloc */


int matrix_d_free(double ***matrix, long zeilen) {
# define CUR_PROC "matrix_d_free"
  long i;
  mes_check_ptr(matrix, return(-1));
  if ( !*matrix) return(0);
  for (i = zeilen - 1; i >=  0; i--) 
    m_free((*matrix)[i]);
  m_free(*matrix);
  return (0);
# undef CUR_PROC
} /* matrix_d_free */



/*============================================================================*/



double** matrix_d_alloc_copy(int zeilen, int spalten, double **copymatrix) {
#define CUR_PROC "matrix_d_alloc_copy"
  double **matrix;
  int i, j;
  if (!m_calloc(matrix, zeilen)) {mes_proc(); goto STOP;}
  for (i = 0; i < zeilen; i++) {
    if (!m_calloc(matrix[i], spalten)) {mes_proc(); goto STOP;}
    for (j = 0; j < spalten; j++)
      matrix[i][j] = copymatrix[i][j];
  }
  return matrix;
STOP:
  matrix_d_free(&matrix, zeilen);
  return NULL;
#undef CUR_PROC
} /* matrix_d_alloc_copy */

/*============================================================================*/

int** matrix_i_alloc(int zeilen, int spalten) {
#define CUR_PROC "matrix_i_alloc"
  int **matrix;
  int i;
  if (!m_calloc(matrix, zeilen)) {mes_proc(); goto STOP;}
  for (i = 0; i < zeilen; i++)
    if (!m_calloc(matrix[i], spalten)) {mes_proc(); goto STOP;}
  return matrix;
STOP:
  matrix_i_free(&matrix, zeilen);
  return NULL;
#undef CUR_PROC
} /* matrix_i_alloc */

/*============================================================================*/

int matrix_i_free(int ***matrix, long zeilen) {
# define CUR_PROC "matrix_i_free"
  long i;
  mes_check_ptr(matrix, return(-1));
  if ( !*matrix ) return(0);
  for (i = 0; i < zeilen; i++) 
    m_free((*matrix)[i]);
  m_free(*matrix);
  return (0);
# undef CUR_PROC
} /* matrix_i_free */

/*============================================================================*/

void matrix_d_print(FILE *file, double **matrix, int zeilen, int spalten, 
		    char *tab, char *separator, char *ending) {
  int i;
  for (i = 0; i < zeilen; i++)
    vector_d_print(file, matrix[i], spalten, tab, separator, ending);
} /* matrix_d_print */

/*============================================================================*/

void matrix_d_print_prec(FILE *file, double **matrix, int zeilen, int spalten, 
			 int width, int prec, char *tab, char *separator, 
			 char *ending) {
  int i;
  for (i = 0; i < zeilen; i++)
    vector_d_print_prec(file, matrix[i], spalten, width, prec, 
			tab, separator, ending);
} /* matrix_d_print_prec */

/*============================================================================*/

void matrix_i_print(FILE *file, int **matrix, int zeilen, int spalten,
		    char *tab, char *separator, char *ending) {
  int i;
  for (i = 0; i < zeilen; i++)
    vector_i_print(file, matrix[i], spalten, tab, separator, ending);
} /* matrix_i_print */