qr.c

来自「混沌分析工具」· C语言 代码 · 共 54 行

void qcalc(B,b,Qout)
double **B,**Qout;
int b;
{
      double *y,sum,*norm;
      int loop,column,row;

      y = dvector(1,b);
      norm = dvector(1,b);

      /* first column is just divided by normalisation thing */

      sum = 0.0;
      for (loop=1;loop<=b;++loop) sum+=B[loop][1]*B[loop][1];
      sum = pow(sum,0.5);
      for (loop=1;loop<=b;++loop) Qout[loop][1]=B[loop][1]/sum;

      for (column=2;column<=b;++column){
	    for (loop=1;loop<column;++loop){
		  norm[loop]=0.0;
		  for (row=1;row<=b;++row) norm[loop]+=Qout[row][loop]*B[row][column];
	    }
	    for (row=1;row<=b;++row){
		  sum=0.0;
		  for (loop=1;loop<column;++loop) sum+=norm[loop]*Qout[row][loop];
		  y[row]=B[row][column]-sum;
	    }

	    sum=0.0;
	    for (row=1;row<=b;++row) sum+=y[row]*y[row];
	    sum = pow(sum,0.5);

	    for (row=1;row<=b;++row) Qout[row][column]=y[row]/sum;
      }
}


void qrd (A,a,Q,R)
double **A,**Q,**R;
int a;
{
      double **Qinv;

      Qinv = dmatrix(1,a,1,a);

      qcalc(A,a,Q);

      transpose(Q,Qinv,a,a);

      multiply_matrix(R,Qinv,A,a,a,a,a);
}