1 /* clinalg - C interface to basic linear algebra routines. */
2 /* Copyright (c) 1990, by Luke Tierney */
12 extern double rcondest();
14 extern void choldecomp();
15 extern int crludcmp();
16 extern int crlubksb();
18 extern void qrdecomp();
19 extern void make_rotation();
23 extern int fit_spline();
24 extern int kernel_smooth();
31 return((x
< y
) ? x
: y
);
37 return((x
> y
) ? x
: y
);
40 /************************************************************************/
42 /** Machine Epsilon Determination **/
44 /************************************************************************/
48 static int calculated
= FALSE
;
49 static double epsilon
= 1.0;
52 while (1.0 + epsilon
/ 2.0 != 1.0) epsilon
= epsilon
/ 2.0;
57 /************************************************************************/
59 /** Lisp Interfaces to Linear Algebra Routines **/
61 /************************************************************************/
64 chol_decomp_front(PTR mat
, int n
, PTR dpars
)
66 double *dp
= (double *) dpars
;
67 choldecomp((double **) mat
, n
, *dp
, dp
+ 1);
71 lu_decomp_front(PTR mat
, int n
, PTR iv
, int mode
, PTR dp
)
73 return(crludcmp((char **) mat
, n
, (int *) iv
, mode
, (double *) dp
));
77 lu_solve_front(PTR a
, int n
, PTR indx
, PTR b
, int mode
)
79 return(crlubksb((char **) a
, n
, (int *) indx
, (char *) b
, mode
));
83 lu_inverse_front(PTR pmat
, int 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
;
95 singular
= crludcmp(mat
, n
, iv
, mode
, &d
);
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
];
121 sv_decomp_front(PTR mat
, int m
, int n
, PTR w
, PTR v
)
123 return(svdcmp((char **) mat
, m
, n
, (char *) w
, (char **) v
));
127 qr_decomp_front(PTR mat
, int m
, int n
, PTR v
, PTR jpvt
, int pivot
)
129 qrdecomp((char **) mat
, m
, n
, (char **) v
, (char *) jpvt
, pivot
);
133 rcondest_front(PTR mat
, int n
)
135 return(rcondest((char **) mat
, n
));
139 make_rotation_front(int 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
, int n
, PTR w
, PTR z
, PTR fv1
)
149 eigen(&n
, &n
, (char *) a
, (char *) w
, (char *) z
, (char *) fv1
, &ierr
);
154 fft_front(int n
, PTR x
, PTR work
, int isign
)
156 cfft(n
, (char *) x
, (char *) work
, isign
);
160 base_lowess_front(PTR x
, PTR y
, int n
, double f
,
161 int 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
));
168 range_to_rseq(int n
, PTR px
, int ns
, PTR pxs
)
171 double xmin
, xmax
, *x
, *xs
;
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(int n
, PTR x
, PTR y
, int 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
, int n
, double width
, PTR xs
, PTR ys
, int ns
, int code
)
196 if (code
== 0) ktype
= 'U';
197 else if (code
== 1) ktype
= 'T';
198 else if (code
== 2) ktype
= 'G';
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
, int n
, double width
,
207 PTR xs
, PTR ys
, int ns
, int code
)
211 if (code
== 0) ktype
= 'U';
212 else if (code
== 1) ktype
= 'T';
213 else if (code
== 2) ktype
= 'G';
216 return(kernel_smooth((char *) x
, (char *) y
, n
, width
, nil
, nil
,
217 (char *) xs
, (char *) ys
, ns
, ktype
));