invspd.c

这是用matlab编写的支持向量机的函数工具箱

/* invspd.c: Do an "in place" inversion of a real square symmetric positive
 * definite matrix "a" of size "n" by "n" and return log of it's determinant.
 *
 * The function only looks at elements on or above the diagonal. Computes and
 * stores the lower triangular Cholesky factor in "a" (1/6n^3 flops) and does
 * inversion using forward (1/2n^3 flops) and backward (1/6n^3 flops)
 * substitution. On return, the upper diagonal matrix contains the inverse of
 * "a". See: Golub and Van Loan, "Matrix computations", 2nd edition, Johns
 * Hopkins University Press, 1989.
 *
 * (c) Copyright 1996 Carl Edward Rasmussen */

#include <math.h>
#include <stdio.h>
#define real double

real invspd(real **a, int n) 
{
  int  i, j ,k;
  real s, *d, *x;

  d = (real *) malloc((size_t) n*sizeof(real)); 
  x = (real *) malloc((size_t) n*sizeof(real)); 

  for (i=0; i<n; i++)                          /* do Cholesky decomposition */
    for (j=i; j<n; j++) {
      s = a[i][j];
      for (k=i-1; k>=0; k--) s -= a[i][k]*a[j][k];
      if (i == j) {
        if (s <= 0.0) {
          fprintf(stderr,
             "Error: Matrix for inversion is not positive definite... bye!\n");
          exit(-1);
	}
        d[i] = sqrt(s);
      } else a[j][i] = s/d[i];
    }

  for (i=0; i<n; i++) {                                    /* for each colum */
    for (j=0; j<n; j++) {                         /* do forward substitution */
      s = (i == j) ? 1.0 : 0.0;                            /* of unit matrix */
      for (k=j-1; k>=0; k--) s -= a[j][k]*x[k];
      x[j] = s/d[j];
    }
    for (j=n-1; j>=i; j--) {                     /* do backward substitution */
      s = x[j]; 
      for (k=j+1; k<n; k++) s -= a[k][j]*x[k];
      a[i][j] = x[j] = s/d[j];
    }
  } 

  s = 0.0; for (i=0; i<n; i++) s += log(d[i]);          /* compute log det a */
  free(x); free(d);
  return 2.0*s;
}