home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / src / src2.lzh / XLisp-Stat / qrdecomp.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  4KB  |  207 lines
  1. /* adapted from DQRDC of LINPACK */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, 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.  
  8. #include "xlisp.h"
  9. #include "osdef.h"
  10. #ifdef ANSI
  11. #include "xlsproto.h"
  12. #else
  13. #include "xlsfun.h"
  14. #endif ANSI
  15.  
  16. #ifdef ANSI
  17. double NORM2(int,int,int,double **),DOT(int,int,int,int,double **);
  18. void AXPY(int,int,int,int,double,double **),
  19.      SCALE(int,int,int,double,double **),SWAP(int,int,int,double **);
  20. #else
  21. double NORM2(),DOT();
  22. void AXPY(),SCALE(),SWAP();
  23. #endif ANSI
  24.  
  25. #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
  26.  
  27. static double NORM2(i, j, n, x)
  28.      int i, j, n;
  29.      double **x;
  30. {
  31.   int k;
  32.   double maxx, sum, temp;
  33.  
  34.   for (k = i, maxx = 0.0; k < n; k++) {
  35.     temp = fabs(x[k][j]);
  36.     if (maxx < temp) maxx = temp;
  37.   }
  38.   if (maxx == 0.0) return(0.0);
  39.   else {
  40.     for (k = i, sum = 0.0; k < n; k++) {
  41.       temp = x[k][j] / maxx;
  42.       sum += temp * temp;
  43.     }
  44.     return(maxx * sqrt(sum));
  45.   }
  46. }
  47.  
  48. static double DOT(i, j, k, n, x)
  49.      int i, j, k;
  50.      double **x;
  51. {
  52.   int l;
  53.   double sum;
  54.  
  55.   for (l = i, sum = 0.0; l < n; l++) sum += x[l][j] * x[l][k];
  56.   return(sum);
  57. }
  58.  
  59. static void AXPY(i, j, k, n, a, x)
  60.      int i, j, k, n;
  61.      double a, **x;
  62. {
  63.   int l;
  64.   for (l = i; l < n; l++) x[l][k] = a * x[l][j] + x[l][k];
  65. }
  66.  
  67. static void SCALE(i, j, n, a, x)
  68.      int i, j, n;
  69.      double a, **x;
  70. {
  71.   int k;
  72.   for (k = i; k < n; k++) x[k][j] *= a;
  73. }
  74.  
  75. static void SWAP(i, j, n, a)
  76.   int i, j, n;
  77.   double **a;
  78. {
  79.   int k;
  80.   double temp;
  81.   for (k = 0; k < n; k++) {
  82.     temp = a[k][i];
  83.     a[k][i] = a[k][j];
  84.     a[k][j] = temp;
  85.   }
  86. }
  87.  
  88. void qrdecomp(x,n,p,v,jpvt,pivot)
  89.      int n, p, pivot;
  90.      /* int * */ IVector jpvt; /* changed JKL */
  91.      /* double ** */ RMatrix x, v;
  92. {
  93.   int i,j,k,jp,l,lp1,lup,maxj;
  94.   double maxnrm,tt,*qraux,*work;
  95.   double nrmxl,t;
  96.  
  97.   if (n < 0) return;
  98.   work = v[0];
  99.   qraux = rvector(p);
  100.  
  101.   /*
  102.    * compute the norms of the free columns.
  103.    */
  104.   if (pivot)
  105.     for (j = 0; j < p; j++) {
  106.       jpvt[j] = j;
  107.       qraux[j] = NORM2(0, j, n, x);
  108.       work[j] = qraux[j];
  109.     }
  110.   /*
  111.    * perform the householder reduction of x.
  112.    */
  113.   lup = (n < p) ? n : p;
  114.   for (l = 0; l < lup; l++) {
  115.     if (pivot) {
  116.       /*
  117.        * locate the column of largest norm and bring it
  118.        * into the pivot position.
  119.        */
  120.       maxnrm = 0.0;
  121.       maxj = l;
  122.       for (j = l; j < p; j++) 
  123.     if (qraux[j]>maxnrm) {
  124.       maxnrm = qraux[j];
  125.       maxj = j;
  126.     }
  127.       if (maxj!=l) {
  128.     SWAP(l, maxj, n, x);
  129.     qraux[maxj] = qraux[l];
  130.     work[maxj] = work[l];
  131.     jp = jpvt[maxj];
  132.     jpvt[maxj] = jpvt[l];
  133.     jpvt[l] = jp;
  134.       }
  135.     }
  136.     qraux[l] = 0.0;
  137.     if (l != n-1) {
  138.       /*
  139.        * compute the householder transformation for column l.
  140.        */
  141.       nrmxl = NORM2(l, l, n, x);
  142.       if (nrmxl != 0.0) {
  143.     if (x[l][l] != 0.0)
  144.       nrmxl = SIGN(nrmxl, x[l][l]);
  145.     SCALE(l, l, n, 1.0 / nrmxl, x);
  146.     x[l][l] = 1.0+x[l][l];
  147.     /*
  148.      * apply the transformation to the remaining columns,
  149.      * updating the norms.
  150.      */
  151.     lp1 = l+1;
  152.     for (j = lp1; j < p; j++) {
  153.       t = -DOT(l, l, j, n, x) / x[l][l];
  154.       AXPY(l, l, j, n, t, x);
  155.       if (pivot)
  156.         if (qraux[j]!=0.0) {
  157.           tt = 1.0-(fabs(x[l][j])/qraux[j])*(fabs(x[l][j])/qraux[j]);
  158.           if (tt < 0.0) tt = 0.0;
  159.           t = tt;
  160.           tt = 1.0+0.05*tt*(qraux[j]/work[j])*(qraux[j]/work[j]);
  161.           if (tt!=1.0)
  162.         qraux[j] = qraux[j]*sqrt(t);
  163.           else {
  164.         qraux[j] = NORM2(l+1, j, n, x);
  165.         work[j] = qraux[j];
  166.           }
  167.         }
  168.     }
  169.     /*
  170.      * save the transformation.
  171.      */
  172.     qraux[l] = x[l][l];
  173.     x[l][l] = -nrmxl;
  174.       }
  175.     }
  176.   }
  177.  
  178.   /* copy over the upper triangle of a */
  179.   for (i = 0; i < p; i++) {
  180.     for (j = 0; j < i; j++) v[i][j] = 0.0;
  181.     for (j = i; j < p; j++) {
  182.       v[i][j] = x[i][j];
  183.     }
  184.   }
  185.     
  186.   /* accumulate the Q transformation -- assumes p <= n */
  187.   for (i = 0; i < p; i++) {
  188.     x[i][i] = qraux[i];
  189.     for (k = 0; k < i; k++) x[k][i] = 0.0;
  190.   }
  191.   for (i = p - 1; i >= 0; i--) {
  192.     if (i == n - 1) x[i][i] = 1.0;
  193.     else {
  194.       for (k = i; k < n; k++) x[k][i] = -x[k][i];
  195.       x[i][i] += 1.0;
  196.     }
  197.     for (j = i - 1; j >= 0; j--) {
  198.       if (x[j][j] != 0.0) {
  199.     t = -DOT(j, j, i, n, x) / x[j][j];
  200.     AXPY(j, j, i, n, t, x);
  201.       }
  202.     }
  203.   }
  204.  
  205.   free_vector((Vector)qraux); /* cast added  JKL */
  206. }
  207.