File: m_chol.c

package info (click to toggle)
cufflinks 1.3.0-2
  • links: PTS, VCS
  • area: non-free
  • in suites: wheezy
  • size: 3,864 kB
  • sloc: cpp: 48,999; ansic: 12,297; sh: 3,381; python: 432; makefile: 209
file content (72 lines) | stat: -rw-r--r-- 1,233 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
/*
 *   Copyright (c) 1996-2001 Lucent Technologies.
 *   See README file for details.
 */

#include <math.h>
#include "mutil.h"

void chol_dec(A,n)
double *A;
int n;
{ int i, j, k;
  for (j=0; j<n; j++)
  { k = n*j+j;
    for (i=0; i<j; i++) A[k] -= A[n*i+j]*A[n*i+j];
    if (A[k]<=0)
    { for (i=j; i<n; i++) A[n*j+i] = 0.0; }
    else
    { A[k] = sqrt(A[k]);
      for (i=j+1; i<n; i++)
      { for (k=0; k<j; k++)
          A[n*j+i] -= A[n*k+i]*A[n*k+j];
        A[n*j+i] /= A[n*j+j];
      }
    }
  }
  for (j=0; j<n; j++)
    for (i=j+1; i<n; i++) A[n*i+j] = 0.0;
}

int chol_solve(A,v,p)
double *A, *v;
int p;
{ int i, j;

  for (i=0; i<p; i++)
  { for (j=0; j<i; j++) v[i] -= A[j*p+i]*v[j];
    v[i] /= A[i*p+i];
  }
  for (i=p-1; i>=0; i--)
  { for (j=i+1; j<p; j++) v[i] -= A[i*p+j]*v[j];
    v[i] /= A[i*p+i];
  }
  return(p);
}

int chol_hsolve(A,v,p)
double *A, *v;
int p;
{ int i, j;

  for (i=0; i<p; i++)
  { for (j=0; j<i; j++) v[i] -= A[j*p+i]*v[j];
    v[i] /= A[i*p+i];
  }
  return(p);
}

double chol_qf(A,v,p)
double *A, *v;
int p;
{ int i, j;
  double sum;
 
  sum = 0.0;
  for (i=0; i<p; i++)
  { for (j=0; j<i; j++) v[i] -= A[j*p+i]*v[j];
    v[i] /= A[i*p+i];
    sum += v[i]*v[i];
  }
  return(sum);
}