Här kommer jag att ha olika matriser på a, q och r. Men dessa matriser allokeras bara EN gång.
Kod: Markera allt
void cut_internal_qr(matrix* a, int start_n, int stop_n, int start_m, int stop_m, float* b, int out_columns);
float dot_internal_qr(float* a, float* b, int n);
float norm_internal_qr(float* a, int n);
/*===========================================================================
* qr
* Get Q and R matrix from matrix by using Modified Gram-Schmidt method
* Input: Matrix
* Return: void
* Works: OK
*=========================================================================*/
void qr(matrix* a, matrix* q, matrix* r) {
// Get info about matrix a
int n = a->row;
int m = a->column;
float* ptr = a->data;
// Get data
float* ptr_q = q->data;
float* ptr_r = r->data;
float c1[n];
float c2[n];
memset(c1, 0, sizeof(c1));
memset(c2, 0, sizeof(c2));
for(int k = 0; k < m; k++){
// Insert vector
for(int j = 0; j < n; j++)
*((ptr_q + j*n) + k) = *((ptr + j*n) + k);
// Do the dot product
for(int i = 0; i < k; i++){
cut_internal_qr(q, 0, n-1, i, i, c1, 1);
cut_internal_qr(q, 0, n-1, k, k, c2, 1);
*((ptr_r + i*m)+k) = dot_internal_qr(c1, c2, n);
// Insert vector again
for(int j = 0; j < n; j++)
*((ptr_q + j*n) + k) = *((ptr_q + j*n) + k) - *((ptr_r + i*m)+k)*(*((ptr_q + j*n) +i));
}
cut_internal_qr(q, 0, n-1, k, k, c1, 1);
*((ptr_r+k*m)+k) = -norm_internal_qr(c1, n); // Important with negative
// Insert vector again
for(int j = 0; j < n; j++){
if(*((ptr_r + k*m) + k) == 0){
*((ptr_r + k*m) + k) = pow(2.2204, -16); // Same as eps command in MATLAB
}
*((ptr_q + j*n) + k) = (*((ptr_q + j*n) + k)) / (*((ptr_r + k*m) + k));
}
}
}
void cut_internal_qr(matrix* a, int start_n, int stop_n, int start_m, int stop_m, float* b, int out_columns) {
int in_columns = a->column;
float* data = a->data + start_n * in_columns + start_m;
// Instead of having two for loops, we just copy the whole row at once.
for (int i = start_n; i < stop_n + 1; i++) {
memcpy(b, data, sizeof(float) * out_columns);
b += out_columns;
data += in_columns;
}
}
float dot_internal_qr(float* a, float* b, int n) {
// Reset
float sum = 0;
// Multiply each row
for (int i = 0; i < n; ++i) {
sum += (*(a + i)) * (*(b + i));
}
return sum;
}
float norm_internal_qr(float* a, int n) {
float sum = 0; // Initial
for (int i = 0; i < n; i++)
sum += (*(a + i)) * (*(a + i));
return sqrt(sum);
}
/* MATLAB CODE
* function [Q,R] = mgs(X)
% Modified Gram-Schmidt. [Q,R] = mgs(X);
% G. W. Stewart, "Matrix Algorithms, Volume 1", SIAM, 1998.
[n,p] = size(X);
Q = zeros(n,p);
R = zeros(p,p);
for k = 1:p
Q(:,k) = X(:,k);
for i = 1:k-1
R(i,k) = Q(:,i)'*Q(:,k);
Q(:,i)'*Q(:,k)
Q(:,k) = Q(:,k) - R(i,k)*Q(:,i);
end
R(k,k) = -norm(Q(:,k))';
Q(:,k) = Q(:,k)/R(k,k);
end
end
*/
