divqr.c

这是一个用于矩阵QR分解的程序

#include "divQR.h"
#include <string.h>
int qr(double *a,int Row,int Col,double *Rel_Q,double *Rel_R)//m行n列
{
    int i,j,k,jj;
	int temps,templ,tempp;
	double tempu,tempalpha,tempw,tempt;

	memcpy(Rel_R,a,Row*Col*8);

	if(Row<Col)
	{
		return 0;
	}
	//Rel_Q[]=I初值
	//memset(Rel_Q,0x00,Row*Row*8);//double数组 赋值函数?? 
	for(i=0;i<Row*Row;i++)
	{
		if(i%(Row+1))
		{
			Rel_Q[i]=0.0;
		}
		else
		{
			Rel_Q[i]=1.0;
		}
	}
	//s	
	temps=(Row==Col ? Row-1:Col);	////待考查
   //k=0---temps-1
	for(k=0;k<temps;k++)
	{
		tempu=0.0;
		templ=k*Col+k;
		for(i=k;i<Row;i++)
		{
		   tempw=fabs(Rel_R[i*Col+k]);
		   if(tempw>tempu) 
		   {
			   tempu=tempw;
		   }
		}

	///tempalpha
		tempalpha=0.0;
		for(i=k;i<Row;i++)
		{
			tempt=Rel_R[i*Col+k]/tempu;
			tempalpha=tempalpha+tempt*tempt;
		}
		if(Rel_R[templ]>0.0)
		{
			tempu=-tempu;
		}
		tempalpha=tempu*sqrt(tempalpha);

		if(fabs(tempalpha)+1.0==1.0)
		{
			 return 0;
		}
		tempu=sqrt(2.0*tempalpha*(tempalpha-Rel_R[templ]));
		if((tempu+1.0)!=1.0)
		{
			Rel_R[templ]=(Rel_R[templ]-tempalpha)/tempu;
			for(i=k+1;i<=Row-1;i++)
			{
				tempp=i*Col+k;
				Rel_R[tempp]=Rel_R[tempp]/tempu;
			}
			for(j=0;j<=Row-1;j++)
			{
				tempt=0.0;
				for(jj=k;jj<=Row-1;jj++)
				{
					tempt=tempt+Rel_R[jj*Col+k]*Rel_Q[jj*Row+j];
				}
				for(i=k;i<=Row-1;i++)
				{
					tempp=i*Row+j;
					Rel_Q[tempp]=Rel_Q[tempp]-2.0*tempt*Rel_R[i*Col+k];
				}
			}
			for(j=k+1;j<=Col-1;j++)
			{
				tempt=0.0;
				for(jj=k;jj<=Row-1;jj++)
				{
					tempt=tempt+Rel_R[jj*Col+k]*Rel_R[jj*Col+j];
				}
				for(i=k;i<=Row-1;i++)
				{
					tempp=i*Col+j;
					Rel_R[tempp]=Rel_R[tempp]-2.0*tempt*Rel_R[i*Col+k];
				}
			}
			Rel_R[templ]=tempalpha;
			for(i=k+1;i<=Row-1;i++)
			{
				Rel_R[i*Col+k]=0.0;
			}
				
		}

	}
///Q转置
	for(i=0;i<=Row-2;i++)
	{
		for(j=i+1;j<=Row-1;j++)
		{
			tempp=i*Row+j;
			templ=j*Row+i;
			tempt=Rel_Q[tempp];
		    Rel_Q[tempp]=Rel_Q[templ];
			Rel_Q[templ]=tempt;
		}
	}
	return 1;
}