title;title2;
proc iml;
size = {100 200 400 600 800 1000};
levels = repeat(T(0:3), ncol(size));
sizes =  colvec(T(repeat(size, 4)));
results = levels || sizes || j(nrow(levels),1,.);
cnt = 1:4;

do i = 1 to ncol(size);
   X = j(size[i], size[i]);
   call randseed(1);
   call randgen(x, "Normal"); /* example data */

   n = nrow(X);
   p = ncol(X);
    
   t = time();
   S=j(p,p,.);            /* allocate matrix to hold results */
   do j = 1 to p;
      do k = 1 to p;
         u=0;     
         do m = 1 to n;
            u = u +  X[m,j]*X[m,k];
         end;
         S[j,k] = u;      /* the (j,k)th result */
      end;
   end;
   t0 = time()-t;         /* level-0 BLAS */

   /* Step 2: The sum of the product of two vectors is their dot product. */
   t = time();
   S=j(p,p,.);
   do j = 1 to p;
      do k = 1 to p;
         S[j,k] = X[ ,j]` * X[,k];
      end;
   end;
   t1 = time() - t;

   /* Step 3: Convert to vector-matrix product */
   t = time();
   S=j(p,p,.);
   do j=1 to p;
      S[j, ] = X[ ,j]` * X;
   end;
   t2 = time() - t;

   /* Step 4: A loop over all rows is equivalent to multiplication by X` */
   t = time();
   S = X` * X;
   t3 = time() - t;

   /* store results */
   results[cnt,3] = t0//t1//t2//t3;
   cnt = cnt + 4;
end;

print results;
create Results from results[c={"Level" "Size" "Time"}];
append from results;
close Results;

quit;

data All;
set Results(in=R0)
    Results(where=(Level>0) in=R1)
    Results(where=(Level>1) in=R2);
if R0 then PanelID = 1;
else if R1 then PanelID = 2;
else PanelID = 3;
run;

proc format;
 value Trial 1="All Levels"
             2="Zoom: Level 1 - Level 3"
             3="Zoom: Level 2 - Level 3";
run;

title "Comparison of Different Computational Methods";
proc sgpanel data=All;
format PanelID Trial.;
panelby PanelID / columns=1 novarname uniscale=column;
series x=Size y=Time / group=Level;
run;
    

title2 "Level 0 - Level 3";
proc sgplot data=Results;
series x=Size y=Time / group=Level;
run;

title2 "Level 1 - Level 3";
proc sgplot data=Results(where=(Level>0));
series x=Size y=Time / group=Level;
run;

title2 "Level 2 - Level 3";
proc sgplot data=Results(where=(Level>1));
series x=Size y=Time / group=Level;
run;

data LogResults;
set Results;
LogTime = log10(Time);
label LogTime="Log10(Time)";
CubeRtTime = Time**(1/3);
label CubeRtTime="(Time)^(1/3)";
run;

proc sgplot data=LogResults;
series x=Size y=CubeRtTime / group=Level;
run;