DSPF_dp_qrd_solver_cmplx

**

本博客目的是解释QR分解求解线性方程组（复数形式）的源码

**
一：储备知识
（1）该函数是在完成矩阵QR分解（householder法）的基础之上进行的。
故double *restrict Q和double *restrict R都是已知的，double *restrict b也是已知的（restrict用于 声明只会通过该指针对其指向的内存空间进行读写操作，以便编译器能进行更好的优化。）
（2）共轭转置矩阵
即先取共轭（把虚部取反），后进行转置。

(有个小错误，应该是4+2i)

（3）复数的有关运算
即复数的乘除法（略）

（4）高斯消元法
高斯消元法详解

简述就是：
将Ax=b按照从上至下、从左至右的顺序化为上三角方程组，中间过程不对矩阵进行交换。
（存在的限制详见上述博客）

二：过程简述

谭博宇同学的QR分解求解线性方程组（实数形式）
三：源码解释

（1）：求y（y=Q*b）

//（restrict用于 声明只会通过该指针对其指向的内存空间进行读写操作，以便编译器能进行更好的优化。）
int DSPF_dp_qrd_solver_cmplx(const int Nrows,const int Ncols,double *restrict Q,double *restrict R,double *restrict b,double *restrict y,double *restrict x) {

  short row,col,loop_cnt,Nrows2,Ncols2;//
  double xreal,ximag,yreal,yimag,zreal,zimag,sum_real,sum_imag;
//  _nassert()：简单判断
  _nassert(Nrows>0);
  _nassert(Ncols>0);
  _nassert((int)Q % 8 == 0);
  _nassert((int)R % 8 == 0);
  _nassert((int)b % 8 == 0);
  _nassert((int)y % 8 == 0);
  _nassert((int)x % 8 == 0);
  //第一步：求y
  /* generate y=Q'*b */
  Nrows2=2*Nrows;//存实部和虚部
  Ncols2=2*Ncols;
//把#pragma MUST_ITERATE(, , )放在循环体之前，告知开发板循环次数，改善软件流水
//MUST_ITERATE告诉编译器循环的属性，但是这些属性必须是真实的，不然程序可能运行出错。此指令主要用于优化C函数循环，一般情况下，只要有循环都最好带上此指令
//
//#pragma MUST_ITERATE(min, max, multiple)；其中multiple参数必须有，循环执行次数必是multiple的整数倍。
//这个信息对编译器使用软件流水技术非常重要
#pragma MUST_ITERATE(1,,)
  for (row=0;row<Nrows;row++) {
	sum_real=0;
	sum_imag=0;
#pragma MUST_ITERATE(1,,)
	for (col=0;col<Nrows;col++) 
  {
	  xreal= Q[2*row  +col*Nrows2];
	  ximag=-Q[2*row+1+col*Nrows2];
	  yreal=b[2*col  ];
	  yimag=b[2*col+1];
	  sum_real+=xreal*yreal-ximag*yimag;
	  sum_imag+=xreal*yimag+ximag*yreal;
	}
	y[2*row  ]=sum_real;
	y[2*row+1]=sum_imag;
  }

（2）：求x（高斯消元法）

这部分和实数形式如出一辙，看代码就很好理解
complex_dp_div函数也就是进行复数的除法运算

/* use backward substitution to solve x=inv(R)*y */


//loop_cnt取两者较小值
//但是如果Nrows小的话说明可能有无穷多解那么就把(Ncols-Nrows)这几个多出来的x值取0就行了
  if (Nrows>=Ncols) 
  	loop_cnt=Ncols;
  else             
  	loop_cnt=Nrows;
  	
   //memset():常用于较大的对结构体或数组的清零操作。
  memset(x,0,sizeof(double)*Ncols2);
  
#pragma MUST_ITERATE(1,,)
  for (row=loop_cnt-1;row>=0;row--)
   {
	sum_real=0;
	sum_imag=0;
#pragma MUST_ITERATE(1,,)
	for (col=row+1;col<loop_cnt;col++)
	 {
	  xreal=R[2*col  +row*Ncols2];
	  ximag=R[2*col+1+row*Ncols2];
	  yreal=x[2*col  ];
	  yimag=x[2*col+1];
	  //当然一开始sum_real,sum_imag都为零，后面求出xn_real,xn_imag的值了便不为零（理解高斯消元法就很方便了）
	  sum_real+=xreal*yreal-ximag*yimag;
	  sum_imag+=xreal*yimag+ximag*yreal;
	}
	xreal=y[2*row  ]-sum_real;
	ximag=y[2*row+1]-sum_imag;
	yreal=R[2*row  +row*2*Ncols];
	yimag=R[2*row+1+row*2*Ncols];
	//复数的除法运算
	complex_dp_div(xreal,ximag,yreal,yimag,&zreal,&zimag);
	x[2*row  ]=zreal;
	x[2*row+1]=zimag;
  }

  return 0;
}

static void complex_dp_div(double x_real,double x_imag,double y_real,double y_imag,double *z_real,double *z_imag) 
{
//#define ENABLE_NR 1
//条件编译的用法：略
  double inv_mag_sq;
#ifdef ENABLE_NR
  double x,y;
#endif
#ifndef ENABLE_NR
  inv_mag_sq=1/(y_real*y_real+y_imag*y_imag);
#else
  x=y_real*y_real+y_imag*y_imag;
  y=_rcpdp(x);//取倒数
  //多次操作用以提高精度（项目对精度要求很高）
  y=y*(2.0-x*y);
  y=y*(2.0-x*y);
  y=y*(2.0-x*y);
  inv_mag_sq=y;
#endif
  *z_real=(x_real*y_real+x_imag*y_imag)*inv_mag_sq;
  *z_imag=(x_imag*y_real-x_real*y_imag)*inv_mag_sq;

}

下面举个小李子：

①：

sum_real=0;
	sum_imag=0;
#pragma MUST_ITERATE(1,,)
	for (col=row+1;col<loop_cnt;col++)
	 {
	  xreal=R[2*col  +row*Ncols2];
	  ximag=R[2*col+1+row*Ncols2];
	  yreal=x[2*col  ];
	  yimag=x[2*col+1];
	  //当然一开始sum_real,sum_imag都为零，后面求出xn_real,xn_imag的值了便不为零（理解高斯消元法就很方便了）
	  sum_real+=xreal*yreal-ximag*yimag;
	  sum_imag+=xreal*yimag+ximag*yreal;

②：

xreal=y[2*row  ]-sum_real;
ximag=y[2*row+1]-sum_imag;

③：

complex_dp_div(xreal,ximag,yreal,yimag,&zreal,&zimag);
	x[2*row  ]=zreal;
	x[2*row+1]=zimag;