QR-разложение для вектора в C - PullRequest
Новая
       19

QR-разложение для вектора в C

0 голосов
/ 29 октября 2018

Сначала я написал функцию для разложения QR, используя алгоритм Householder, и она работала для случайных матриц. Но сейчас я пытаюсь заставить мою функцию для вычисления Q * b работать. Алгоритм по сути тот же, но теперь мы умножаем на вектор b в конце, а не на матрицу A. Тем не менее, после 8 часов работы я не смог понять, что не так в моем коде. Я только начинающий в Си и буду очень признателен за некоторые подсказки. 

void qrhh2(int m, int n, double **A, double *b) { // Replace input vector b with Q* b, Q in A = QR
    double **V; // Initialize outer product matrix V and W for W = VA
    double *e; // Init basis vector

    for (int k=0; k<m; k++) { // Create operable column vector
        double *x = calloc(m-k, sizeof(double));
        for (int l=k; l<m; l++) {
            x[l-k] = A[l][k];
        }

        double *v = calloc(m-k, sizeof(double)); // Fill v with 0s
        e = calloc(m-k, sizeof(double *)); // Fill basis vector e with 0s
        e[0] = 1; // Complete making basis vector e
        for (int j=k; j<m; j++) {
            v[j-k] = sign(x[0]) * norm(m-k, x)*e[j-k] + x[j-k];
        }

        normalize(m-k, v);

        /* Fill V, W with 0s, initialize matrices V, W */
        V = calloc(m-k, sizeof(double *));
        for (int i=k; i<m; i++) V[i-k] = malloc((m-k) * sizeof(double));
        for (int i=k; i<m; i++) { // Create matrix V by outer product v*v^*
                for (int j=k; j<m; j++) {
                V[i-k][j-k] = v[i-k] * v[j-k];
            }

        }

        double *c = malloc((m-k)*sizeof(double));
        for (int i=0; i<m-k; i++) {
            c[i] = b[k+i];
        }
        double *r = mat_vec_mult(m-k, m-k, V, c);


        for (int i=k; i<m; i++) { // Create R from A=QR
                b[i] = b[i] -  2*r[i-k];
            }
        free(V);
    }
}

Пожалуйста, дайте мне знать, если я также должен опубликовать функции, используемые в этом коде, такие как нормализация, а также код, который я использую для его тестирования. Я получаю правильную первую запись нового вектора b, но не остальные записи.

...