lib/clinalg.c

   1 /* - C interface to basic linear algebra routines.             */
   2 /* Copyright (c) 1990, by Luke Tierney                                 */
   3
   4 #include "linalg.h"
   5
   6 #ifdef INTPTR
   7 typedef int PTR;
   8 #else
   9 typedef char *PTR;
  10 #endif
  11
  12 extern double rcondest();
  13
  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();
  25
  26
  27
  28 long
  29 min (long x, long y)
  30 {
  31   return((x < y) ? x : y);
  32 }
  33
  34 long
  35 max (long x, long y)
  36 {
  37   return((x > y) ? x : y);
  38 }
  39
  40 /************************************************************************/
  41 /**                                                                    **/
  42 /**                  Machine Epsilon Determination                     **/
  43 /**                                                                    **/
  44 /************************************************************************/
  45
  46 double macheps()
  47 {
  48   static int calculated = FALSE;
  49   static double epsilon = 1.0;
  50
  51   if (! calculated)
  52     while (1.0 + epsilon / 2.0 != 1.0) epsilon = epsilon / 2.0;
  53   calculated = TRUE;
  54   return(epsilon);
  55 }
  56
  57 /************************************************************************/
  58 /**                                                                    **/
  59 /**           Lisp Interfaces to Linear Algebra Routines               **/
  60 /**                                                                    **/
  61 /************************************************************************/
  62
  63 void
  64 chol_decomp_front(PTR mat, size_t n, PTR dpars)
  65 {
  66   double *dp = (double *) dpars;
  67   choldecomp((double **) mat, n, *dp, dp + 1);
  68 }
  69
  70 int
  71 lu_decomp_front(PTR mat, size_t n, PTR iv, int mode, PTR dp)
  72 {
  73   return(crludcmp((char **) mat, n, (int *) iv, mode, (double *) dp));
  74 }
  75
  76 int
  77 lu_solve_front(PTR a, size_t n, PTR indx, PTR b, int mode)
  78 {
  79   return(crlubksb((char **) a, n, (int *) indx, (char *) b, mode));
  80 }
  81
  82 int
  83 lu_inverse_front(PTR pmat, size_t n, PTR piv, PTR pv, int mode, PTR pinv)
  84 {
  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;
  94
  95   singular = crludcmp(mat, n, iv, mode, &d);
  96
  97   if (! singular) {
  98     rinv = (RMatrix) inv;
  99     cinv = (CMatrix) inv;
 100     rv = (RVector) v;
 101     cv = (CVector) v;
 102
 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];
 107       }
 108
 109       singular = singular || crlubksb(mat, n, iv, v, mode);
 110
 111       for (i = 0; i < n; i++) {
 112         if (mode == RE) rinv[i][j] = rv[i];
 113         else cinv[i][j] = cv[i];
 114       }
 115     }
 116   }
 117   return(singular);
 118 }
 119
 120 int
 121 sv_decomp_front(PTR mat, size_t m, size_t n, PTR w, PTR v)
 122 {
 123   return(svdcmp((char **) mat, m, n, (char *) w, (char **) v));
 124 }
 125
 126 void
 127 qr_decomp_front(PTR mat, size_t m, size_t n, PTR v, PTR jpvt, int pivot)
 128 {
 129   qrdecomp((char **) mat, m, n, (char **) v, (char *) jpvt, pivot);
 130 }
 131
 132 double
 133 rcondest_front(PTR mat, size_t n)
 134 {
 135   return(rcondest((char **) mat, n));
 136 }
 137
 138 void
 139 make_rotation_front(size_t n, PTR rot, PTR x, PTR y, int use_alpha, double alpha)
 140 {
 141   make_rotation(n, (char **) rot, (char *) x, (char *) y, use_alpha, alpha);
 142 }
 143
 144 int
 145 eigen_front(PTR a, size_t n, PTR w, PTR z, PTR fv1)
 146 {
 147   int ierr;
 148
 149   eigen(&n, &n, (char *) a, (char *) w, (char *) z, (char *) fv1, &ierr);
 150   return(ierr);
 151 }
 152
 153 void
 154 fft_front(size_t n, PTR x, PTR work, int isign)
 155 {
 156   cfft(n, (char *) x, (char *) work, isign);
 157 }
 158
 159 int
 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)
 162 {
 163   return(lowess((char *) x, (char *) y, n, f, nsteps, delta,
 164                 (char *) ys, (char *) rw, (char *) res));
 165 }
 166
 167 void
 168 range_to_rseq(size_t n, PTR px, long ns, PTR pxs)
 169 {
 170   size_t i;
 171   double xmin, xmax, *x, *xs;
 172
 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];
 178   }
 179   for (i = 0; i < ns; i++)
 180     xs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (ns - 1));
 181 }
 182
 183 int
 184 spline_front(size_t n, PTR x, PTR y, size_t ns,
 185              PTR xs, PTR ys, PTR work)
 186 {
 187   return(fit_spline(n, (char *) x, (char *) y,
 188                     ns, (char *) xs, (char *) ys, (char *) work));
 189 }
 190
 191 int
 192 kernel_dens_front(PTR x, size_t n, double width, PTR xs, PTR ys, size_t ns, int code)
 193 {
 194   int ktype;
 195
 196   if (code == 0) ktype = 'U';
 197   else if (code == 1) ktype = 'T';
 198   else if (code == 2) ktype = 'G';
 199   else ktype = 'B';
 200
 201   return(kernel_smooth((char *) x, nil, n, width, nil, nil,
 202                        (char *) xs, (char *) ys, ns, ktype));
 203 }
 204
 205 int
 206 kernel_smooth_front(PTR x, PTR y, size_t n, double width,
 207                     PTR xs, PTR ys, size_t ns, int code)
 208 {
 209   int ktype;
 210
 211   if (code == 0) ktype = 'U';
 212   else if (code == 1) ktype = 'T';
 213   else if (code == 2) ktype = 'G';
 214   else ktype = 'B';
 215
 216   return(kernel_smooth((char *) x, (char *) y, n, width, nil, nil,
 217                        (char *) xs, (char *) ys, ns, ktype));
 218 }