下面是一段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>