Proc GLM v.s. IML in ANOVA

data thrust;
    input speed100 speed250 speed400 speed550;
    if _n_ le 15 then feed=0.015;
    if _n_ le 10 then feed=0.010;
    if _n_ le 5 then feed=0.005;
    cards;
     121 98 83 58
     124 108 81 59
     104 87 88 60
     124 94 90 66
     110 91 86 56
     329 291 281 265
     331 265 278 265
     324 295 275 269
     338 288 276 260
     332 297 287 251
     640 569 551 487
     600 575 552 481
     612 565 570 487
     620 573 546 500
     623 588 569 497
    ;
run;

data thrust2;
    set thrust;
    array speedarray(4) speed:;
    do i = 1 to 4;
       force = speedarray(i);
       output;
    end;
    drop speed:;
run;

data thrust3;
    set thrust2;
    if i = 1 then speed = 100;
    if i = 2 then speed = 250;
    if i = 3 then speed = 400;
    if i = 4 then speed = 550;
    drop i;
run;

proc sort data=thrust3 out=thrust4;
    by feed speed;
run;

proc glm data=thrust3;
    title 'GLM solution';
    class feed speed;
    model force= feed speed feed*speed;
run;

proc iml;
    title 'IML solution';
    use thrust4;
    read all var{force} into y;
    Xzero=J(60,1);
    Xfeed=I(3)@j(4,1)@j(5,1);
    Xspeed=j(3,1)@I(4)@j(5,1);
    Xint=I(3)@I(4)@j(5,1);
    Xfull=Xzero||Xfeed||Xspeed||Xint;

    SSE=y`*y-y`*Xfull*ginv(Xfull`*Xfull)*Xfull`*y;
    SSTrtFeed=y`*((I(3)-(1/3)*J(3,3))@((1/4)*J(4,4))@((1/5)*J(5,5)))*y;
    SSTrtSpeed=y`*(((1/3)*J(3,3))@(I(4)-(1/4)*J(4,4))@((1/5)*J(5,5)))*y;
    SSInt=y`*((I(3)-(1/3)*J(3,3))@(I(4)-(1/4)*J(4,4))@((1/5)*J(5,5)))*y;
    SSTotal=y`*(I(60)-(1/60)*J(60,60))*y;

    print SSTrtFeed SSTrtSpeed SSInt SSE SSTotal;

    MSTrtFeed=SSTrtFeed/2; MSTrtSpeed=SSTrtSpeed/3; MSInt=SSInt/6; MSE=SSE/48; 
    print MSTrtFeed MSTrtSpeed MSInt MSE ;

    F_TrtFeed=MSTrtFeed/MSE; F_TrtSpeed=MSTrtSpeed/MSE; F_Int=MSInt/MSE;
    Print  F_TrtFeed F_TrtSpeed F_Int;

    P_TrtFeed=1-probf(F_TrtFeed,2,48); 
    P_TrtSpeed=1-probf(F_TrtSpeed,3,48); 
    P_Int=1-probf(F_Int,6,48);
    Print P_TrtFeed P_TrtSpeed P_Int;
run;quit;