bearing skrev:Precis, resultatet blev ju exakt som jag skrev.
Är det inte meningen att man ska skriva in förväntade resultat någonstans i det där programmet Valgrind?
Blir ju inte mycket till verifiering annars. Det blir som att en biltillverkares testförare kör bilen från 0-100 utan att mäta tiden tiden för accelerationen och jämföra med kravet.
Vadå exakt som du skrev?
Jag har skrivit flera hundra rader C kod och Valgrind kunde inte hitta ett enda fel på koden, alltså minnesläcka i stacken.
Däremot om jag medvetet gjorde felet så säger den till.
Här är ett exempel på One-Side SVD med Jacobi rotation. Vad kan man använda denna till då?
Jo, har man mycket mätdata som man samlar i en kvadratisk matris så kan man göra regression, systemidentifiering eller klassificering.
Har man en icke-kvadratisk matris så måste man gå på en QR-metod. Men QR-metoden är dyr och tar mycket tid. Då är det bättre med Golub och Reins algoritm, vilket är de fakto standard. One-Side SVD algorithmen är extremt snabb, men passar bara kvadratisk matriser.
Kod: Markera allt
void svd_jacobi_one_sided(float *A, int row, float *U, float *S, float *V) {
// i and j are the indices of the point we've chosen to zero out
float al, b, c, l, t, cs, sn, tmp, sign;
int i, j, p, k;
// Create the identity matrix
memset(U, 0, row*row*sizeof(float));
memset(V, 0, row*row*sizeof(float));
memset(S, 0, row*sizeof(float));
for(i = 0; i < row; i++){
*(U + row*i + i) = 1;
*(V + row*i + i) = 1;
}
// Apply MAX_ITERATION_COUNT_SVD times. That should be good enough
for (p = 0; p < MAX_ITERATION_COUNT_SVD; p++) {
// For all pairs i < j
for (i = 0; i < row; i++) {
for (j = i + 1; j < row; j++) {
al = b = c = l = t = cs = sn = tmp = sign = 0.0;
// Find the 2x2 submatrix
for (k = 0; k < row; k++) {
al += *(A + row*k + i) * *(A + row*k + i);
b += *(A + row*k + j) * *(A + row*k + j);
c += *(A + row*k + i) * *(A + row*k + j);
}
// Compute Jacobi rotation
l = (b - al) / (2.0 * c);
sign = 1.0;
if (l < 0.0)
sign = -1.0;
t = sign / ((sign * l) + sqrtf(1.0 + l * l));
cs = 1.0 / sqrtf(1.0 + t * t);
sn = cs * t;
// Change columns i and j only
for (k = 0; k < row; k++) {
tmp = *(A + row*k + i);
*(A + row*k + i) = cs * tmp - sn * *(A + row*k + j);
*(A + row*k + j) = sn * tmp + cs * *(A + row*k + j);
}
// Update the right singular vectors
for (k = 0; k < row; k++) {
tmp = *(V + row*k + i);
*(V + row*k + i) = cs * tmp - sn * *(V + row*k + j);
*(V + row*k + j) = sn * tmp + cs * *(V + row*k + j);
}
}
}
}
// Find the singular values by adding up squares of each column, then taking square root of each column
for (j = 0; j < row; j++) {
for (i = 0; i < row; i++) {
*(S + j) += *(A + row*i + j) * *(A + row*i + j);
}
tmp = *(S + j);
*(S + j) = sqrtf(tmp);
}
// Sort the singular values largest to smallest, and the right matrix accordingly
for (p = 0; p < (row - 1); p++) {
for (j = 0; j < row - p - 1; j++) {
if (*(S + j) < *(S + j + 1)) {
tmp = *(S + j);
*(S + j) = *(S + j + 1);
*(S + j + 1) = tmp;
// Rearrange columns of u and v accordingly
for (i = 0; i < row; i++) {
tmp = *(V + row*i + j);
*(V + row*i + j) = *(V + row*i + j + 1);
*(V + row*i + j + 1) = tmp;
tmp = *(A + row*i + j);
*(A + row*i + j) = *(A + row*i + j + 1);
*(A + row*i + j + 1) = tmp;
}
}
}
}
// A is A*V, so in order to get U, we divide by S in each column
for (i = 0; i < row; i++) {
for (j = 0; j < row; j++) {
*(A + row*i + j) = *(A + row*i + j) / *(S + j);
}
}
// Set U to A, since we have been making modifications to A
memcpy(U, A, row*row*sizeof(float));
}
/*
* GNU Octave code:
* >> A
A =
3 4 5
2 5 6
5 6 7
>> [u, s, v] = svd(A)
u =
-0.47384 0.10757 -0.87402
-0.53340 -0.82478 0.18767
-0.70069 0.55512 0.44819
s =
Diagonal Matrix
14.91526 0 0
0 1.58310 0
0 0 0.16940
v =
-0.401719 0.915134 -0.033930
-0.587752 -0.229240 0.775885
-0.702260 -0.331630 -0.629962
>>
*/
Jag har aldrig använt det verktyget förut eller ens kollat om jag har skrivit rätt kod.
