![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to a8.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / contrib / apl / src|
Return to a8.c CVS log | Up to [CSRG BSD Unix] / 43BSD / contrib / apl / src |
1.1 root 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.