blob e097ffa56f7f79c887c3590d3eb5aa1a71861213
1 /* - C interface to basic linear algebra routines. */
2 /* Copyright (c) 1990, by Luke Tierney */
4 #include "linalg.h"
6 #ifdef INTPTR
7 typedef int PTR;
8 #else
9 typedef char *PTR;
10 #endif
12 extern double rcondest();
14 extern void choldecomp();
15 extern int crludcmp();
16 extern int crlubksb();
17 extern int svdcmp();
18 extern void qrdecomp();
19 extern void make_rotation();
20 extern void eigen();
21 extern void cfft();
22 extern int lowess();
23 extern int fit_spline();
24 extern int kernel_smooth();
28 long
29 min (long x, long y)
31 return((x < y) ? x : y);
34 long
35 max (long x, long y)
37 return((x > y) ? x : y);
40 /************************************************************************/
41 /** **/
42 /** Machine Epsilon Determination **/
43 /** **/
44 /************************************************************************/
46 double macheps()
48 static int calculated = FALSE;
49 static double epsilon = 1.0;
51 if (! calculated)
52 while (1.0 + epsilon / 2.0 != 1.0) epsilon = epsilon / 2.0;
53 calculated = TRUE;
54 return(epsilon);
57 /************************************************************************/
58 /** **/
59 /** Lisp Interfaces to Linear Algebra Routines **/
60 /** **/
61 /************************************************************************/
63 void
64 chol_decomp_front(PTR mat, size_t n, PTR dpars)
66 double *dp = (double *) dpars;
67 choldecomp((double **) mat, n, *dp, dp + 1);
70 int
71 lu_decomp_front(PTR mat, size_t n, PTR iv, int mode, PTR dp)
73 return(crludcmp((char **) mat, n, (int *) iv, mode, (double *) dp));
76 int
77 lu_solve_front(PTR a, size_t n, PTR indx, PTR b, int mode)
79 return(crlubksb((char **) a, n, (int *) indx, (char *) b, mode));
82 int
83 lu_inverse_front(PTR pmat, size_t n, PTR piv, PTR pv, int mode, PTR pinv)
85 Matrix mat = (Matrix) pmat, inv = (Matrix) pinv;
86 IVector iv = (IVector) piv;
87 Vector v = (Vector) pv;
88 CMatrix cinv;
89 RMatrix rinv;
90 CVector cv;
91 RVector rv;
92 double d;
93 int i, j, singular;
95 singular = crludcmp(mat, n, iv, mode, &d);
97 if (! singular) {
98 rinv = (RMatrix) inv;
99 cinv = (CMatrix) inv;
100 rv = (RVector) v;
101 cv = (CVector) v;
103 for (j = 0; j < n; j++) {
104 for (i = 0; i < n; i++) {
105 if (mode == RE) rv[i] = rinv[i][j];
106 else cv[i] = cinv[i][j];
109 singular = singular || crlubksb(mat, n, iv, v, mode);
111 for (i = 0; i < n; i++) {
112 if (mode == RE) rinv[i][j] = rv[i];
113 else cinv[i][j] = cv[i];
117 return(singular);
121 sv_decomp_front(PTR mat, size_t m, size_t n, PTR w, PTR v)
123 return(svdcmp((char **) mat, m, n, (char *) w, (char **) v));
126 void
127 qr_decomp_front(PTR mat, size_t m, size_t n, PTR v, PTR jpvt, int pivot)
129 qrdecomp((char **) mat, m, n, (char **) v, (char *) jpvt, pivot);
132 double
133 rcondest_front(PTR mat, size_t n)
135 return(rcondest((char **) mat, n));
138 void
139 make_rotation_front(size_t n, PTR rot, PTR x, PTR y, int use_alpha, double alpha)
141 make_rotation(n, (char **) rot, (char *) x, (char *) y, use_alpha, alpha);
145 eigen_front(PTR a, size_t n, PTR w, PTR z, PTR fv1)
147 int ierr;
149 eigen(&n, &n, (char *) a, (char *) w, (char *) z, (char *) fv1, &ierr);
150 return(ierr);
153 void
154 fft_front(size_t n, PTR x, PTR work, int isign)
156 cfft(n, (char *) x, (char *) work, isign);
160 base_lowess_front(PTR x, PTR y, size_t n, double f,
161 size_t nsteps, double delta, PTR ys, PTR rw, PTR res)
163 return(lowess((char *) x, (char *) y, n, f, nsteps, delta,
164 (char *) ys, (char *) rw, (char *) res));
167 void
168 range_to_rseq(size_t n, PTR px, long ns, PTR pxs)
170 size_t i;
171 double xmin, xmax, *x, *xs;
173 x = (double *) px;
174 xs = (double *) pxs;
175 for (xmax = xmin = x[0], i = 1; i < n; i++) {
176 if (x[i] > xmax) xmax = x[i];
177 if (x[i] < xmin) xmin = x[i];
179 for (i = 0; i < ns; i++)
180 xs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (ns - 1));
184 spline_front(size_t n, PTR x, PTR y, size_t ns,
185 PTR xs, PTR ys, PTR work)
187 return(fit_spline(n, (char *) x, (char *) y,
188 ns, (char *) xs, (char *) ys, (char *) work));
192 kernel_dens_front(PTR x, size_t n, double width, PTR xs, PTR ys, size_t ns, int code)
194 int ktype;
196 if (code == 0) ktype = 'U';
197 else if (code == 1) ktype = 'T';
198 else if (code == 2) ktype = 'G';
199 else ktype = 'B';
201 return(kernel_smooth((char *) x, nil, n, width, nil, nil,
202 (char *) xs, (char *) ys, ns, ktype));
206 kernel_smooth_front(PTR x, PTR y, size_t n, double width,
207 PTR xs, PTR ys, size_t ns, int code)
209 int ktype;
211 if (code == 0) ktype = 'U';
212 else if (code == 1) ktype = 'T';
213 else if (code == 2) ktype = 'G';
214 else ktype = 'B';
216 return(kernel_smooth((char *) x, (char *) y, n, width, nil, nil,
217 (char *) xs, (char *) ys, ns, ktype));