divqr1.c

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

#include "divQR.h"
int qr(double *a,int m,int n,double *q)//m行n列
{
    int i,j,k,l,nn,p,jj;
	double u,alpha,w,t;

	if(m<n)
	{
		return 0;
	}
	//q[]=I初值
	for(i=0;i<=m-1;i++)
	{
	  for(j=0;j<=m-1;j++)
	  {
		  l=i*m+j;
		  q[l]=0.0;
		  if(i==j)
		  {
			 q[l]=1.0;
		  }
	  }
	}
	//s	
	nn=n;
	if(m==n)
	{
	   nn=m-1;
	}
   //k=0---nn-1
	for(k=0;k<=nn-1;k++)
	{
	   u=0.0;
	   l=k*n+k;
	   for(i=k;i<=m-1;i++)
	   {
		   w=fabs(a[i*n+k]);
		   if(w>u) 
		   {
			   u=w;
		   }
	   }

	///alpha
	   alpha=0.0;
	   for(i=k;i<=m-1;i++)
	   {
	      t=a[i*n+k]/u;
		  alpha=alpha+t*t;
		}
		if(a[l]>0.0)
		{
			u=-u;
		}
		alpha=u*sqrt(alpha);

		if(fabs(alpha)+1.0==1.0)
		{
			 return 0;
		}
		u=sqrt(2.0*alpha*(alpha-a[l]));
		if((u+1.0)!=1.0)
		{
			a[l]=(a[l]-alpha)/u;
			for(i=k+1;i<=m-1;i++)
			{
				p=i*n+k;
				a[p]=a[p]/u;
			}
			for(j=0;j<=m-1;j++)
			{
				t=0.0;
				for(jj=k;jj<=m-1;jj++)
				{
					t=t+a[jj*n+k]*q[jj*m+j];
				}
				for(i=k;i<=m-1;i++)
				{
					p=i*m+j;
					q[p]=q[p]-2.0*t*a[i*n+k];
				}
			}
			for(j=k+1;j<=n-1;j++)
			{
				t=0.0;
				for(jj=k;jj<=m-1;jj++)
				{
					t=t+a[jj*n+k]*a[jj*n+j];
				}
				for(i=k;i<=m-1;i++)
				{
					p=i*n+j;
					a[p]=a[p]-2.0*t*a[i*n+k];
				}
			}
			a[l]=alpha;
			for(i=k+1;i<=m-1;i++)
			{
				a[i*n+k]=0.0;
			}
				
		}

	}

///Q转置
	for(i=0;i<=m-2;i++)
	{
		for(j=i+1;j<=m-1;j++)
		{
			p=i*m+j;
			l=j*m+i;
			t=q[p];
			q[p]=q[l];
			q[l]=t;
		}
	}
	return 1;
}