QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

tbeyn does not crash anymore but the result is wrong 

Former-commit-id: 3288074674f1a3f61b43d905a96c865fa32faccd
This commit is contained in:
Marek Nečada 2019-08-23 12:20:51 +03:00
parent 26e7d9acee
commit e898fd5ad8
1 changed files with 12 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
qpms/beyn.c
View File

@ -156,8 +156,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function,
//A0->SVD(Sigma, &V0Full, &W0TFull); //A0->SVD(Sigma, &V0Full, &W0TFull);
QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
// FIXME not supported by GSL; use LAPACKE_zgesdd QPMS_ENSURE(V0Full->size1 >= V0Full->size2, "m = %zd, l = %zd, what the hell?");
//gsl_linalg_complex_SV_decomp(V0Full, W0TFull, Sigma, work);
QPMS_ENSURE_SUCCESS(LAPACKE_zgesdd(LAPACK_ROW_MAJOR, // A = U*Σ*conjg(V') QPMS_ENSURE_SUCCESS(LAPACKE_zgesdd(LAPACK_ROW_MAJOR, // A = U*Σ*conjg(V')
'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/, 'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/,
V0Full->size1 /* m, number of rows */, V0Full->size1 /* m, number of rows */,
@ -166,7 +165,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function,
V0Full->tda /*lda*/, V0Full->tda /*lda*/,
Sigma->data /*s*/, Sigma->data /*s*/,
NULL /*u; not used*/, NULL /*u; not used*/,
0 /*ldu; not used*/, M /*ldu; should not be really used but lapacke requires it for some obscure reason*/,
(lapack_complex_double *)W0TFull->data /*vt*/, (lapack_complex_double *)W0TFull->data /*vt*/,
W0TFull->tda /*ldvt*/ W0TFull->tda /*ldvt*/
)); ));
@ -186,24 +185,21 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function,
} }
// set V0, W0T = matrices of first K right, left singular vectors // set V0, W0T = matrices of first K right, left singular vectors
gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(M,K); gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(M,K);
gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,L); gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,L);
gsl_vector_complex *TempM = gsl_vector_complex_calloc(M); for (int k = 0; k < K; ++k) {
gsl_vector_complex *TempL = gsl_vector_complex_calloc(L); gsl_vector_complex_view tmp;
for(int k=0; k<K; k++){ tmp = gsl_matrix_complex_column(V0Full, k);
// It should be rows and cols like this, right..? gsl_matrix_complex_set_col(V0, k, &(tmp.vector));
gsl_matrix_complex_get_row(TempM, V0Full,k); tmp = gsl_matrix_complex_row(W0TFull, k);
gsl_matrix_complex_set_row(V0, k, TempM); gsl_matrix_complex_set_row(W0T, k, &(tmp.vector));
gsl_matrix_complex_get_col(TempL,W0TFull,k);
gsl_matrix_complex_set_col(W0T,k, TempL);
} }
gsl_matrix_complex_free(V0Full); gsl_matrix_complex_free(V0Full);
gsl_matrix_complex_free(W0TFull); gsl_matrix_complex_free(W0TFull);
gsl_vector_complex_free(TempM);
gsl_vector_complex_free(TempL);
// B <- V0' * A1 * W0 * Sigma^-1 // B <- V0' * A1 * W0 * Sigma^-1
@ -262,7 +258,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function,
S->tda /* lda */, S->tda /* lda */,
(lapack_complex_double *) Lambda->data /* w */, (lapack_complex_double *) Lambda->data /* w */,
NULL /* vl */, NULL /* vl */,
0 /* ldvl */, M /* ldvl, not used by for some reason needed */,
(lapack_complex_double *)(B->data)/* vr */, (lapack_complex_double *)(B->data)/* vr */,
B->tda/* ldvr */ B->tda/* ldvr */
)); ));