|
|
1.1 ! root 1: static char Sccsid[] = "a7.c @(#)a7.c 1.1 10/1/82 Berkeley "; ! 2: #include "apl.h" ! 3: ! 4: ex_iprod() ! 5: { ! 6: register i, j; ! 7: struct item *p, *q, *r; ! 8: int param[10], ipr1(); ! 9: data (*fn)(); ! 10: ! 11: param[0] = exop[*pcp++]; ! 12: param[1] = exop[*pcp++]; ! 13: p = fetch2(); ! 14: q = sp[-2]; ! 15: if(p->type != DA || q->type != DA) ! 16: error("iprod T"); ! 17: /* ! 18: * extend scalars to match corresponding arg ! 19: */ ! 20: if(scalar(p)) { ! 21: if(scalar(q)){ ! 22: r = newdat(DA, 0, 1); ! 23: fn = param[1]; ! 24: r->datap[0] = (*fn)(p->datap[0], q->datap[0]); ! 25: goto out; ! 26: } ! 27: r = extend(DA, q->dim[0], p->datap[0]); ! 28: pop(); ! 29: *sp++ = p = r; ! 30: } ! 31: if(scalar(q)){ ! 32: r = extend(DA,p->dim[p->rank - 1], q->datap[0]); ! 33: free(sp[-2]); ! 34: sp[-2] = q = r; ! 35: } ! 36: bidx(p); ! 37: idx.rank--; ! 38: param[2] = idx.dim[idx.rank]; ! 39: if((param[2] != q->dim[0])) ! 40: /* && (param[2] != 1) */ ! 41: /* && (q->dim[0] != 1) */ ! 42: error("inner prod C"); ! 43: param[3] = q->size/param[2]; ! 44: for(i=1; i<q->rank; i++) ! 45: idx.dim[idx.rank++] = q->dim[i]; ! 46: r = newdat(DA, idx.rank, size()); ! 47: copy(IN, idx.dim, r->dim, idx.rank); ! 48: param[4] = 0; ! 49: param[5] = 0; ! 50: param[6] = p->datap; ! 51: param[7] = q->datap; ! 52: param[8] = r->datap; ! 53: param[9] = p->size; ! 54: forloop(ipr1, param); ! 55: out: ! 56: pop(); ! 57: pop(); ! 58: /* ! 59: * KLUDGE (we need the dim[0]'s for above stuff to work) ! 60: */ ! 61: if(r->rank == 1 && r->size == 1) ! 62: r->rank = 0; ! 63: *sp++ = r; ! 64: } ! 65: ! 66: ipr1(param) ! 67: int param[]; ! 68: { ! 69: register i, dk; ! 70: int lk, a, b; ! 71: data *dp1, *dp2, *dp3; ! 72: data (*f1)(), (*f2)(), d; ! 73: ! 74: f1 = param[0]; ! 75: f2 = param[1]; ! 76: dk = param[2]; ! 77: lk = param[3]; ! 78: a = param[4]; ! 79: b = param[5]; ! 80: dp1 = param[6]; ! 81: dp2 = param[7]; ! 82: dp3 = param[8]; ! 83: a += dk; ! 84: b += (dk * lk); ! 85: for(i=0; i<dk; i++) { ! 86: a--; ! 87: b -= lk; ! 88: d = (*f2)(dp1[a], dp2[b]); ! 89: if(i == 0) ! 90: datum = d; else ! 91: datum = (*f1)(d, datum); ! 92: } ! 93: *dp3++ = datum; ! 94: param[8] = dp3; ! 95: param[5]++; ! 96: if(param[5] >= lk) { ! 97: param[5] = 0; ! 98: param[4] += dk; ! 99: if(param[4] >= param[9]) ! 100: param[4] = 0; ! 101: } ! 102: } ! 103: ! 104: ex_oprod() ! 105: { ! 106: register i, j; ! 107: register data *dp; ! 108: struct item *p, *q, *r; ! 109: data *dp1, *dp2; ! 110: data (*f)(); ! 111: ! 112: f = (data *)exop[*pcp++]; ! 113: p = fetch2(); ! 114: q = sp[-2]; ! 115: if(p->type != DA || q->type != DA) ! 116: error("oprod T"); ! 117: /* ! 118: * collapse 1 element vectors to scalars ! 119: * ! 120: if(scalar(p)) ! 121: p->rank = 0; ! 122: if(scalar(q)) ! 123: q->rank = 0; ! 124: */ ! 125: bidx(p); ! 126: for(i=0; i<q->rank; i++) ! 127: idx.dim[idx.rank++] = q->dim[i]; ! 128: r = newdat(DA, idx.rank, size()); ! 129: copy(IN, idx.dim, r->dim, idx.rank); ! 130: dp = r->datap; ! 131: dp1 = p->datap; ! 132: for(i=0; i<p->size; i++) { ! 133: datum = *dp1++; ! 134: dp2 = q->datap; ! 135: for(j=0; j<q->size; j++) ! 136: *dp++ = (*f)(datum, *dp2++); ! 137: } ! 138: pop(); ! 139: pop(); ! 140: *sp++ = r; ! 141: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.