4.3BSD - /usr/src/new/apl/src/a8.c

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

Defined functions

Defined variables