Householder Qr Zerlegung



  • Hi! Ich habe hier einen Code für die QR Zerlegung.
    Er läuft und soweit ist er mir klar bis auf eine Sache:
    In der qrdcmp Funktion habe ich eine scale Variable (Zeile 115). Kann mir wer sagen, warum man sie einführt oder was die bringen soll?

    #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;
    }
    

    Vielen Dank soweit,
    Thorsten



  • Wäre das nicht besser im Unterforum "Mathematik und Physik" aufgehoben?

    Und nur weil C keine Einrückung verlangt, heißt es nicht, dass man darauf verzichten soll.


Log in to reply