应热心回帖同学要求，我把之前发的SAS代码的等同C语言代码发上来，大家再帮我想想法子，已经试验过，SAS和C的这两段代码的运行结果是一摸一样的。
之前的SAS代码：
do j=1 to k;
s0=vk[j];
zvz=z[,j]*vi*z[,j]; //z[,j]表示矩阵的第j列，z[,j]表示z[,j]的转置
yvz=r*vi*z[,j]; dum=max((yvz**2-zvz)/zvz**2+s0,1e-30); vk[j]=dum; v = v + z[,j]*z[,j]*(vk[j]-s0);
vi= invupdt(vi,z[,j],vk[j]-s0);
vi=inv(v); //加是屏蔽代码，因为这句代码跟它前面一句在这里是一样作用的代码，但是效率低很多
end;
相应的C代码：
for(int j=0;j<XCOL;j++) //这句代码对应SAS里面的 do j=1 to k;XCOL的值就相当于k
{
for(int ione=0;ione<XROWXROW;ione++) //这里的XCOL和XROW已经#define过了，就是一个矩阵的行和列，为已知值
{
VItemp1[ione]=VI[ione]; //VI[XROWXROW]在这段代码之前已经附了初值
}
s0=VK[j]; //对应SAS里面的 s0=vk[j];
double ZVZ=0.0;
double YVZ=0.0;
double ZTVI[XROW]; //ZTVI[XROW]相当于SAS里面的z[,j] for(int i=0;i<XROW;i++) { ZTVI[i]=0.0; } for(int i=0;i<XROW;i++) { for(int k=0;k<XROW;k++) { ZTVI[i]=ZTVI[i]+XMatrix[j+k*XCOL]*VI[i+k*XROW]; //相当于SAS里面的z[,j]*vi
}
}
for(int i=0;i<XROW;i++)
{
ZVZ=ZVZ+ZTVI[i]XMatrix[j+iXCOL]; //相当于SAS里面的zvz
}
for(int i=0;i<XROW;i++)
{
RTVI[i]=0.0; //相当于SAS里面的r } for(int i=0;i<XROW;i++) { for(int k=0;k<XROW;k++) { RTVI[i]=RTVI[i]+R[k]*VI[i+k*XROW]; //相当于SAS里面的r*vi
}
}
for(int i=0;i<XROW;i++)
{
YVZ=YVZ+RTVI[i]XMatrix[j+iXCOL]; //相当于SAS里面的yvz
}
double dum=0.0;
dum=Maxmium((YVZYVZ-ZVZ)/(ZVZZVZ)+s0,1e-30); //之前已经定义过一个求两个数中最大值的函数
VK[j]=dum;
for(int i=0;i<XROW;i++)
{
for(int k=0;k<XROW;k++)
{
V[i*XROW+k]=V[i*XROW+k]+XMatrix[i*XCOL+j]XMatrix[kXCOL+j](VK[j]-s0);
//相当于SAS里面的v = v + z[,j]z[,j]`(vk[j]-s0);
}
}
for(int i=0;i<XROWXROW;i++)
{
VI[i]=V[i];
}
if((VK[j]-s0)!=0)
{
inverse(VI,XROW); //相当于SAS里面的vi=inv(v);这里调用了CULA的包
for(int ione=0;ione<XROWXROW;ione++)
{
VItemp1[ione]=VI[ione];
}
}
else
{
for(int ione=0;ione<XROWXROW;ione++) // 这里的if()else()选择对应了SAS里面的vi= invupdt(vi,z[,j],vk[j]-s0);
{
VI[ione]=VItemp1[ione];
}
}
}