home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 13 / MA_Cover_13.bin / source / c / apl / a8.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-11-23  |  4.0 KB  |  234 lines
  1. #include "apl.h"
  2.  
  3. ex_mdom()
  4. {
  5.     data *dp;
  6.     int a, i, j;
  7.     struct item *p, *q;
  8.  
  9.     p = fetch1();
  10.     if(p->rank != 2) error("mdom C");
  11.     a = p->dim[0];
  12.     q = newdat(DA, 2, a*a);
  13.     q->dim[0] = a;
  14.     q->dim[1] = a;
  15.     *sp++ = q;
  16.     dp = q->datap;
  17.     for(i=0; i<a; i++) {
  18.         for(j=0; j<a; j++) {
  19.             datum = zero;
  20.             if(i == j) datum = one;
  21.             *dp++ = datum;
  22.         }
  23.     }
  24.     ex_ddom();
  25. }
  26.  
  27. ex_ddom()
  28. {
  29.     struct item *p, *q;
  30.     int a, b, m, n, o, *in;
  31.     data *d1, *dmn, *dn1, *dn2, *dn3, *dm;
  32.     char *al;
  33.  
  34.     p = fetch2();
  35.     q = sp[-2];
  36.     if(p->type != DA || q->type != DA) error("domino T");
  37.     if((p->rank != 1 && p->rank != 2) || q->rank != 2) error("domino C");
  38.     m = q->dim[0];
  39.     n = q->dim[1];
  40.     if(m < n || m != p->dim[0]) error("domino R");
  41.     o = 1;
  42.     if(p->rank == 2) o = p->dim[1];
  43.     al = alloc(n*(SINT + SDAT*m + SDAT*3) + SDAT*m);
  44.     if(al == 0) error("domino M");
  45.     dmn = (data *)al;
  46.     dn1 = dmn + m*n;
  47.     dn2 = dn1 + n;
  48.     dn3 = dn2 + n;
  49.     dm = dn3 + n;
  50.     in = (int *)(dm + m);
  51.     d1 = q->datap;
  52.     for(b=0; b<m; b++) {
  53.         for(a=0; a<n; a++) dmn[a*m+b] = *d1++;
  54.     }
  55.     a = lsq(dmn, dn1, dn2, dn3, dm, in, m, n, o, p->datap, q->datap);
  56.     aplfree(dmn);
  57.     if(a) error("domino D");
  58.     sp--;
  59.     pop();
  60.     *sp++ = p;
  61.     p->dim[0] = n;
  62.     p->size = n*o;
  63. }
  64.  
  65. lsq(dmn, dn1, dn2, dn3, dm, in, m, n, p, d1, d2)
  66. data *dmn, *dn1, *dn2, *dn3, *dm;
  67. data *d1, *d2;
  68. int *in;
  69. {
  70.     data *dp1, *dp2;
  71.     double f1, f2, f3, f4;
  72.     int i, j, k, l;
  73.  
  74.     dp1 = dmn;
  75.     dp2 = dn1;
  76.     for(j=0; j<n; j++) {
  77.         f1 = 0.;
  78.         for(i=0; i<m; i++) {
  79.             f2 = *dp1++;
  80.             f1 += f2*f2;
  81.         }
  82.         *dp2++ = f1;
  83.         in[j] = j;
  84.     }
  85.     for(k=0; k<n; k++) {
  86.         f1 = dn1[k];
  87.         l = k;
  88.         dp1 = dn1 + k+1;
  89.         for(j=k+1; j<n; j++) {
  90.             if(f1 < *dp1++) {
  91.                 f1 = dp1[-1];
  92.                 l = j;
  93.             }
  94.         }
  95.         if(l != k) {
  96.             i = in[k];
  97.             in[k] = in[l];
  98.             in[l] = i;
  99.             dn1[l] = dn1[k];
  100.             dn1[k] = f1;
  101.             dp1 = dmn + k*m;
  102.             dp2 = dmn + l*m;
  103.             for(i=0; i<m; i++) {
  104.                 f1 = *dp1;
  105.                 *dp1++ = *dp2;
  106.                 *dp2++ = f1;
  107.             }
  108.         }
  109.         f1 = 0.;
  110.         dp1 = dmn + k*m + k;
  111.         f3 = *dp1;
  112.         for(i=k; i<m; i++) {
  113.             f2 = *dp1++;
  114.             f1 += f2*f2;
  115.         }
  116.         if(f1 == 0.) return(1);
  117.         f2 = sqrt(f1);
  118.         if(f3 >= 0.) f2 = -f2;
  119.         dn2[k] = f2;
  120.         f1 = 1./(f1 - f3*f2);
  121.         dmn[k*m+k] = f3 - f2;
  122.         for(j=k+1; j<n; j++) {
  123.             f2 = 0.;
  124.             dp1 = dmn + k*m + k;
  125.             dp2 = dmn + j*m + k;
  126.             for(i=k; i<m; i++) f2 += *dp1++ * *dp2++;
  127.             dn3[j] = f1*f2;
  128.         }
  129.         for(j=k+1; j<n; j++) {
  130.             dp1 = dmn + k*m + k;
  131.             dp2 = dmn + j*m + k;
  132.             f1 = dn3[j];
  133.             for(i=k; i<m; i++) *dp2++ -= *dp1++ * f1;
  134.             f1 = dmn[j*m + k];
  135.             dn1[j] -= f1;
  136.         }
  137.     }
  138.     for(k=0; k<p; k++) {
  139.         dp1 = dm;
  140.         l = k;
  141.         for(i=0; i<m; i++) {
  142.             *dp1++ = d1[l];
  143.             l += p;
  144.         }
  145.         solve(m, n, dmn, dn2, in, dm, dn3);
  146.         l = k;
  147.         dp1 = dm;
  148.         for(i=0; i<m; i++) {
  149.             f1 = -d1[l];
  150.             l += p;
  151.             dp2 = dn3;
  152.             for(j=0; j<n; j++) f1 += d2[i*n+j] * *dp2++;
  153.             *dp1++ = f1;
  154.         }
  155.         solve(m, n, dmn, dn2, in, dm, dn1);
  156.         f4 = 0.;
  157.         f3 = 0.;
  158.         dp1 = dn3;
  159.         dp2 = dn1;
  160.         for(i=0; i<n; i++) {
  161.             f1 = *dp1++;
  162.             f4 += f1*f1;
  163.             f1 = *dp2++;
  164.             f3 += f1*f1;
  165.         }
  166.         if(f3 > 0.0625*f4) return(1);
  167.     loop:
  168.         if(intflg) return(1);
  169.         dp1 = dn3;
  170.         dp2 = dn1;
  171.         for(i=0; i<n; i++) *dp1++ += *dp2++;
  172.         if(f3 <= 4.81e-35*f4) goto out;
  173.         l = k;
  174.         dp1 = dm;
  175.         for(i=0; i<m; i++) {
  176.             f1 = -d1[l];
  177.             l += p;
  178.             dp2 = dn3;
  179.             for(j=0; j<n; j++) f1 += d2[i*n+j] * *dp2++;
  180.             *dp1++ = f1;
  181.         }
  182.         solve(m, n, dmn, dn2, in, dm, dn1);
  183.         f2 = f3;
  184.         f3 = 0.;
  185.         dp1 = dn1;
  186.         for(i=0; i<n; i++) {
  187.             f1 = *dp1++;
  188.             f3 += f1*f1;
  189.         }
  190.         if(f3 < 0.0625*f2) goto loop;
  191.     out:
  192.         l = k;
  193.         dp1 = dn3;
  194.         for(i=0; i<n; i++) {
  195.             d1[l] = *dp1++;
  196.             l += p;
  197.         }
  198.     }
  199.     return(0);
  200. }
  201.  
  202. solve(m, n, dmn, dn2, in, dm, dn1)
  203. data *dmn, *dn1, *dn2, *dm;
  204. int *in;
  205. {
  206.     data *dp1, *dp2;
  207.     int i, j, k;
  208.     double f1, f2;
  209.  
  210.     for(j=0; j<n; j++) {
  211.         f1 = 0.;
  212.         dp1 = dmn + j*m + j;
  213.         dp2 = dm + j;
  214.         f2 = *dp1;
  215.         for(i=j; i<m; i++) f1 += *dp1++ * *dp2++;
  216.         f1 /= f2*dn2[j];
  217.         dp1 = dmn + j*m + j;
  218.         dp2 = dm + j;
  219.         for(i=j; i<m; i++) *dp2++ += f1 * *dp1++;
  220.     }
  221.     dp1 = dm + n;
  222.     dp2 = dn2 + n;
  223.     dn1[in[n-1]] = *--dp1 / *--dp2;
  224.     for(i=n-2; i>=0; i--) {
  225.         f1 = -*--dp1;
  226.         k = (i+1)*m + i;
  227.         for(j=i+1; j<n; j++) {
  228.             f1 += dmn[k] * dn1[in[j]];
  229.             k += m;
  230.         }
  231.         dn1[in[i]] = -f1 / *--dp2;
  232.     }
  233. }
  234.