Return to a8.c CVS log

Up to [CSRG BSD Unix] / 43BSD / contrib / apl / src

Annotation of 43BSD/contrib/apl/src/a8.c, revision 1.1

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.