File: m_jacob.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 (118 lines) | stat: -rw-r--r-- 2,443 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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
/*
 *   Copyright (c) 1996-2001 Lucent Technologies.
 *   See README file for details.
 */

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

#define DEF_METH JAC_EIGD

int jac_reqd(int p) { return(2*p*(p+1)); }

double *jac_alloc(J,p,wk)
jacobian *J;
int p;
double *wk;
{ if (wk==NULL)
    wk = (double *)calloc(2*p*(p+1),sizeof(double));
  J->Z = wk; wk += p*p;
  J->Q = wk; wk += p*p;
  J->wk= wk; wk += p;
  J->dg= wk; wk += p;
  return(wk);
}

void jacob_dec(J, meth)
jacobian *J;
int meth;
{ int i, j, p;

  if (J->st != JAC_RAW) return;

  J->sm = J->st = meth;
  switch(meth)
  { case JAC_EIG:
      eig_dec(J->Z,J->Q,J->p);
      return;
    case JAC_EIGD:
      p = J->p;
      for (i=0; i<p; i++)
        J->dg[i] = (J->Z[i*(p+1)]<=0) ? 0.0 : 1/sqrt(J->Z[i*(p+1)]);
      for (i=0; i<p; i++)
        for (j=0; j<p; j++)
          J->Z[i*p+j] *= J->dg[i]*J->dg[j];
      eig_dec(J->Z,J->Q,J->p);
      J->st = JAC_EIGD;
      return;
    case JAC_CHOL:
      chol_dec(J->Z,J->p);
      return;
    default: printf("jacob_dec: unknown method %d",meth);
  }
}

int jacob_solve(J,v) /* (X^T W X)^{-1} v */
jacobian *J;
double *v;
{ int i, rank;

  if (J->st == JAC_RAW) jacob_dec(J,DEF_METH);

  switch(J->st)
  { case JAC_EIG:
      return(eig_solve(J,v));
    case JAC_EIGD:
      for (i=0; i<J->p; i++) v[i] *= J->dg[i];
      rank = eig_solve(J,v);
      for (i=0; i<J->p; i++) v[i] *= J->dg[i];
      return(rank);
    case JAC_CHOL:
      return(chol_solve(J->Z,v,J->p));
  }
  printf("jacob_solve: unknown method %d",J->st);
  return(0);
}

int jacob_hsolve(J,v) /*  J^{-1/2} v */
jacobian *J;
double *v;
{ int i;

  if (J->st == JAC_RAW) jacob_dec(J,DEF_METH);

  switch(J->st)
  { case JAC_EIG:
      return(eig_hsolve(J,v));
    case JAC_EIGD: /* eigenvalues on corr matrix */
      for (i=0; i<J->p; i++) v[i] *= J->dg[i];
      return(eig_hsolve(J,v));
    case JAC_CHOL:
      return(chol_hsolve(J->Z,v,J->p));
  }
  printf("jacob_hsolve: unknown method %d",J->st);
  return(0);
}

double jacob_qf(J,v)  /* vT J^{-1} v */
jacobian *J;
double *v;
{ int i;

  if (J->st == JAC_RAW) jacob_dec(J,DEF_METH);

  switch (J->st)
  { case JAC_EIG:
      return(eig_qf(J,v));
    case JAC_EIGD:
      for (i=0; i<J->p; i++) v[i] *= J->dg[i];
      return(eig_qf(J,v));
    case JAC_CHOL:
      return(chol_qf(J->Z,v,J->p));
    default:
      printf("jacob_qf: invalid method\n");
      return(0.0);
  }
}