lib/clinalg.c

   1 /* clinalg - 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 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); }
  16
  17 /************************************************************************/
  18 /**                                                                    **/
  19 /**                  Machine Epsilon Determination                     **/
  20 /**                                                                    **/
  21 /************************************************************************/
  22
  23 double macheps()
  24 {
  25   static int calculated = FALSE;
  26   static double epsilon = 1.0;
  27
  28   if (! calculated)
  29     while (1.0 + epsilon / 2.0 != 1.0) epsilon = epsilon / 2.0;
  30   calculated = TRUE;
  31   return(epsilon);
  32 }
  33
  34 /************************************************************************/
  35 /**                                                                    **/
  36 /**           Lisp Interfaces to Linear Algebra Routines               **/
  37 /**                                                                    **/
  38 /************************************************************************/
  39
  40 chol_decomp_front(mat, n, dpars)
  41      PTR mat, dpars;
  42      int n;
  43 {
  44   double *dp = (double *) dpars;
  45   choldecomp((double **) mat, n, *dp, dp + 1);
  46 }
  47
  48 int lu_decomp_front(mat, n, iv, mode, dp)
  49      PTR mat, iv, dp;
  50      int n, mode;
  51 {
  52   return(crludcmp((char **) mat, n, (int *) iv, mode, (double *) dp));
  53 }
  54
  55 int lu_solve_front(a, n, indx, b, mode)
  56      PTR a, indx, b;
  57      int n, mode;
  58 {
  59   return(crlubksb((char **) a, n, (int *) indx, (char *) b, mode));
  60 }
  61
  62 int lu_inverse_front(pmat, n, piv, pv, mode, pinv)
  63      PTR pmat, piv, pv, pinv;
  64      int n, mode;
  65 {
  66   Matrix mat = (Matrix) pmat, inv = (Matrix) pinv;
  67   IVector iv = (IVector) piv;
  68   Vector v = (Vector) pv;
  69   CMatrix cinv;
  70   RMatrix rinv;
  71   CVector cv;
  72   RVector rv;
  73   double d;
  74   int i, j, singular;
  75
  76   singular = crludcmp(mat, n, iv, mode, &d);
  77
  78   if (! singular) {
  79     rinv = (RMatrix) inv;
  80     cinv = (CMatrix) inv;
  81     rv = (RVector) v;
  82     cv = (CVector) v;
  83
  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];
  88       }
  89
  90       singular = singular || crlubksb(mat, n, iv, v, mode);
  91
  92       for (i = 0; i < n; i++) {
  93         if (mode == RE) rinv[i][j] = rv[i];
  94         else cinv[i][j] = cv[i];
  95       }
  96     }
  97   }
  98   return(singular);
  99 }
 100
 101 sv_decomp_front(mat, m, n, w, v)
 102      PTR mat, w, v;
 103      int m, n;
 104 {
 105   return(svdcmp((char **) mat, m, n, (char *) w, (char **) v));
 106 }
 107
 108 qr_decomp_front(mat, m, n, v, jpvt, pivot)
 109      PTR mat, v, jpvt;
 110      int m, n, pivot;
 111 {
 112   qrdecomp((char **) mat, m, n, (char **) v, (char *) jpvt, pivot);
 113 }
 114
 115 double rcondest_front(mat, n)
 116      PTR mat;
 117      int n;
 118 {
 119   return(rcondest((char **) mat, n));
 120 }
 121
 122 make_rotation_front(n, rot, x, y, use_alpha, alpha)
 123      int n, use_alpha;
 124      PTR rot, x, y;
 125      double alpha;
 126 {
 127   make_rotation(n, (char **) rot, (char *) x, (char *) y, use_alpha, alpha);
 128 }
 129
 130 int eigen_front(a, n, w, z, fv1)
 131      PTR a, w, z, fv1;
 132      int n;
 133 {
 134   int ierr;
 135
 136   eigen(&n, &n, (char *) a, (char *) w, (char *) z, (char *) fv1, &ierr);
 137   return(ierr);
 138 }
 139
 140 fft_front(n, x, work, isign)
 141      int n, isign;
 142      PTR x, work;
 143 {
 144   cfft(n, (char *) x, (char *) work, isign);
 145 }
 146
 147 int base_lowess_front(x, y, n, f, nsteps, delta, ys, rw, res)
 148      PTR x, y, ys, rw, res;
 149      int n, nsteps;
 150      double f, delta;
 151 {
 152   return(lowess((char *) x, (char *) y, n, f, nsteps, delta,
 153                 (char *) ys, (char *) rw, (char *) res));
 154 }
 155
 156 range_to_rseq(n, px, ns, pxs)
 157      int n, ns;
 158      PTR px, pxs;
 159 {
 160   int i;
 161   double xmin, xmax, *x, *xs;
 162
 163   x = (double *) px;
 164   xs = (double *) pxs;
 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];
 168   }
 169   for (i = 0; i < ns; i++)
 170     xs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (ns - 1));
 171 }
 172
 173 int spline_front(n, x, y, ns, xs, ys, work)
 174      PTR x, y, xs, ys, work;
 175      int n, ns;
 176 {
 177   return(fit_spline(n, (char *) x, (char *) y,
 178                     ns, (char *) xs, (char *) ys, (char *) work));
 179 }
 180
 181 kernel_dens_front(x, n, width, xs, ys, ns, code)
 182      PTR x, xs, ys;
 183      int n, ns, code;
 184      double width;
 185 {
 186   int ktype;
 187
 188   if (code == 0) ktype = 'U';
 189   else if (code == 1) ktype = 'T';
 190   else if (code == 2) ktype = 'G';
 191   else ktype = 'B';
 192
 193   return(kernel_smooth((char *) x, nil, n, width, nil, nil,
 194                        (char *) xs, (char *) ys, ns, ktype));
 195 }
 196
 197 int kernel_smooth_front(x, y, n, width, xs, ys, ns, code)
 198      PTR x, y, xs, ys;
 199      int n, ns, code;
 200      double width;
 201 {
 202   int ktype;
 203
 204   if (code == 0) ktype = 'U';
 205   else if (code == 1) ktype = 'T';
 206   else if (code == 2) ktype = 'G';
 207   else ktype = 'B';
 208
 209   return(kernel_smooth((char *) x, (char *) y, n, width, nil, nil,
 210                        (char *) xs, (char *) ys, ns, ktype));
 211 }