OPTIONS LS=75 PS=60 NODATE;   
                              
DATA DIAG;                    
                              
INPUT  n   p  ajj;            
h2=0.25;                       
r=.5;                         
                              

di=(   ( ( 1+(n-1)*r) /n  )   + (p -1) * ajj * h2 ) / p;   
                             
*INPUT n   p  ajj    rel $;  
CARDS;
         2  20   0.25
         1  30   0.25
         1 135   0.25
         3   1   1.00
;
PROC PRINT; RUN;         
 
                             
                             
DATA OD;                     
INPUT G1 G2 G3 G4;           
CARDS; 
        0      0.25   0.125   0.25
        0.25   0      0.125   0.25
        0.125  0.125  0       0
        0.25   0.25   0       0
;
PROC IML; START;
h2=0.25;
P={17852 18751 18000 20000};   
PMEAN={14000 14000 14000 14000};
P=P-PMEAN;          
                                                                                                                                                                
aij={0.5,   0.5,   0.25,    0.5};   
RHS=aij*h2;
USE DIAG;
READ all      into DII; print DII;
di = DII (|,NCOL(DII) |);
DI = DIAG(di`);
PRINT di;                  
                                                                                                              
USE OD;
READ all      into OD;
OD=OD*h2;    PRINT OD;     
VP=OD+DI;    PRINT VP;     
b=INV(VP)*RHS;
PRINT b;
INDEX=P*b;
PRINT INDEX;
r_IH= SQRT(b`*aij);PRINT r_IH; 
FINISH; RUN;