PROC IML;
                                    

rG={   1.00     0.433     -0.401,   
       0.433    1.00       0.529,
      -0.401    0.529       1.00};
                                    

rP={   1.00    0.416      0.058,
       0.416   1.00       0.104,
       0.058   0.104      1.00};

sP={4.17, 1.24, 23.58};              
H2={0.210, 0.378, 0.157};              
sG=SQRT(sP#H2#sP);

V0=Hdir(sP, I(3))* rP * Hdir(sP, I(3));  
G0=Hdir(sG, I(3))* rG * Hdir(sG, I(3));  
V00=V0;
G00=G0;
 

W= { 3.97,                              
    -6.36,
    1};

INV_V0= INV(V0);                          
b = INV(V0)*G0*w;   PRINT b;
                         

G0w=G0*w;                           PRINT G0w; 
V_I = b`*V0*b;                      PRINT V_I; 
V_H = w`*G0*w;                      PRINT V_H; 
R_IH  = SQRT (V_I/V_H);             PRINT R_IH;

DLT_A = G0*b*(1 / SQRT(V_I) );      PRINT DLT_A;
                                                

 DUMY=J(2,2,0);k2={0,0};  PRINT DUMY k2;  
                                          

 y=G0w//k2;          PRINT y;    
 G1=G0(| ,1|);       PRINT G1;   
 G2=G0(| ,2|);       PRINT G2;   
 G02=G1||G2;         PRINT G02;  
 

 C1=V0||G02;                     
 C2=G02`||DUMY;                  
 C =C1//C2;          PRINT C;    
 

 INV_C=INV(C);
 b0=INV_C*y;         PRINT INV_C;PRINT b0;  
V_I = b0[1:3,]`*V00*b0[1:3,];           PRINT V_I; 
DLT_A =G0*b0[1:3,]*(1 / SQRT(V_I) );      PRINT DLT_A;

 kc=G0*b;                             
 yc=G0w//kc;            PRINT kc yc;  
 

 O22=J(3,3,0);          PRINT O22;    
 CC1=V0||G0;
 CC2=G0||O22;                         
 CC =CC1//CC2;          PRINT CC;     
                                      

 INV_CC=INV(CC);
 b02=INV_CC*yc;         PRINT INV_CC;PRINT b02;   
                                

 
 ratio={0, 0, 3 };   
 kp=ratio#kc;                   
 PRINT ratio kc kp;
 yp=G0w//kp;            PRINT G0w; PRINT kp yp;
 

bd1=INV_CC*yp;                                           
bd2=INV(G0)*kp;        PRINT INV_CC;PRINT bd1 bd2;       
bd3=bd2*(0.099/ 6.1628537);   print bd3; * 解2，解3 成比例非唯一解;

RUN;