/* vector processing routines

from Handbook of C Tools for  Scientists and Engineers by L. Baker

often simplified versions of BLAS routines for vectors
with contiguous storage

dot(a,b,n)		returns double value of dot product of two vectors a,b
				of n elements
pv(v,n)			prints vector of n elements

mv(a,x,y,m,n)	y=ax where a is m x n matrix
normv(v,n)		normalize a vector (scale its length to one)
sqnor(v,n)		square of the length of vector v of n elements 
mvt(a,x,y,m,n)	like mv except y=a^x a^ =transpose of a
resid(r,a,x,n)	r=ax-x, a matrix all other vectors.

vs(v,s,n)		scale vector by multiplying each element by s
vset(v,s,n)		set vector v to scalar value s for each element
vcopy(x,y,n)	y=x, vectors 
vv(a,b,c,s,n)	a=b+c*s, s scalar, a,b,c,vectors 

*/

#include "ftoc.h"

double dot(a,b,n) int n; float a[],b[];
{
int i;
double sum;
if(n<=0)return(0.);
sum=0.;
DOFOR(i,n)sum+=a[i]*b[i];
return(sum);
}

pv(v,n) int n;float v[];
{
int btm,top,i;
btm=0;
top=0;
while (btm<n)
	{
	top=min(btm+6,n);
	printf(" printing vector from %d to %d\n",btm,(top-1));
	for(i=btm;i<top;i++)printf(" %e",v[i]);
	printf("\n");
	btm+=6;
	}
return;
}

mv(a,x,y,m,coln) int coln,m;
float x[],y[],a[];
/* y=ax a(m,n) m rows n columns [row of length n]*/
{
int i,j,k;
float sum;
DOFOR(i,m)
	{
	sum=0.;
	DOFOR(j,coln) sum+= x[j]*a INDEX(i,j);
	y[i]=sum;
	}
return;
}

normv(v,n) int n; float v[];
{
double x,sqrt(),sqnor();
x=sqnor(v,n);
if(x!=0.)x=1./sqrt(x);
vs(v,x,n);
return;
}

double sqnor(x,n) float x[]; int n;
{
int i;
double ans;
ans=0.;
if(n<=0)return(ans);
DOFOR(i,n)ans+= x[i]*x[i];
return(ans);
}

mvt(a,x,y,m,coln) int coln,m;
float a[],x[],y[];
{
/* y= a^x  a(m,n) a m rows n columns [n elements/row]
a^ n rows m columns*/
float sum;
int i,j;
DOFOR(i,coln)
	{
	sum=0.;
	DOFOR(j,m)sum+= x[j]*a INDEX(j,i);
	y[i]=sum;
	}
return;
}

resid(a,x,y,r,n) int n;
float a[],x[],y[],r[];
{
int i;
mv(a,x,r,n,n);
DOFOR(i,n)r[i]-=y[i];
return;
}

vs(v,s,n) int n; float v[],s;
{
int i;
DOFOR(i,n)v[i]*=s;
return;
}

vset(x,s,n) int n; float s,x[];
{
int i;
DOFOR(i,n)x[i]=s;
return;
}


vcopy(x,y,n) int n; float x[],y[];
{
int i;
DOFOR(i,n)y[i]=x[i];
return;
}

vv(a,b,c,s,n) int n; float a[],b[],c[],s;
{
int i;
DOFOR(i,n) a[i]=b[i]+s*c[i];
return;
}