src/c/linalg.c

   1 /* linalg - Lisp interface for basic linear algebra routines           */
   2 /* XLISP-STAT 2.1 Copyright (c) 1990-1995, by Luke Tierney             */
   3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
   4 /* You may give out copies of this software; for conditions see the    */
   5 /* file COPYING included with this distribution.                       */
   6
   7 #include "linalg.h"
   8
   9 #define checktvecsize(x,n) if ((n) > gettvecsize(x)) xlbadtype(x)
  10 #define checknonneg(n) if ((n) < 0) xlbadtype(cvfixnum((FIXTYPE) (n)))
  11 #define checksquarematrix(x) \
  12   if (! matrixp(x) || numrows(x) != numcols(x)) xlbadtype(x);
  13
  14 #define CXDAT(x) ((dcomplex *) gettvecdata(x))
  15 #define REDAT(x) ((double *) gettvecdata(x))
  16 #define INDAT(x) ((int *) gettvecdata(x))
  17
  18 #define IN 0
  19 #define RE 1
  20 #define CX 2
  21
  22
  23 /************************************************************************/
  24 /**                                                                    **/
  25 /**            Argument Checking and Conversion Functions              **/
  26 /**                                                                    **/
  27 /************************************************************************/
  28
  29 int anycomplex P1C(LVAL, x)
  30 {
  31   LVAL data;
  32
  33   data = compounddataseq(x);
  34
  35   switch (ntype(data)) {
  36   case CONS:
  37     for (; consp(data); data = cdr(data))
  38       if (complexp(car(data)))
  39         return TRUE;
  40     return FALSE;
  41   case VECTOR:
  42     {
  43       int i, n;
  44       n = getsize(data);
  45       for (i = 0; i < n; i++)
  46         if (complexp(getelement(data, i)))
  47           return TRUE;
  48       return FALSE;
  49     }
  50   case TVEC:
  51     switch (gettvectype(data)) {
  52     case CD_CXFIXTYPE:
  53     case CD_CXFLOTYPE:
  54     case CD_COMPLEX:
  55     case CD_DCOMPLEX:
  56       return TRUE;
  57     default:
  58       return FALSE;
  59     }
  60   default:
  61     return FALSE;
  62   }
  63 }
  64
  65 LVAL xsanycomplex(V)
  66 {
  67   while (moreargs())
  68     if (anycomplex(xlgetarg()))
  69       return s_true;
  70   return NIL;
  71 }
  72
  73 LOCAL LVAL newmatrix P2C(int, m, int, n)
  74 {
  75   LVAL dim, result;
  76
  77   xlsave1(dim);
  78   checknonneg(m);
  79   checknonneg(n);
  80   dim = integer_list_2(m, n);
  81   result = mkarray(dim, NIL, NIL, s_true);
  82   xlpop();
  83
  84   return result;
  85 }
  86
  87 LOCAL LVAL getlinalgdata P4C(int, off, int, n, LVAL, arg, int, type)
  88 {
  89   LVAL x;
  90
  91   x = darrayp(arg) ? getdarraydata(arg) : arg;
  92   if (! tvecp(x))
  93     xlbadtype(arg);
  94   if (off < 0 || n < 0 || gettvecsize(x) < off + n)
  95     xlerror("incompatible with access indices", x);
  96   switch (type) {
  97   case IN:
  98     if (gettvectype(x) != CD_INT)
  99       xlbadtype(x);
 100     break;
 101   case RE:
 102     switch(gettvectype(x)) {
 103     case CD_FLOTYPE:
 104     case CD_DOUBLE:
 105       break;
 106     default:
 107       xlbadtype(x);
 108     }
 109     break;
 110   case CX:
 111     switch(gettvectype(x)) {
 112     case CD_CXFLOTYPE:
 113     case CD_DCOMPLEX:
 114       break;
 115     default:
 116       xlbadtype(x);
 117     }
 118     break;
 119   }
 120   return x;
 121 }
 122
 123 #define getlinalgivec(off,n,arg) (INDAT(getlinalgdata(off,n,arg,IN)) + (off))
 124 #define getlinalgdvec(off,n,arg) (REDAT(getlinalgdata(off,n,arg,RE)) + (off))
 125 #define getlinalgzvec(off,n,arg) (CXDAT(getlinalgdata(off,n,arg,CX)) + (off))
 126
 127 LOCAL VOID transposeinto P4C(LVAL, x, int, m, int, n, LVAL, y)
 128 {
 129   int i, j, in, jm;
 130
 131   x = compounddataseq(x);
 132   y = compounddataseq(y);
 133   if (! vectorp(x) && ! tvecp(x) && ! stringp(x)) xlbadtype(x);
 134   if (! vectorp(y) && ! tvecp(y) && ! stringp(y)) xlbadtype(y);
 135   checknonneg(n);
 136   checknonneg(m);
 137   checktvecsize(x, n * m);
 138   checktvecsize(y, n * m);
 139
 140   for (i = 0, in = 0; i < m; i++, in += n)
 141     for (j = 0, jm = 0; j < n; j++, jm += m)
 142       settvecelement(y, jm + i, gettvecelement(x, in + j));
 143 }
 144
 145 LVAL xstransposeinto(V)
 146 {
 147   LVAL x, y;
 148   int m, n;
 149
 150   x = xlgetarg();
 151   m = getfixnum(xlgafixnum());
 152   n = getfixnum(xlgafixnum());
 153   y = xlgetarg();
 154   xllastarg();
 155
 156   transposeinto(x, m, n, y);
 157
 158   return y;
 159 }
 160
 161 LOCAL LVAL gen2linalg P5C(LVAL, arg, int, m, int, n, LVAL, type, int , trans)
 162 {
 163   LVAL x, y;
 164   int mn;
 165
 166   x = compounddataseq(arg);
 167   mn = n * m;
 168
 169   xlsave1(y);
 170   y = mktvec(mn, type);
 171   if (trans)
 172     transposeinto(x, m, n, y);
 173   else
 174     xlreplace(y, x, 0, mn, 0, mn);
 175   xlpop();
 176   return y;
 177 }
 178
 179 LOCAL LVAL linalg2genvec P2C(LVAL, x, int, n)
 180 {
 181   LVAL y;
 182
 183   if (! tvecp(x)) xlbadtype(x);
 184   if (n <= 0 || gettvecsize(x) < n) xlfail("bad dimensions");
 185
 186   xlsave1(y);
 187   y = newvector(n);
 188   xlreplace(y, x, 0, n, 0, n);
 189   xlpop();
 190   return y;
 191 }
 192
 193 LOCAL LVAL linalg2genmat P4C(LVAL, arg, int, m, int, n, int, trans)
 194 {
 195   LVAL x, y;
 196   int mn;
 197
 198   x = compounddataseq(arg);
 199   mn = m * n;
 200   if (! tvecp(x)) xlbadtype(arg);
 201   if (n <= 0 || m <= 0 || gettvecsize(x) < mn) xlfail("bad dimensions");
 202
 203   xlsave1(y);
 204   y = newmatrix(m, n);
 205   if (trans)
 206     transposeinto(x, n, m, y);
 207   else
 208     xlreplace(getdarraydata(y), x, 0, mn, 0, mn);
 209   xlpop();
 210   return y;
 211 }
 212
 213 LVAL xsgen2linalg(V)
 214 {
 215   LVAL x, type;
 216   int m, n, trans;
 217
 218   x = xlgetarg();
 219   m = getfixnum(xlgafixnum());
 220   n = getfixnum(xlgafixnum());
 221   type = xlgetarg();
 222   trans = moreargs() ? ! null(xlgetarg()) : FALSE;
 223   xllastarg();
 224
 225   return gen2linalg(x, m, n, type, trans);
 226 }
 227
 228 LVAL xslinalg2gen(V)
 229 {
 230   LVAL x, d;
 231   int trans;
 232
 233   x = xlgetarg();
 234   d = xlgetarg();
 235   trans = moreargs() ? ! null(xlgetarg()) : FALSE;
 236   xllastarg();
 237
 238   if (fixp(d))
 239     return linalg2genvec(x, getfixnum(d));
 240   else if (consp(d) && consp(cdr(d)) && fixp(car(d)) && fixp(car(cdr(d))))
 241     return linalg2genmat(x, getfixnum(car(d)), getfixnum(car(cdr(d))), trans);
 242   else
 243     xlbadtype(d);
 244   return NIL;
 245 }
 246
 247 LOCAL VOID checkldim P2C(int, lda, int, n)
 248 {
 249   if (lda < 0 || n < 0 || n < lda)
 250     xlfail("bad dimensions");
 251 }
 252
 253
 254 /****************************************************************************/
 255 /*                                                                          */
 256 /*                         LU Decomposition Routines                        */
 257 /*                                                                          */
 258 /****************************************************************************/
 259
 260 LVAL xslpdgeco(V)
 261 {
 262   LVAL a, ipvt, z;
 263   double *da, rcond, *dz;
 264   int lda, offa, n, *dipvt;
 265
 266   a = xlgetarg();
 267   offa = getfixnum(xlgafixnum());
 268   lda = getfixnum(xlgafixnum());
 269   n = getfixnum(xlgafixnum());
 270   ipvt = xlgetarg();
 271   z = xlgetarg();
 272   xllastarg();
 273
 274   checkldim(lda, n);
 275   da = getlinalgdvec(offa, lda * n, a);
 276   dipvt = getlinalgivec(0, n, ipvt);
 277   dz = getlinalgdvec(0, n, z);
 278
 279   linpack_dgeco(da, lda, n, dipvt, &rcond, dz);
 280
 281   return cvflonum((FLOTYPE) rcond);
 282 }
 283
 284 LVAL xslpdgedi(V)
 285 {
 286   LVAL a, ipvt, det, work;
 287   double *da, *ddet, *dwork;
 288   int lda, offa, n, *dipvt, job, i, ilda;
 289
 290   a = xlgetarg();
 291   offa = getfixnum(xlgafixnum());
 292   lda = getfixnum(xlgafixnum());
 293   n = getfixnum(xlgafixnum());
 294   ipvt = xlgetarg();
 295   det = xlgetarg();
 296   work = xlgetarg();
 297   job = getfixnum(xlgafixnum());
 298   xllastarg();
 299
 300   checkldim(lda, n);
 301   da = getlinalgdvec(offa, lda * n, a);
 302   dipvt = getlinalgivec(0, n, ipvt);
 303   ddet = (job / 10 != 0) ? getlinalgdvec(0, 2, det) : NULL;
 304   dwork = getlinalgdvec(0, n, work);
 305
 306   if (job % 10 != 0)
 307     for (i = 0, ilda = 0; i < n; i++, ilda += lda)
 308       if (da[ilda + i] == 0.0)
 309         xlfail("matrix is (numerically) singular");
 310
 311   linpack_dgedi(da, lda, n, dipvt, ddet, dwork, job);
 312
 313   return NIL;
 314 }
 315
 316 LVAL xslpdgefa(V)
 317 {
 318   LVAL a, ipvt;
 319   double *da;
 320   int lda, offa, n, *dipvt, info;
 321
 322   a = xlgetarg();
 323   offa = getfixnum(xlgafixnum());
 324   lda = getfixnum(xlgafixnum());
 325   n = getfixnum(xlgafixnum());
 326   ipvt = xlgetarg();
 327   xllastarg();
 328
 329   checkldim(lda, n);
 330   da = getlinalgdvec(offa, lda * n, a);
 331   dipvt = getlinalgivec(0, n, ipvt);
 332
 333   linpack_dgefa(da, lda, n, dipvt, &info);
 334
 335   return cvfixnum((FIXTYPE) info);
 336 }
 337
 338 LVAL xslpdgesl(V)
 339 {
 340   LVAL a, ipvt, b;
 341   double *da, *db;
 342   int lda, offa, n, *dipvt, job, i, ilda;
 343
 344   a = xlgetarg();
 345   offa = getfixnum(xlgafixnum());
 346   lda = getfixnum(xlgafixnum());
 347   n = getfixnum(xlgafixnum());
 348   ipvt = xlgetarg();
 349   b = xlgetarg();
 350   job = getfixnum(xlgafixnum());
 351   xllastarg();
 352
 353   checkldim(lda, n);
 354   da = getlinalgdvec(offa, lda * n, a);
 355   dipvt = getlinalgivec(0, n, ipvt);
 356   db = getlinalgdvec(0, n, b);
 357
 358   for (i = 0, ilda = 0; i < n; i++, ilda += lda)
 359     if (da[ilda + i] == 0.0)
 360       xlfail("matrix is (numerically) singular");
 361
 362   linpack_dgesl(da, lda, n, dipvt, db, job);
 363
 364   return NIL;
 365 }
 366
 367 LVAL xslpzgeco(V)
 368 {
 369   LVAL a, ipvt, z;
 370   dcomplex *da, *dz;
 371   double rcond;
 372   int lda, offa, n, *dipvt;
 373
 374   a = xlgetarg();
 375   offa = getfixnum(xlgafixnum());
 376   lda = getfixnum(xlgafixnum());
 377   n = getfixnum(xlgafixnum());
 378   ipvt = xlgetarg();
 379   z = xlgetarg();
 380   xllastarg();
 381
 382   checkldim(lda, n);
 383   da = getlinalgzvec(offa, lda * n, a);
 384   dipvt = getlinalgivec(0, n, ipvt);
 385   dz = getlinalgzvec(0, n, z);
 386
 387   linpack_zgeco(da, lda, n, dipvt, &rcond, dz);
 388
 389   return cvflonum((FLOTYPE) rcond);
 390 }
 391
 392 LVAL xslpzgedi(V)
 393 {
 394   LVAL a, ipvt, det, work;
 395   dcomplex *da, *ddet, *dwork;
 396   int lda, offa, n, *dipvt, job, i, ilda;
 397
 398   a = xlgetarg();
 399   offa = getfixnum(xlgafixnum());
 400   lda = getfixnum(xlgafixnum());
 401   n = getfixnum(xlgafixnum());
 402   ipvt = xlgetarg();
 403   det = xlgetarg();
 404   work = xlgetarg();
 405   job = getfixnum(xlgafixnum());
 406   xllastarg();
 407
 408   checkldim(lda, n);
 409   da = getlinalgzvec(offa, lda * n, a);
 410   dipvt = getlinalgivec(0, n, ipvt);
 411   ddet = (job / 10 != 0) ? getlinalgzvec(0, 2, det) : NULL;
 412   dwork = getlinalgzvec(0, n, work);
 413
 414   if (job % 10 != 0)
 415     for (i = 0, ilda = 0; i < n; i++, ilda += lda)
 416       if (da[ilda + i].r == 0.0 && da[ilda + i].i == 0.0)
 417         xlfail("matrix is (numerically) singular");
 418
 419   linpack_zgedi(da, lda, n, dipvt, ddet, dwork, job);
 420
 421   return NIL;
 422 }
 423
 424 LVAL xslpzgefa(V)
 425 {
 426   LVAL a, ipvt;
 427   dcomplex *da;
 428   int lda, offa, n, *dipvt, info;
 429
 430   a = xlgetarg();
 431   offa = getfixnum(xlgafixnum());
 432   lda = getfixnum(xlgafixnum());
 433   n = getfixnum(xlgafixnum());
 434   ipvt = xlgetarg();
 435   xllastarg();
 436
 437   checkldim(lda, n);
 438   da = getlinalgzvec(offa, lda * n, a);
 439   dipvt = getlinalgivec(0, n, ipvt);
 440
 441   linpack_zgefa(da, lda, n, dipvt, &info);
 442
 443   return cvfixnum((FIXTYPE) info);
 444 }
 445
 446 LVAL xslpzgesl(V)
 447 {
 448   LVAL a, ipvt, b;
 449   dcomplex *da, *db;
 450   int lda, offa, n, *dipvt, job, i, ilda;
 451
 452   a = xlgetarg();
 453   offa = getfixnum(xlgafixnum());
 454   lda = getfixnum(xlgafixnum());
 455   n = getfixnum(xlgafixnum());
 456   ipvt = xlgetarg();
 457   b = xlgetarg();
 458   job = getfixnum(xlgafixnum());
 459   xllastarg();
 460
 461   checkldim(lda, n);
 462   da = getlinalgzvec(offa, lda * n, a);
 463   dipvt = getlinalgivec(0, n, ipvt);
 464   db = getlinalgzvec(0, n, b);
 465
 466   for (i = 0, ilda = 0; i < n; i++, ilda += lda)
 467     if (da[ilda + i].r == 0.0 && da[ilda + i].i == 0.0)
 468       xlfail("matrix is (numerically) singular");
 469
 470   linpack_zgesl(da, lda, n, dipvt, db, job);
 471
 472   return NIL;
 473 }
 474
 475 /****************************************************************************/
 476 /*                                                                          */
 477 /*                         SV Decomposition Routines                        */
 478 /*                                                                          */
 479 /****************************************************************************/
 480
 481 LVAL xslpdsvdc(V)
 482 {
 483   LVAL x, s, e, u, v, work;
 484   int n, p, job, info, jobu, jobv, ncu, offx, ldx, offu, ldu, offv, ldv;
 485   double *dx, *ds, *de, *du, *dv, *dwork;
 486
 487   x = xlgetarg();
 488   offx = getfixnum(xlgafixnum());
 489   ldx = getfixnum(xlgafixnum());
 490   n = getfixnum(xlgafixnum());
 491   p = getfixnum(xlgafixnum());
 492   s = xlgetarg();
 493   e = xlgetarg();
 494   u = xlgetarg();
 495   offu = getfixnum(xlgafixnum());
 496   ldu = getfixnum(xlgafixnum());
 497   v = xlgetarg();
 498   offv = getfixnum(xlgafixnum());
 499   ldv = getfixnum(xlgafixnum());
 500   work = xlgetarg();
 501   job = getfixnum(xlgafixnum());
 502   xllastarg();
 503
 504   jobu = job % 100 / 10;
 505   ncu = jobu > 1 ? (n < p ? n : p) : n;
 506   jobv = job % 10;
 507
 508   checkldim(ldx, n);
 509   if (jobu) checkldim(ldu, n);
 510   if (jobv) checkldim(ldv, p);
 511
 512   dx = getlinalgdvec(offx, ldx * p, x);
 513   ds = getlinalgdvec(0, n + 1 < p ? n + 1 : p, s);
 514   de = getlinalgdvec(0, p, e);
 515   du = jobu > 0 ? getlinalgdvec(offu, ldu * ncu, u) : NULL;
 516   dv = jobv > 0 ? getlinalgdvec(offv, ldv * p, v) : NULL;
 517   dwork = getlinalgdvec(0, n, work);
 518
 519   linpack_dsvdc(dx, ldx, n, p, ds, de, du, ldu, dv, ldv, dwork, job, &info);
 520   return info ? cvfixnum((FIXTYPE) (info - 1)) : NIL;
 521 }
 522
 523 LVAL xslpzsvdc(V)
 524 {
 525   LVAL x, s, e, u, v, work;
 526   int n, p, job, info, jobu, jobv, ncu, offx, ldx, offu, ldu, offv, ldv;
 527   dcomplex *dx, *ds, *de, *du, *dv, *dwork;
 528
 529   x = xlgetarg();
 530   offx = getfixnum(xlgafixnum());
 531   ldx = getfixnum(xlgafixnum());
 532   n = getfixnum(xlgafixnum());
 533   p = getfixnum(xlgafixnum());
 534   s = xlgetarg();
 535   e = xlgetarg();
 536   u = xlgetarg();
 537   offu = getfixnum(xlgafixnum());
 538   ldu = getfixnum(xlgafixnum());
 539   v = xlgetarg();
 540   offv = getfixnum(xlgafixnum());
 541   ldv = getfixnum(xlgafixnum());
 542   work = xlgetarg();
 543   job = getfixnum(xlgafixnum());
 544   xllastarg();
 545
 546   jobu = job % 100 / 10;
 547   ncu = jobu > 1 ? (n < p ? n : p) : n;
 548   jobv = job % 10;
 549
 550   checkldim(ldx, n);
 551   if (jobu) checkldim(ldu, n);
 552   if (jobv) checkldim(ldv, p);
 553
 554   dx = getlinalgzvec(offx, ldx * p, x);
 555   ds = getlinalgzvec(0, n + 1 < p ? n + 1 : p, s);
 556   de = getlinalgzvec(0, p, e);
 557   du = jobu > 0 ? getlinalgzvec(offu, ldu * ncu, u) : NULL;
 558   dv = jobv > 0 ? getlinalgzvec(offv, ldv * p, v) : NULL;
 559   dwork = getlinalgzvec(0, n, work);
 560
 561   linpack_zsvdc(dx, ldx, n, p, ds, de, du, ldu, dv, ldv, dwork, job, &info);
 562   return info ? cvfixnum((FIXTYPE) (info - 1)) : NIL;
 563 }
 564
 565
 566 /****************************************************************************/
 567 /*                                                                          */
 568 /*                         QR Decomposition Routines                        */
 569 /*                                                                          */
 570 /****************************************************************************/
 571
 572 LVAL xslpdqrdc(V)
 573 {
 574   LVAL x, a, j, w, r, q;
 575   double *dx, *da, *dw, *dr, *dq;
 576   int offx, ldx, n, p, *dj, job;
 577
 578   x = xlgetarg();
 579   offx = getfixnum(xlgafixnum());
 580   ldx = getfixnum(xlgafixnum());
 581   n = getfixnum(xlgafixnum());
 582   p = getfixnum(xlgafixnum());
 583   a = xlgetarg();
 584   j = xlgetarg();
 585   w = xlgetarg();
 586   job = getfixnum(xlgafixnum());
 587   r = moreargs() ? xlgetarg() : NIL;
 588   q = moreargs() ? xlgetarg() : NIL;
 589   xllastarg();
 590
 591   if (p > n && (! null(r) || ! null(q)))
 592     xlfail("more columns than rows");
 593   checkldim(ldx, n);
 594
 595   dx = getlinalgdvec(offx, ldx * p, x);
 596   da = getlinalgdvec(0, p, a);
 597   dj = job != 0 ? getlinalgivec(0, p, j) : NULL;
 598   dw = job != 0 ? getlinalgdvec(0, p, w) : NULL;
 599   dr = null(r) ? NULL : getlinalgdvec(0, p * p, r);
 600   dq = null(q) ? NULL : getlinalgdvec(0, n * p, q);
 601
 602   linpack_dqrdc(dx, ldx, n, p, da, dj, dw, job);
 603
 604   if (! null(r)) {
 605     int i, j, ip, jn;
 606
 607     /* copy the upper triangle of X to R in row major order */
 608     for (i = 0, ip = 0; i < p; i++, ip += p) {
 609       for (j = 0; j < i; j++)
 610         dr[ip + j] = 0.0;
 611       for (j = i, jn = j * n; j < p; j++, jn += n)
 612         dr[ip + j] = dx[jn + i];
 613     }
 614   }
 615
 616   if (! null(q)) {
 617     int i, j, ip, jn, jp;
 618     double t;
 619
 620     /* copy X into Q in row major order */
 621     for (i = 0, ip = 0; i < n; i++, ip += p)
 622       for (j = 0, jn = 0; j < p; j++, jn += n)
 623         dq[ip + j] = dx[jn + i];
 624
 625     /* accumulate the Q transformation */
 626     for (i = 0, ip = 0; i < p; i++, ip += p) {
 627       dq[ip + i] = da[i];
 628       for (j = i + 1; j < p; j++)
 629         dq[ip + j] = 0.0;
 630     }
 631
 632     for (i = p - 1, ip = i * p; i >= 0; i--, ip -= p) {
 633       if (i == n - 1)
 634         dq[ip + i] = 1.0;
 635       else {
 636         for (j = i, jp = ip; j < n; j++, jp += p)
 637           dq[jp + i] = -dq[jp + i];
 638         dq[ip + i] += 1.0;
 639       }
 640       for (j = i - 1, jp = ip - p; j >= 0; j--, jp -= p) {
 641         if (dq[jp + j] != 0.0) {
 642           t = -blas_ddot(n - j, dq + jp + j, p, dq + jp + i, p) / dq[jp + j];
 643           blas_daxpy(n - j, t, dq + jp + j, p, dq + jp + i, p);
 644         }
 645       }
 646     }
 647   }
 648
 649   return NIL;
 650 }
 651
 652 LVAL xslpdqrsl(V)
 653 {
 654   LVAL x, qraux, y, qy, qty, b, rsd, xb;
 655   int n, k, job, info, cqy, cqty, cb, cr, cxb, offx, ldx;
 656   double *dx, *dqraux, *dy, *dqy, *dqty, *db, *drsd, *dxb;
 657
 658   x = xlgetarg();
 659   offx = getfixnum(xlgafixnum());
 660   ldx = getfixnum(xlgafixnum());
 661   n = getfixnum(xlgafixnum());
 662   k = getfixnum(xlgafixnum());
 663   qraux = xlgetarg();
 664   y = xlgetarg();
 665   qy = xlgetarg();
 666   qty = xlgetarg();
 667   b = xlgetarg();
 668   rsd = xlgetarg();
 669   xb = xlgetarg();
 670   job = getfixnum(xlgafixnum());
 671   xllastarg();
 672
 673   cqy = job / 10000 != 0;
 674   cqty = job % 10000 != 0;
 675   cb = job % 1000 / 100 != 0;
 676   cr = job % 100 / 10 != 0;
 677   cxb = job % 10 != 0;
 678
 679   if (k > n)
 680     xlfail("more columns than rows");
 681   checkldim(ldx, n);
 682
 683   dx = getlinalgdvec(offx, ldx * k, x);
 684   dqraux = getlinalgdvec(0, k, qraux);
 685   dy = getlinalgdvec(0, n, y);
 686   dqy = cqy ? getlinalgdvec(0, n, qy) : NULL;
 687   dqty = cqty || cb || cr || cxb ? getlinalgdvec(0, n, qty) : NULL;
 688   db = cb ? getlinalgdvec(0, k, b) : NULL;
 689   drsd = cr ? getlinalgdvec(0, n, rsd) : NULL;
 690   dxb = cxb ? getlinalgdvec(0, n, xb) : NULL;
 691
 692   linpack_dqrsl(dx, ldx, n, k, dqraux,
 693                 dy, dqy, dqty, db, drsd, dxb, job, &info);
 694   return info ? cvfixnum((FIXTYPE) (info - 1)) : NIL;
 695 }
 696
 697 LVAL xslpzqrdc(V)
 698 {
 699   LVAL x, a, j, w, r, q;
 700   dcomplex *dx, *da, *dw, *dr, *dq;
 701   int offx, ldx, n, p, *dj, job;
 702
 703   x = xlgetarg();
 704   offx = getfixnum(xlgafixnum());
 705   ldx = getfixnum(xlgafixnum());
 706   n = getfixnum(xlgafixnum());
 707   p = getfixnum(xlgafixnum());
 708   a = xlgetarg();
 709   j = xlgetarg();
 710   w = xlgetarg();
 711   job = getfixnum(xlgafixnum());
 712   r = moreargs() ? xlgetarg() : NIL;
 713   q = moreargs() ? xlgetarg() : NIL;
 714   xllastarg();
 715
 716   if (p > n && (! null(r) || ! null(q)))
 717     xlfail("more columns than rows");
 718   checkldim(ldx, n);
 719
 720   dx = getlinalgzvec(offx, ldx * p, x);
 721   da = getlinalgzvec(0, p, a);
 722   dj = job != 0 ? getlinalgivec(0, p, j) : NULL;
 723   dw = job != 0 ? getlinalgzvec(0, p, w) : NULL;
 724   dr = null(r) ? NULL : getlinalgzvec(0, p * p, r);
 725   dq = null(q) ? NULL : getlinalgzvec(0, n * p, q);
 726
 727   linpack_zqrdc(dx, ldx, n, p, da, dj, dw, job);
 728
 729   if (! null(r)) {
 730     int i, j, ip, jn;
 731
 732     /* copy the upper triangle of X to R in row major order */
 733     for (i = 0, ip = 0; i < p; i++, ip += p) {
 734       for (j = 0; j < i; j++)
 735         dr[ip + j].r = 0.0,
 736         dr[ip + j].i = 0.0;
 737       for (j = i, jn = j * n; j < p; j++, jn += n)
 738         dr[ip + j] = dx[jn + i];
 739     }
 740   }
 741
 742   if (! null(q)) {
 743     int i, j, ip, jn, jp;
 744     dcomplex t;
 745
 746     /* copy X into Q in row major order */
 747     for (i = 0, ip = 0; i < n; i++, ip += p)
 748       for (j = 0, jn = 0; j < p; j++, jn += n)
 749         dq[ip + j] = dx[jn + i];
 750
 751     /* accumulate the Q transformation */
 752     for (i = 0, ip = 0; i < p; i++, ip += p) {
 753       dq[ip + i] = da[i];
 754       for (j = i + 1; j < p; j++)
 755         dq[ip + j].r = 0.0,
 756         dq[ip + j].i = 0.0;
 757     }
 758
 759     for (i = p - 1, ip = i * p; i >= 0; i--, ip -= p) {
 760       if (i == n - 1)
 761         dq[ip + i].r = 1.0,
 762         dq[ip + i].i = 0.0;
 763       else {
 764         for (j = i, jp = ip; j < n; j++, jp += p)
 765           dq[jp + i].r = -dq[jp + i].r,
 766           dq[jp + i].i = -dq[jp + i].i;
 767         dq[ip + i].r += 1.0;
 768       }
 769       for (j = i - 1, jp = ip - p; j >= 0; j--, jp -= p) {
 770         if (dq[jp + j].r != 0.0 || dq[jp + j].i != 0.0) {
 771           blas_zdotc(&t, n - j, dq + jp + j, p, dq + jp + i, p);
 772           t.r = -t.r,
 773           t.i = -t.i;
 774           z_div(&t, &t, &dq[jp + j]);
 775           blas_zaxpy(n - j, &t, dq + jp + j, p, dq + jp + i, p);
 776         }
 777       }
 778     }
 779   }
 780
 781   return NIL;
 782 }
 783
 784 LVAL xslpzqrsl(V)
 785 {
 786   LVAL x, qraux, y, qy, qty, b, rsd, xb;
 787   int n, k, job, info, cqy, cqty, cb, cr, cxb, offx, ldx;
 788   dcomplex *dx, *dqraux, *dy, *dqy, *dqty, *db, *drsd, *dxb;
 789
 790   x = xlgetarg();
 791   offx = getfixnum(xlgafixnum());
 792   ldx = getfixnum(xlgafixnum());
 793   n = getfixnum(xlgafixnum());
 794   k = getfixnum(xlgafixnum());
 795   qraux = xlgetarg();
 796   y = xlgetarg();
 797   qy = xlgetarg();
 798   qty = xlgetarg();
 799   b = xlgetarg();
 800   rsd = xlgetarg();
 801   xb = xlgetarg();
 802   job = getfixnum(xlgafixnum());
 803   xllastarg();
 804
 805   cqy = job / 10000 != 0;
 806   cqty = job % 10000 != 0;
 807   cb = job % 1000 / 100 != 0;
 808   cr = job % 100 / 10 != 0;
 809   cxb = job % 10 != 0;
 810
 811   if (k > n)
 812     xlfail("more columns than rows");
 813   checkldim(ldx, n);
 814
 815   dx = getlinalgzvec(offx, ldx * k, x);
 816   dqraux = getlinalgzvec(0, k, qraux);
 817   dy = getlinalgzvec(0, n, y);
 818   dqy = cqy ? getlinalgzvec(0, n, qy) : NULL;
 819   dqty = cqty || cb || cr || cxb ? getlinalgzvec(0, n, qty) : NULL;
 820   db = cb ? getlinalgzvec(0, k, b) : NULL;
 821   drsd = cr ? getlinalgzvec(0, n, rsd) : NULL;
 822   dxb = cxb ? getlinalgzvec(0, n, xb) : NULL;
 823
 824   linpack_zqrsl(dx, ldx, n, k, dqraux,
 825                 dy, dqy, dqty, db, drsd, dxb, job, &info);
 826   return info ? cvfixnum((FIXTYPE) (info - 1)) : NIL;
 827 }
 828
 829
 830 /************************************************************************/
 831 /**                                                                    **/
 832 /**                    Cholesky Decomposition                          **/
 833 /**                                                                    **/
 834 /************************************************************************/
 835
 836 LVAL xschol_decomp(V)
 837 {
 838   LVAL a, da, val;
 839   int n;
 840   double maxoffl, maxadd;
 841
 842   a = xlgadarray();
 843   maxoffl = moreargs() ? makefloat(xlgetarg()) : 0.0;
 844   xllastarg();
 845
 846   checksquarematrix(a);
 847   n = numrows(a);
 848
 849   xlstkcheck(2);
 850   xlsave(da);
 851   xlsave(val);
 852
 853   da = gen2linalg(a, n, n, s_c_double, FALSE);
 854   choldecomp(REDAT(da), n, maxoffl, &maxadd);
 855
 856   val = consa(cvflonum((FLOTYPE) maxadd));
 857   val = cons(linalg2genmat(da, n, n, FALSE), val);
 858
 859   xlpopn(2);
 860
 861   return val;
 862 }
 863
 864
 865 /************************************************************************/
 866 /**                                                                    **/
 867 /**                 Rotation Matrix Construction                       **/
 868 /**                                                                    **/
 869 /************************************************************************/
 870
 871 LVAL xsmake_rotation(V)
 872 {
 873   LVAL x, y, dx, dy, val;
 874   double alpha=0.0;
 875   int n, use_alpha = FALSE;
 876
 877   x = xlgetarg();
 878   y = xlgetarg();
 879   if (moreargs()) {
 880     use_alpha = TRUE;
 881     alpha = makefloat(xlgetarg());
 882   }
 883   xllastarg();
 884
 885   xlstkcheck(3);
 886   xlsave(dx);
 887   xlsave(dy);
 888   xlsave(val);
 889
 890   dx = coerce_to_tvec(x, s_c_double);
 891   dy = coerce_to_tvec(y, s_c_double);
 892   n = gettvecsize(dx);
 893
 894   if (gettvecsize(dy) != n)
 895     xlfail("sequences not the same length");
 896
 897   val = mktvec(n * n, s_c_double);
 898   make_rotation(n, REDAT(val), REDAT(dx), REDAT(dy), use_alpha, alpha);
 899   val = linalg2genmat(val, n, n, FALSE);
 900
 901   xlpopn(3);
 902
 903   return val;
 904 }
 905
 906
 907 /************************************************************************/
 908 /**                                                                    **/
 909 /**                        EISPACK Routines                            **/
 910 /**                                                                    **/
 911 /************************************************************************/
 912
 913 LVAL xseispackch(V)
 914 {
 915   int nm, n, matz, ierr;
 916   LVAL ar, ai, w, zr, zi, fv1, fv2, fm1;
 917   double *dar, *dai, *dw, *dzr, *dzi, *dfv1, *dfv2, *dfm1;
 918
 919   nm = getfixnum(xlgafixnum());
 920   n = getfixnum(xlgafixnum());
 921   ar = xlgetarg();
 922   ai = xlgetarg();
 923   w = xlgetarg();
 924   matz = getfixnum(xlgafixnum());
 925   zr = xlgetarg();
 926   zi = xlgetarg();
 927   fv1 = xlgetarg();
 928   fv2 = xlgetarg();
 929   fm1 = xlgetarg();
 930   xllastarg();
 931
 932   checkldim(nm, n);
 933   dar = getlinalgdvec(0, nm * n, ar);
 934   dai = getlinalgdvec(0, nm * n, ai);
 935   dw = getlinalgdvec(0, n, w);
 936   dzr = (matz != 0) ? getlinalgdvec(0, nm * n, zr) : NULL;
 937   dzi = (matz != 0) ? getlinalgdvec(0, nm * n, zi) : NULL;
 938   dfv1 = getlinalgdvec(0, n, fv1);
 939   dfv2 = getlinalgdvec(0, n, fv2);
 940   dfm1 = getlinalgdvec(0, 2 * n, fm1);
 941
 942   eispack_ch(nm, n, dar, dai, dw, matz, dzr, dzi, dfv1, dfv2, dfm1, &ierr);
 943   return (ierr == 0) ? NIL : cvfixnum((FIXTYPE) ierr);
 944 }
 945
 946 LVAL xseispackrs(V)
 947 {
 948   int nm, n, matz, ierr;
 949   LVAL a, w, z, fv1, fv2;
 950   double *da, *dw, *dz, *dfv1, *dfv2;
 951
 952   nm = getfixnum(xlgafixnum());
 953   n = getfixnum(xlgafixnum());
 954   a = xlgetarg();
 955   w = xlgetarg();
 956   matz = getfixnum(xlgafixnum());
 957   z = xlgetarg();
 958   fv1 = xlgetarg();
 959   fv2 = xlgetarg();
 960   xllastarg();
 961
 962   checkldim(nm, n);
 963   da = getlinalgdvec(0, nm * n, a);
 964   dw = getlinalgdvec(0, n, w);
 965   dz = (matz != 0) ? getlinalgdvec(0, nm * n, z) : NULL;
 966   dfv1 = getlinalgdvec(0, n, fv1);
 967   dfv2 = getlinalgdvec(0, n, fv2);
 968
 969
 970   eispack_rs(nm, n, da, dw, matz, dz, dfv1, dfv2, &ierr);
 971   return (ierr == 0) ? NIL : cvfixnum((FIXTYPE) ierr);
 972 }
 973
 974
 975 /************************************************************************/
 976 /**                                                                    **/
 977 /**                             A x + y                                **/
 978 /**                                                                    **/
 979 /************************************************************************/
 980
 981 LVAL xsaxpy(V)
 982 {
 983   LVAL result, next, tx, a, x, y;
 984   int i, j, m, n, start, end, lower;
 985   double val;
 986
 987   a = getdarraydata(xlgamatrix());
 988   x = xlgaseq();
 989   y = xlgaseq();
 990   lower = (moreargs() && xlgetarg() != NIL) ? TRUE : FALSE;
 991
 992   n = seqlen(x);
 993   m = seqlen(y);
 994   if (lower && m != n)
 995     xlfail("dimensions do not match");
 996
 997   xlsave1(result);
 998   result = mklist(m, NIL);
 999   for (i = 0, start = 0, next = result;
1000        i < m;
1001        i++, start += n, next = cdr(next)) {
1002     val = makefloat(getnextelement(&y, i));
1003     end = (lower) ? i +1 : n;
1004     for (j = 0, tx = x; j < end; j++) {
1005       val += makefloat(getnextelement(&tx, j))
1006         * makefloat(gettvecelement(a, start + j));
1007     }
1008     rplaca(next, cvflonum((FLOTYPE) val));
1009   }
1010   xlpop();
1011   return(result);
1012 }
1013
1014
1015 /************************************************************************/
1016 /**                                                                    **/
1017 /**                      Fast Fourier Transform                        **/
1018 /**                                                                    **/
1019 /************************************************************************/
1020
1021 LVAL xsfft(V)
1022 {
1023   LVAL data, result, x, work;
1024   int n, isign;
1025
1026   data = xlgaseq();
1027   isign = (moreargs() && xlgetarg() != NIL) ? -1.0 : 1.0;
1028   xllastarg();
1029
1030   /* check and convert the data */
1031   n = seqlen(data);
1032   if (n <= 0)
1033     xlfail("not enough data");
1034
1035   xlstkcheck(2);
1036   xlsave(x);
1037   xlsave(work);
1038   x = gen2linalg(data, n, 1, s_c_dcomplex, FALSE);
1039   work = mktvec(4 * n + 15, s_c_double);
1040
1041   cfft(n, REDAT(x), REDAT(work), isign);
1042
1043   result = listp(x) ? coerce_to_list(x) : coerce_to_tvec(x, s_true);
1044   xlpopn(2);
1045
1046   return result;
1047 }
1048
1049
1050 /************************************************************************/
1051 /**                                                                    **/
1052 /**               Smoothing and Interpolation Routines                 **/
1053 /**                                                                    **/
1054 /************************************************************************/
1055
1056 #define NS_DEFAULT 30
1057
1058 LVAL xsgetsmdata(V)
1059 {
1060   LVAL s1, s2, arg;
1061   LVAL x, y, xs, ys;
1062   int n, ns, i, supplied, is_reg;
1063   double xmin, xmax, *dx, *dxs;
1064
1065   s1 = xlgaseq();
1066   s2 = xlgetarg();
1067   arg = xlgetarg();
1068   is_reg = ! null(xlgetarg());
1069   xllastarg();
1070
1071   if (is_reg && ! seqp(s2))
1072     xlbadtype(s2);
1073   if (! seqp(arg) && ! fixp(arg))
1074     xlbadtype(arg);
1075
1076   ns = (fixp(arg)) ? getfixnum(arg) : seqlen(arg);
1077   supplied = (seqp(arg) && ns >= 1) ? TRUE : FALSE;
1078   if (ns < 1) ns = NS_DEFAULT;
1079
1080   n = seqlen(s1);
1081   if (n <= 0)
1082     xlfail("sequence too short");
1083   if (is_reg && seqlen(s2) != n)
1084     xlfail("sequences not the same length");
1085
1086   xlstkcheck(4);
1087   xlsave(x);
1088   xlsave(y);
1089   xlsave(xs);
1090   xlsave(ys);
1091
1092   x = gen2linalg(s1, n, 1, s_c_double, FALSE);
1093   y = is_reg ? gen2linalg(s2, n, 1, s_c_double, FALSE) : NIL;
1094   xs = supplied ?
1095     gen2linalg(arg, ns, 1, s_c_double, FALSE) : mktvec(ns, s_c_double);
1096   ys = mktvec(ns, s_c_double);
1097
1098   if (! supplied) {
1099     dx = REDAT(x);
1100     dxs = REDAT(xs);
1101     for (xmax = xmin = dx[0], i = 1; i < n; i++) {
1102       if (dx[i] > xmax) xmax = dx[i];
1103       if (dx[i] < xmin) xmin = dx[i];
1104     }
1105     for (i = 0; i < ns; i++)
1106       dxs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (ns - 1));
1107   }
1108
1109   xlnumresults = 0;
1110   xlresults[xlnumresults++] = cvfixnum((FIXTYPE) n);
1111   xlresults[xlnumresults++] = x;
1112   xlresults[xlnumresults++] = y;
1113   xlresults[xlnumresults++] = cvfixnum((FIXTYPE) ns);
1114   xlresults[xlnumresults++] = xs;
1115   xlresults[xlnumresults++] = ys;
1116   xlpopn(4);
1117   return xlresults[0];
1118 }
1119
1120 LVAL xsbasespline(V)
1121 {
1122   LVAL x, y, xs, ys, work;
1123   double *dx, *dy, *dxs, *dys, *dwork;
1124   int n, ns, error;
1125
1126   n = getfixnum(xlgafixnum());
1127   x = xlgetarg();
1128   y = xlgetarg();
1129   ns = getfixnum(xlgafixnum());
1130   xs = xlgetarg();
1131   ys = xlgetarg();
1132   work = xlgetarg();
1133   xllastarg();
1134
1135   dx = getlinalgdvec(0, n, x);
1136   dy = getlinalgdvec(0, n, y);
1137   dxs = getlinalgdvec(0, ns, xs);
1138   dys = getlinalgdvec(0, ns, ys);
1139   dwork = getlinalgdvec(0, 2 * n, work);
1140
1141   error = fit_spline(n, dx, dy, ns, dxs, dys, dwork);
1142
1143   return error ? s_true : NIL;
1144 }
1145
1146 LVAL xsbasekernelsmooth(V)
1147 {
1148   LVAL x, y, xs, ys, targ;
1149   int n, ns, error, ktype;
1150   double *dx, *dy, *dxs, *dys, width;
1151
1152   n = getfixnum(xlgafixnum());
1153   x = xlgetarg();
1154   y = xlgetarg();
1155   ns = getfixnum(xlgafixnum());
1156   xs = xlgetarg();
1157   ys = xlgetarg();
1158   width = makefloat(xlgetarg());
1159   targ = xlgasymbol();
1160   xllastarg();
1161
1162   dx = getlinalgdvec(0, n, x);
1163   dy = null(y) ? NULL : getlinalgdvec(0, n, y);
1164   dxs = getlinalgdvec(0, ns, xs);
1165   dys = getlinalgdvec(0, ns, ys);
1166
1167   switch (getstring(getpname(targ))[0]) {
1168   case 'U': ktype = 'U'; break;
1169   case 'T': ktype = 'T'; break;
1170   case 'G': ktype = 'G'; break;
1171   default:  ktype = 'B'; break;
1172   }
1173
1174   error = kernel_smooth(dx, dy, n, width, NULL, NULL, dxs, dys, ns, ktype);
1175
1176   return error ? s_true : NIL;
1177 }
1178
1179 LVAL xsbaselowess(V)
1180 {
1181   LVAL x, y, ys, rw, res;
1182   double *dx, *dy, *dys, *drw, *dres;
1183   int n, nsteps, error;
1184   double f, delta;
1185
1186   x = xlgetarg();
1187   y = xlgetarg();
1188   n = getfixnum(xlgafixnum());
1189   f = makefloat(xlgetarg());
1190   nsteps = getfixnum(xlgafixnum());
1191   delta = makefloat(xlgetarg());
1192   ys = xlgetarg();
1193   rw = xlgetarg();
1194   res = xlgetarg();
1195   xllastarg();
1196
1197   dx = getlinalgdvec(0, n, x);
1198   dy = getlinalgdvec(0, n, y);
1199   dys = getlinalgdvec(0, n, ys);
1200   drw = getlinalgdvec(0, n, rw);
1201   dres = getlinalgdvec(0, n, res);
1202
1203   error = lowess(dx, dy, n, f, nsteps, delta, dys, drw, dres);
1204
1205   return error ? s_true : NIL;
1206 }
1207
1208 static LVAL add_contour_point P10C(int, m,
1209                                    int, i,
1210                                    int, j,
1211                                    int,  k,
1212                                    int, l,
1213                                    double *, x,
1214                                    double *, y,
1215                                    double *, z,
1216                                    double, v,
1217                                    LVAL, result)
1218 {
1219   LVAL pt;
1220   double p, q;
1221   double zij = z[i * m + j];
1222   double zkl = z[k * m + l];
1223
1224   if ((zij <= v && v < zkl) || (zkl <= v && v < zij)) {
1225     xlsave(pt);
1226     pt = mklist(2, NIL);
1227     p = (v - zij) / (zkl - zij);
1228     q = 1.0 - p;
1229     rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
1230     rplaca(cdr(pt), cvflonum((FLOTYPE) (q * y[j] + p * y[l])));
1231     result = cons(pt, result);
1232     xlpop();
1233   }
1234   return(result);
1235 }
1236
1237 LVAL xssurface_contour(V)
1238 {
1239   LVAL s1, s2, mat, result;
1240   LVAL x, y, z;
1241   double *dx, *dy, *dz;
1242   double v;
1243   int i, j, n, m;
1244
1245   s1 = xlgaseq();
1246   s2 = xlgaseq();
1247   mat = xlgamatrix();
1248   v = makefloat(xlgetarg());
1249   xllastarg();
1250
1251   n = seqlen(s1);
1252   m = seqlen(s2);
1253   if (n != numrows(mat) || m != numcols(mat))
1254     xlfail("dimensions do not match");
1255
1256   xlstkcheck(4);
1257   xlsave(x);
1258   xlsave(y);
1259   xlsave(z);
1260   xlsave(result);
1261
1262   x = gen2linalg(s1,  n, 1, s_c_double, FALSE); dx = REDAT(x);
1263   y = gen2linalg(s2,  m, 1, s_c_double, FALSE); dy = REDAT(y);
1264   z = gen2linalg(mat, n, m, s_c_double, FALSE); dz = REDAT(z);
1265   result = NIL;
1266
1267   for (i = 0; i < n - 1; i++) {
1268     for (j = 0; j < m - 1; j++) {
1269       result = add_contour_point(m, i,   j,   i,   j+1, dx, dy, dz, v, result);
1270       result = add_contour_point(m, i,   j+1, i+1, j+1, dx, dy, dz, v, result);
1271       result = add_contour_point(m, i+1, j+1, i+1, j,   dx, dy, dz, v, result);
1272       result = add_contour_point(m, i+1, j,   i,   j,   dx, dy, dz, v, result);
1273     }
1274   }
1275   xlpopn(4);
1276
1277   return(result);
1278 }
1279
1280 /************************************************************************/
1281 /**                                                                    **/
1282 /**                  Machine Epsilon Determination                     **/
1283 /**                                                                    **/
1284 /************************************************************************/
1285
1286 double macheps(V)
1287 {
1288   static int calculated = FALSE;
1289   static double epsilon = 1.0;
1290
1291   if (! calculated)
1292     while (1.0 + epsilon / 2.0 != 1.0) epsilon = epsilon / 2.0;
1293   calculated = TRUE;
1294   return(epsilon);
1295 }