Householder Qr Zerlegung

[Thorsten Ransch]

Hi!
Ich habe hier einen Code für eine QR Householder Zerlegung. Er funktioniert auch. Ich vestehe jedoch die scale Variable nicht.
Für was braucht man die hier oder besser, was macht sie hier? (Zeile 115)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define SIGN(a,b) ((b) > 0.0 ? fabs(a) : - fabs(a))

static double maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1 = (a),maxarg2 = (b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))

static double sqrarg;
#define SQR(a) ((sqrarg = (a)) == 0.0 ? 0.0 : sqrarg * sqrarg)

//Global Variables
double **a; //Matrix whose SVD needs to be found
double *c;
double *d;

//Function
int qrdcmp(double **a, int m, int n, double *c, double *d);

int main (void)
{
int m, n;
int i, j;

printf("Enter the parameters for a square matrix(A)\n");
scanf("%d%d", &m, &n);

//Assigning A matrix
a = malloc(sizeof (double*)*m); //allocate M rows to dmatrix

for (i =0; i < m; i++)
{a[i] = malloc (sizeof (double)*n);} // Now for each row, allocate N actual floats

// From this step we have a matrix of M rows and N columns worth of floats
printf("Enter the elements of an array\n");
//2D array in effect
for (i =0; i < m; i++)
{
for (j =0; j<n; j++)
{
scanf("%lf", &a[i][j]);
}
}

//Assigning w matrix n

printf("\n");
printf("A matrix\n");
for (i =0; i < m; i++)
{
for (j =0; j< n; j++)
{
printf("%9.4lf\t",a[i][j]);
}
printf("\n");
}

c = malloc(sizeof(double)*n);
d = malloc(sizeof(double)*n);

for (i=0; i<n; i++)
{
c[i] = 0.0;
d[i] = 0.0;
}

qrdcmp (a,m,n,c,d);
double R[m][n];
for (i =0; i < m; i++)
{ for (j=0; j <n; j++)
{
if (i == j)
{R[i][i] = d[i];}

else if ( i < j)
{R[i][j] = a[i][j];}

else
{R[i][j] = 0;}
}
}

printf("\n");
printf("R decomposition of a matrix:\n");
for (i =0; i <m; i++)
{
for (j =0; j<n; j++)
{
printf("%9.4lf\t",R[i][j]);
}
printf("\n");
}

getch();
}

int qrdcmp(double **a, int m, int n, double *c, double *d)
/*Constructs the QR decomposition of a[1..n][1..n]. The upper triangular matrix R is returned in the upper triangle of a, except for the diagonal elements of R which are returned in
d[1..n]. The orthogonal matrix Q is represented as a product of n - 1 Householder matrices
Q1
. . . Qn-1
, where Qj = 1 - uj ? uj /cj . The ith component of uj is zero for i = 1, . . . , j - 1
while the nonzero components are returned in a[i][j] for i = j, . . . , n. sing returns as
true (1) if singularity is encountered during the decomposition, but the decomposition is still
completed in this case; otherwise it returns false (0).*/
{
int i,j,k;
double scale,sigma,sum,tau;
scale = sigma = sum = tau =0.0;
i=j=k=0;

int min;
if (m >= n) {min = n;}
else {min = m;}

for (k =0;k < min;k++) {
//printf("K =%d \n N = %d\n", k, n);

//for (i=k;i<n;i++) scale=FMAX(scale,fabs(a[i][k]));
for (i=k;i<m;i++) scale=FMAX(scale,fabs(a[i][k]));
//printf ("scale %lf\n", scale);

if (scale == 0.0) { //Singular case.
c[k]=d[k]=0.0;
}

else
{ //Form Q_k and Q_k * A.
//for (i=k;i<n;i++)
for (i=k;i<m;i++)
{a[i][k] /= scale;}
//printf ("a[%d][%d] = %lf\n",i,k, a[i][k]);}

//orgfor (sum=0.0,i=k;i<n;i++)
for (sum=0.0,i=k;i<m;i++)
{sum += SQR(a[i][k]);}
//printf ("Sum = %lf\n",sum);

sigma=SIGN(sqrt(sum),a[k][k]);
//printf ("sigma = %lf\n", sigma);

a[k][k] += sigma;
//printf ("a[%d][%d]= %lf\n",k ,k, a[k][k]);

c[k]=sigma*a[k][k];
//printf ("c[%d] = %lf\n", k,c[k]);

d[k] = -scale*sigma;
printf ("d[%d] = %lf\n", k,d[k]);

for (j=k+1;j < n;j++) {
//printf("j=%d\n",j);

//orgfor (sum=0.0,i=k;i<n;i++)
for (sum=0.0,i=k;i<m;i++)
{sum += a[i][k]*a[i][j];
//printf("sum = %lf, a[%d][%d] = %lf, a[%d][%d]=%lf\n",sum, i, k, a[i][k], i, j, a[i][j]);
}
//printf("sum = %lf\n",sum);

tau=sum/c[k];
//printf("tau= %lf\n",tau);

//org for (i=k;i<n;i++)
for (i=k;i<m;i++)
{
a[i][j] -= tau*a[i][k];
// printf("a[%d][%d] = %lf\n", i, j, a[i][j]);
}

}
}
}

return 0;
}