下面是一段R程序。
<br />
tau<-function(tau,cx,cyx,wi,lambda,p,q)##tau和cyx是p*q的矩阵，cx是q*q的矩阵，wi是p*p的矩阵##<br />
{<br />
tau1<-tau<br />
tau2<-tau1+1<br />
while(sum(abs(tau1-tau2))>t2)<br />
{<br />
for(i in 1:p)<br />
{<br />
for(j in 1:q)<br />
{<br />
g<-2*((t(cyx)%*%wi)[j,i]+cx[j,j]*wi[i,i]*tau[i,j]-(cx%*%t(tau)%*%wi)[j,i])<br />
tau[i,j]=sign(g)*ifelse((abs(g)-lambda)>0,abs(g)-lambda,0)/(2*cx[j,j]*wi[i,i])</p>
<p> }<br />
}<br />
tau2<-tau<br />
count<-c(count,sum(abs(tau1-tau2)))<br />
tau1<-tau<br />
}<br />
tau<br />
}<br />
但是运行速度太慢，所以想用c编写，但是用c编写的话就要写一个矩阵的乘法函数，以及sign的函数。其实整个c大致已经写好了，但是我不知道该怎么组织，而且我想让p，q是变量。c如下：
<br />
#include<stdio.h><br />
#include<math.h><br />
//////////////////////////*************///////////////////////<br />
#define p 50<br />
#define q 50<br />
///////////////////////////////////////////////////////////</p>
<p>//matrix multiply<br />
void multi1(double a[p][q],double b[p][p],double result[q][p])<br />
{</p>
<p> int i,j,k;<br />
for(i=0;i<q;i++)<br />
{<br />
for(j=0;j<p;j++)<br />
{<br />
for(k=0;k<p;k++)<br />
{<br />
result[i][j]=result[i][j]+a[k][i]*b[k][j];<br />
}<br />
}<br />
}<br />
}<br />
double multi2(double a[q][q],double b[p][q], double c[p][p],int row, int col)<br />
{<br />
int i,j;<br />
double sum=0;<br />
double result[q][p];<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
result[i][j]=0;<br />
}<br />
}<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
result[row][i]=result[row][i]+a[row][j]*b[i][j];<br />
}<br />
}</p>
<p> for(i=0;i<p;i++)<br />
{<br />
sum=sum+result[row][i]*c[i][col];<br />
}<br />
return sum;</p>
<p>}<br />
//give the value in matrix a to matrix b<br />
void copy(double a[p][q],double b[p][q])<br />
{<br />
int i,j;<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
b[i][j]=a[i][j];<br />
}<br />
}<br />
}<br />
//sign function<br />
int sign(double g)<br />
{<br />
if(g>0)<br />
return 1;<br />
else<br />
if(g<0)<br />
return -1;<br />
else<br />
return 0;</p>
<p>}<br />
//max function<br />
double max(double a, double b)<br />
{<br />
return a>=b?a:b;<br />
}</p>
<p>void main()<br />
{<br />
//open file pointers---change the path!!!!<br />
//by the way the R write out need to be : write.table(tau,"C:\\tau.txt",row.names=F,col.names=F,sep=" ")<br />
/////////////////////*******************//////////////////////////////////<br />
FILE *fp1=fopen("c:\\tau.txt","r");<br />
FILE *fp2=fopen("c:\\cyx.txt","r");<br />
FILE *fp3=fopen("c:\\wi.txt","r");<br />
FILE *fp4=fopen("c:\\cx.txt","r");<br />
FILE *fp=fopen("c:\\tau_new.txt","w");<br />
////////////////////////////////////////////////////////////////////////////</p>
<p> double tau[p][q];<br />
double tau1[p][q];<br />
double cyx[p][q];<br />
double wi[p][p];<br />
double cx[q][q];</p>
<p>//中间变量<br />
double cyxwi[q][p];<br />
double tauwi[q][p];<br />
double cxtauwi[q][p];</p>
<p> double t2=0.0001;<br />
double flag=1;<br />
double lambda=0.001;<br />
double term=0;<br />
double g=0;</p>
<p> int i,j;<br />
int time=0;<br />
//read in data<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
fscanf(fp1,"%lf",&tau[i][j]);<br />
fscanf(fp2,"%lf",&cyx[i][j]);<br />
}<br />
}<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<p;j++)<br />
{<br />
fscanf(fp3,"%lf",&wi[i][j]);</p>
<p> }<br />
}<br />
for(i=0;i<q;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
fscanf(fp4,"%lf",&cx[i][j]);</p>
<p> }<br />
}</p>
<p>//initialize<br />
for (i=0;i<q;i++)<br />
{<br />
for(j=0;j<p;j++)<br />
{<br />
cyxwi[i][j]=0;<br />
tauwi[i][j]=0;<br />
cxtauwi[i][j]=0;<br />
}<br />
}</p>
<p> multi1(cyx,wi,cyxwi);</p>
<p> while(flag)<br />
{<br />
time++;<br />
//give tau value to tau1<br />
copy(tau,tau1);<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
cxtauwi[j][i]=multi2(cx,tau,wi,j,i);<br />
g=2*(cyxwi[j][i]+cx[j][j]*wi[i][i]*tau[i][j]-cxtauwi[j][i]);<br />
tau[i][j]=sign(g)*max(fabs(g)-lambda,0)/(2*cx[j][j]*wi[i][i]);<br />
}<br />
}</p>
<p> for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
term=term+fabs(tau[i][j]-tau1[i][j]);<br />
}<br />
}<br />
if(term<=t2)<br />
flag=0;<br />
term=0;<br />
}<br />
/* print the tau matrix result to the screen<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
printf("%lf ",tau[i][j]);<br />
}<br />
printf("\n");<br />
}<br />
//printf loop time<br />
printf("%d",time);<br />
*/<br />
//write out the new tau matrix value<br />
for(i=0;i<p;i++)<br />
{<br />
for(j=0;j<q;j++)<br />
{<br />
fprintf(fp,"%lf",tau[i][j]);<br />
}<br />
fprintf(fp,"\n");<br />
}</p>
<p>}<br />
</p>
