1 /* clinalg - C interface to basic linear algebra routines. */
2 /* Copyright (c) 1990, by Luke Tierney */
12 extern double rcondest();
14 int min (x
, y
) int x
, y
; { return((x
< y
) ? x
: y
); }
15 int max (x
, y
) int x
, y
; { return((x
> y
) ? x
: y
); }
17 /************************************************************************/
19 /** Machine Epsilon Determination **/
21 /************************************************************************/
25 static int calculated
= FALSE
;
26 static double epsilon
= 1.0;
29 while (1.0 + epsilon
/ 2.0 != 1.0) epsilon
= epsilon
/ 2.0;
34 /************************************************************************/
36 /** Lisp Interfaces to Linear Algebra Routines **/
38 /************************************************************************/
40 chol_decomp_front(mat
, n
, dpars
)
44 double *dp
= (double *) dpars
;
45 choldecomp((double **) mat
, n
, *dp
, dp
+ 1);
48 int lu_decomp_front(mat
, n
, iv
, mode
, dp
)
52 return(crludcmp((char **) mat
, n
, (int *) iv
, mode
, (double *) dp
));
55 int lu_solve_front(a
, n
, indx
, b
, mode
)
59 return(crlubksb((char **) a
, n
, (int *) indx
, (char *) b
, mode
));
62 int lu_inverse_front(pmat
, n
, piv
, pv
, mode
, pinv
)
63 PTR pmat
, piv
, pv
, pinv
;
66 Matrix mat
= (Matrix
) pmat
, inv
= (Matrix
) pinv
;
67 IVector iv
= (IVector
) piv
;
68 Vector v
= (Vector
) pv
;
76 singular
= crludcmp(mat
, n
, iv
, mode
, &d
);
84 for (j
= 0; j
< n
; j
++) {
85 for (i
= 0; i
< n
; i
++) {
86 if (mode
== RE
) rv
[i
] = rinv
[i
][j
];
87 else cv
[i
] = cinv
[i
][j
];
90 singular
= singular
|| crlubksb(mat
, n
, iv
, v
, mode
);
92 for (i
= 0; i
< n
; i
++) {
93 if (mode
== RE
) rinv
[i
][j
] = rv
[i
];
94 else cinv
[i
][j
] = cv
[i
];
101 sv_decomp_front(mat
, m
, n
, w
, v
)
105 return(svdcmp((char **) mat
, m
, n
, (char *) w
, (char **) v
));
108 qr_decomp_front(mat
, m
, n
, v
, jpvt
, pivot
)
112 qrdecomp((char **) mat
, m
, n
, (char **) v
, (char *) jpvt
, pivot
);
115 double rcondest_front(mat
, n
)
119 return(rcondest((char **) mat
, n
));
122 make_rotation_front(n
, rot
, x
, y
, use_alpha
, alpha
)
127 make_rotation(n
, (char **) rot
, (char *) x
, (char *) y
, use_alpha
, alpha
);
130 int eigen_front(a
, n
, w
, z
, fv1
)
136 eigen(&n
, &n
, (char *) a
, (char *) w
, (char *) z
, (char *) fv1
, &ierr
);
140 fft_front(n
, x
, work
, isign
)
144 cfft(n
, (char *) x
, (char *) work
, isign
);
147 int base_lowess_front(x
, y
, n
, f
, nsteps
, delta
, ys
, rw
, res
)
148 PTR x
, y
, ys
, rw
, res
;
152 return(lowess((char *) x
, (char *) y
, n
, f
, nsteps
, delta
,
153 (char *) ys
, (char *) rw
, (char *) res
));
156 range_to_rseq(n
, px
, ns
, pxs
)
161 double xmin
, xmax
, *x
, *xs
;
165 for (xmax
= xmin
= x
[0], i
= 1; i
< n
; i
++) {
166 if (x
[i
] > xmax
) xmax
= x
[i
];
167 if (x
[i
] < xmin
) xmin
= x
[i
];
169 for (i
= 0; i
< ns
; i
++)
170 xs
[i
] = xmin
+ (xmax
- xmin
) * ((double) i
) / ((double) (ns
- 1));
173 int spline_front(n
, x
, y
, ns
, xs
, ys
, work
)
174 PTR x
, y
, xs
, ys
, work
;
177 return(fit_spline(n
, (char *) x
, (char *) y
,
178 ns
, (char *) xs
, (char *) ys
, (char *) work
));
181 kernel_dens_front(x
, n
, width
, xs
, ys
, ns
, code
)
188 if (code
== 0) ktype
= 'U';
189 else if (code
== 1) ktype
= 'T';
190 else if (code
== 2) ktype
= 'G';
193 return(kernel_smooth((char *) x
, nil
, n
, width
, nil
, nil
,
194 (char *) xs
, (char *) ys
, ns
, ktype
));
197 int kernel_smooth_front(x
, y
, n
, width
, xs
, ys
, ns
, code
)
204 if (code
== 0) ktype
= 'U';
205 else if (code
== 1) ktype
= 'T';
206 else if (code
== 2) ktype
= 'G';
209 return(kernel_smooth((char *) x
, (char *) y
, n
, width
, nil
, nil
,
210 (char *) xs
, (char *) ys
, ns
, ktype
));