static char Sccsid[] = "a8.c @(#)a8.c 1.1 10/1/82 Berkeley "; #include "apl.h" ex_mdom() { register data *dp; register a; int i, j; struct item *p, *q; p = fetch1(); if(p->rank != 2) error("mdom C"); a = p->dim[0]; q = newdat(DA, 2, a*a); q->dim[0] = a; q->dim[1] = a; *sp++ = q; dp = q->datap; for(i=0; i<a; i++) for(j=0; j<a; j++) { datum = zero; if(i == j) datum = one; *dp++ = datum; } ex_ddom(); } ex_ddom() { struct item *p, *q; register a, b; register data *d1; int m, n, o; data *dmn, *dn1, *dn2, *dn3, *dm; int *in; char *al; p = fetch2(); q = sp[-2]; if(p->type != DA || q->type != DA) error("domino T"); if((p->rank != 1 && p->rank != 2) || q->rank != 2) error("domino C"); m = q->dim[0]; n = q->dim[1]; if(m < n || m != p->dim[0]) error("domino R"); o = 1; if(p->rank == 2) o = p->dim[1]; al = alloc(n*(SINT + SDAT*m + SDAT*3) + SDAT*m); if(al == 0) error("domino M"); dmn = (data *)al; dn1 = dmn + m*n; dn2 = dn1 + n; dn3 = dn2 + n; dm = dn3 + n; in = (int *)(dm + m); d1 = q->datap; for(b=0; b<m; b++) for(a=0; a<n; a++) dmn[a*m+b] = *d1++; a = lsq(dmn, dn1, dn2, dn3, dm, in, m, n, o, p->datap, q->datap); free(dmn); if(a) error("domino D"); sp--; pop(); *sp++ = p; p->dim[0] = n; p->size = n*o; } lsq(dmn, dn1, dn2, dn3, dm, in, m, n, p, d1, d2) data *dmn, *dn1, *dn2, *dn3, *dm; data *d1, *d2; int *in; { register data *dp1, *dp2; double f1, f2, f3, f4; int i, j, k, l; dp1 = dmn; dp2 = dn1; for(j=0; j<n; j++) { f1 = 0.; for(i=0; i<m; i++) { f2 = *dp1++; f1 += f2*f2; } *dp2++ = f1; in[j] = j; } for(k=0; k<n; k++) { f1 = dn1[k]; l = k; dp1 = dn1 + k+1; for(j=k+1; j<n; j++) if(f1 < *dp1++) { f1 = dp1[-1]; l = j; } if(l != k) { i = in[k]; in[k] = in[l]; in[l] = i; dn1[l] = dn1[k]; dn1[k] = f1; dp1 = dmn + k*m; dp2 = dmn + l*m; for(i=0; i<m; i++) { f1 = *dp1; *dp1++ = *dp2; *dp2++ = f1; } } f1 = 0.; dp1 = dmn + k*m + k; f3 = *dp1; for(i=k; i<m; i++) { f2 = *dp1++; f1 += f2*f2; } if(f1 == 0.) return(1); f2 = sqrt(f1); if(f3 >= 0.) f2 = -f2; dn2[k] = f2; f1 = 1./(f1 - f3*f2); dmn[k*m+k] = f3 - f2; for(j=k+1; j<n; j++) { f2 = 0.; dp1 = dmn + k*m + k; dp2 = dmn + j*m + k; for(i=k; i<m; i++) f2 += *dp1++ * *dp2++; dn3[j] = f1*f2; } for(j=k+1; j<n; j++) { dp1 = dmn + k*m + k; dp2 = dmn + j*m + k; f1 = dn3[j]; for(i=k; i<m; i++) *dp2++ -= *dp1++ * f1; f1 = dmn[j*m + k]; dn1[j] -= f1; } } for(k=0; k<p; k++) { dp1 = dm; l = k; for(i=0; i<m; i++) { *dp1++ = d1[l]; l += p; } solve(m, n, dmn, dn2, in, dm, dn3); l = k; dp1 = dm; for(i=0; i<m; i++) { f1 = -d1[l]; l += p; dp2 = dn3; for(j=0; j<n; j++) f1 += d2[i*n+j] * *dp2++; *dp1++ = f1; } solve(m, n, dmn, dn2, in, dm, dn1); f4 = 0.; f3 = 0.; dp1 = dn3; dp2 = dn1; for(i=0; i<n; i++) { f1 = *dp1++; f4 += f1*f1; f1 = *dp2++; f3 += f1*f1; } if(f3 > 0.0625*f4) return(1); loop: if(intflg) return(1); dp1 = dn3; dp2 = dn1; for(i=0; i<n; i++) *dp1++ += *dp2++; if(f3 <= 4.81e-35*f4) goto out; l = k; dp1 = dm; for(i=0; i<m; i++) { f1 = -d1[l]; l += p; dp2 = dn3; for(j=0; j<n; j++) f1 += d2[i*n+j] * *dp2++; *dp1++ = f1; } solve(m, n, dmn, dn2, in, dm, dn1); f2 = f3; f3 = 0.; dp1 = dn1; for(i=0; i<n; i++) { f1 = *dp1++; f3 += f1*f1; } if(f3 < 0.0625*f2) goto loop; out: l = k; dp1 = dn3; for(i=0; i<n; i++) { d1[l] = *dp1++; l += p; } } return(0); } solve(m, n, dmn, dn2, in, dm, dn1) data *dmn, *dn1, *dn2, *dm; int *in; { register data *dp1, *dp2; int i, j, k; double f1, f2; for(j=0; j<n; j++) { f1 = 0.; dp1 = dmn + j*m + j; dp2 = dm + j; f2 = *dp1; for(i=j; i<m; i++) f1 += *dp1++ * *dp2++; f1 /= f2*dn2[j]; dp1 = dmn + j*m + j; dp2 = dm + j; for(i=j; i<m; i++) *dp2++ += f1 * *dp1++; } dp1 = dm + n; dp2 = dn2 + n; dn1[in[n-1]] = *--dp1 / *--dp2; for(i=n-2; i>=0; i--) { f1 = -*--dp1; k = (i+1)*m + i; for(j=i+1; j<n; j++) { f1 += dmn[k] * dn1[in[j]]; k += m; } dn1[in[i]] = -f1 / *--dp2; } }