procedure square(x: ary2;
		 y: ary;
	     var a: ary2s;
	     var g: arys;
	 nrow,ncol: integer);

{ matrix multiplication routine }
{ a= transpose x times x }
{ g= y times x }

var
 i,k,l	: integer;

begin		{ square }
  for k:=1 to ncol do
    begin
	for l:=1 to k do
	  begin
	    a[k,l]:=0.0;
	    for i:=1 to nrow do
	    begin
		  a[k,l]:=a[k,l]+x[i,l]*x[i,k];
		  if k<>l then a[l,k]:=a[k,l]
	    end
	 end;		{ l-loop }
	 g[k]:=0.0;
	 for i:=1 to nrow do
	 g[k]:=g[k]+y[i]*x[i,k]
  end	{ k-loop }
end;	{ square }