/* Copyright (C) 1994 by Thomas Glen Smith.	All Rights Reserved. */
/* innrprdp APL2 V1.0.0 ************************************************
* Called by execdota, execdotc, and execdotd.					 *
* Differs from innrprd in that it expects opera and operb to be		 *
* procedure calls rather than functions.						 *
* Does inner product for all matching data types, or of mixed data	 *
* types where one is convertible to the other, i.e. innrprd will never *
* be invoked with either left or rite of type APLCHAR, and the other	 *
* something else.  If left is of type APLCHAR, rite will be also.	 *
***********************************************************************/
#define INCLUDES APLCB
#include "includes.h"
Aplcb innrprdp(opera,operb,identity,left,rite)
void (*opera)(),(*operb)(),*identity;
Aplcb left,rite;
{
	Dtacopy; Errstop; Innrcom; Matchok;
	int axc,dty,iw,jw,kw,lxc,lax,lbot,lincr,ltop,mw,nw,
		raxc,rbot,rincr,rtop,type;
	char *od,*ip,*jp,*kp,*ld,*mp,*np,*rd,*tp;
	Aplcb out;
	double wa[2],wb[2],wc[2];

	if (!matchok(&left,&rite,APLMASK))
		return(NULL);
	out = innrcom(1,left,rite,&lax,&lxc,&lbot,&ltop,&lincr,
		&raxc,&rbot,&rtop,&rincr,&dty,&od,&ld,&rd);
	if (out == NULL)
		return(NULL);
	axc = raxc;
	for (iw=0; iw<ltop; iw++) {
		ip = ld + iw*axc*lincr*left->aplsize;
		for (jw=0; jw<lbot; jw++) {
			jp = ip + jw*left->aplsize;
			for (kw=0; kw<rtop; kw++) {
				kp = rd + kw*axc*rincr*rite->aplsize;
				for (mw=0; mw<rbot; mw++) {
					mp = kp + (mw + axc*rincr)*rite->aplsize;
					np = jp + axc*lincr*left->aplsize;
					tp = dtacopy(wa, identity, 1, 1, dty);
					for (nw=0; nw<axc; nw++) {
						(*operb)(	np -= lincr*left->aplsize,
								mp -= rincr*rite->aplsize,
								wb );
						tp = dtacopy(wc, wa, 1, 1, dty);
						(*opera)(	wb, wc, wa);
					}
					od = dtacopy(od, wa, 1, 1, dty);
				}
			}
		}
	}
	return(errstop(0,left,rite,out));
}