Nu körde jag SVD på en 36x36 matris....1.867489 sekunder.
Kod: Markera allt
#include <stdio.h>
#include <stdlib.h>
#include "LinearAlgebra/declareFunctions.h"
#include <time.h>
int main() {
/*
* G(s) = 3.2/(2s^2 + 0.7s + 3.1) - Model
*/
float input[72] = { 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5};
float output[72] = { 0.00000, 1.89649, 3.30777, 3.40331, 3.23615, 3.18151, 3.19049, 3.20013, 3.20132, 3.20031, 3.19990,
3.19993, 3.20000, 3.20001, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000,
3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000,
3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000,
3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000,
3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000,
3.20000, 3.20000, 3.20000, 3.20000, 3.20000, 3.20000};
// Create toeplitz matrix
matrix* toe = toeplitz(input, 72);
// Create upper triangular matrix of the toeplitz matrix
matrix* tru = triu(toe, 0);
// Find the inverse of tru
matrix* iv = inv(tru);
// Create vector the horizon - Important! Else, we cannot find the markov-parameters g
matrix* Y = create(output, 1, 72);
// Multiply Y with the iv matrix - Find the markov-parameters g
matrix* G = mul(Y, iv, false);
// Detete some mats
freeMatrix(toe);
freeMatrix(tru);
freeMatrix(iv);
freeMatrix(Y);
/*--------------------------------------*/
// turn vector G into a normal vector because the function hankel only want 1D vector
float g[72];
for(int i = 0; i < 72; i++)
g[i] = *(G->data + i);
// Create hankel matrix H0 and H1
matrix* H0 = hankel(g, 72, 0);
matrix* H1 = hankel(g, 72, 1);
// Cut H1 and H2 to the half - Remember indexing from zero here!
matrix* H0_half = cut(H0, 0, 35, 0, 35); // 36x36 from 72x72
matrix* H1_half = cut(H1, 0, 35, 0, 35);
// Remove H0
freeMatrix(H0);
// Measure the size
/*
int n;
int m;
size(H0_half, &n, &m);
printf("The H0_half have the size %dx%d\n\n,", n, m);
*/
// Do SVD on
matrix* u = initMatrix(36,36); // We know the size from the size() command above
matrix* s = initMatrix(36,36);
matrix* v = initMatrix(36,36);
clock_t start, end;
float cpu_time_used;
start = clock();
svd(H0_half, u, s, v);
end = clock();
cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
printf("Total speed of SVD was %f,", cpu_time_used);
//printMatrix(u);
return EXIT_SUCCESS;
}
Nu får vi inte glömma att GNU Octave gör sina tuffaste beräkningar i Fortran, vilket är snabbare än C.