72 template<
class mat_t,
class mat2_t,
class vec_t,
class vec2_t>
74 mat_t &A, mat2_t &V, vec_t &S, vec2_t &work) {
77 O2SCL_ERR2(
"SVD decomposition of MxN matrix with M<N, ",
121 double fbm1=O2SCL_IX(work,b-1);
122 if (fbm1==0.0 || gsl_isnan(fbm1)) {
135 double fam1=O2SCL_IX(work,a-1);
136 if (fam1==0.0 || gsl_isnan(fam1)) {
146 O2SCL_ERR(
"SV decomposition failed to converge in SV_decomp().",
153 size_t n_block=b-a+1;
158 for(
size_t i=0;i<n_block;i++) {
159 double aa=fabs(O2SCL_IX(S,a+i));
160 if (aa>norm) norm=aa;
164 for(
size_t i=0;i<n_block-1;i++) {
165 double aa=fabs(O2SCL_IX(work,a+i));
166 if (aa>norm) norm=aa;
172 if (norm>GSL_SQRT_DBL_MAX) {
173 scale=norm/GSL_SQRT_DBL_MAX;
175 }
else if (norm<GSL_SQRT_DBL_MIN && norm>0.0) {
176 scale=norm/GSL_SQRT_DBL_MIN;
234 for (
size_t j=0;j<K;j++) {
235 double Sj=O2SCL_IX(S,j);
237 for (
size_t i=0;i<N;i++) {
238 O2SCL_IX2(V,i,j)=-O2SCL_IX2(V,i,j);
247 for (
size_t i=0;i<K;i++) {
249 double S_max=O2SCL_IX(S,i);
252 for (
size_t j=i+1;j<K;j++) {
253 double Sj=O2SCL_IX(S,j);
262 o2scl::vector_swap<vec_t,double>(S,i,i_max);
265 o2scl::matrix_swap_cols<mat_t,double>(M,A,i,i_max);
266 o2scl::matrix_swap_cols<mat2_t,double>(M,V,i,i_max);
288 template<
class mat_t,
class mat2_t,
class mat3_t,
class vec_t,
290 (
size_t M,
size_t N, mat_t &A, mat2_t &X, mat3_t &V, vec_t &S,
296 O2SCL_IX2(V,0,0)=1.0;
309 for (
size_t i=0;i<N;i++) {
320 for (
size_t i=0;i<N;i++) {
321 for (
size_t j=0;j<i;j++) {
322 O2SCL_IX2(X,i,j)=0.0;
325 O2SCL_IX2(X,i,i)=O2SCL_IX2(A,i,i);
326 for (
size_t j=i+1;j<N;j++) {
327 O2SCL_IX2(X,i,j)=O2SCL_IX2(A,i,j);
333 for (
size_t j=N;j-- > 0;) {
335 double tj=O2SCL_IX(S,j);
345 for (
size_t i=0;i<M;i++) {
347 for(
size_t j=0;j<N;j++) {
348 O2SCL_IX(work,j)=0.0;
351 for (
size_t j=0;j<N;j++) {
352 double Lij=O2SCL_IX2(A,i,j);
356 for (
size_t j=0;j<N;j++) {
357 O2SCL_IX2(A,i,j)=O2SCL_IX(work,j);
372 template<
class mat_t,
class mat2_t,
class vec_t,
class vec2_t,
class vec3_t>
374 mat_t &U, mat2_t &V, vec_t &S, vec2_t &b, vec3_t &x) {
382 double *w=
new double[N];
385 (O2SCL_CBLAS_NAMESPACE::o2cblas_RowMajor,
386 O2SCL_CBLAS_NAMESPACE::o2cblas_Trans,M,N,1.0,U,b,0.0,w);
388 for (
size_t i=0;i<N;i++) {
389 double alpha=O2SCL_IX(S,i);
390 if (alpha != 0) alpha=1.0/alpha;
395 (O2SCL_CBLAS_NAMESPACE::o2cblas_RowMajor,
396 O2SCL_CBLAS_NAMESPACE::o2cblas_NoTrans,N,N,1.0,V,w,0.0,x);
437 template<
class mat_t,
class mat2_t,
class vec_t>
452 double dbl_eps=std::numeric_limits<double>::epsilon();
454 double tolerance=10*M*dbl_eps;
457 if (sweepmax<12) sweepmax=12;
460 for(
size_t i=0;i<N;i++) {
461 for(
size_t j=0;j<N;j++) {
462 if (i==j) O2SCL_IX2(Q,i,j)=1.0;
463 else O2SCL_IX2(Q,i,j)=0.0;
470 for (
size_t j=0;j<N;j++) {
472 O2SCL_IX(S,j)=dbl_eps*sj;
477 while (count > 0 && sweep <= sweepmax) {
482 for (
size_t j=0;j<N-1;j++) {
483 for (
size_t k=j+1;k<N;k++) {
491 double abserr_a, abserr_b;
492 int sorted, orthog, noisya, noisyb;
495 for(
size_t ii=0;ii<M;ii++) {
496 p+=O2SCL_IX2(A,ii,j)*O2SCL_IX2(A,ii,k);
509 abserr_a=O2SCL_IX(S,j);
510 abserr_b=O2SCL_IX(S,k);
515 orthog=(fabs (p) <= tolerance*a*b);
519 if (sorted && (orthog || noisya || noisyb)) {
525 if (v == 0 || !sorted) {
529 cosine=sqrt((v+q)/(2.0*v));
530 sine=p/(2.0*v*cosine);
534 for (
size_t i=0;i<M;i++) {
535 const double Aij=O2SCL_IX2(A,i,j);
536 const double Aik=O2SCL_IX2(A,i,k);
537 O2SCL_IX2(A,i,j)=Aij*cosine+Aik*sine;
538 O2SCL_IX2(A,i,k)=-Aij*sine+Aik*cosine;
541 O2SCL_IX(S,j)=fabs(cosine)*abserr_a+fabs(sine)*abserr_b;
542 O2SCL_IX(S,k)=fabs(sine)*abserr_a+fabs(cosine)*abserr_b;
545 for (
size_t i=0;i<N;i++) {
546 const double Qij=O2SCL_IX2(Q,i,j);
547 const double Qik=O2SCL_IX2(Q,i,k);
548 O2SCL_IX2(Q,i,j)=Qij*cosine+Qik*sine;
549 O2SCL_IX2(Q,i,k)=-Qij*sine+Qik*cosine;
561 double prev_norm=-1.0;
563 for (
size_t j=0;j<N;j++) {
571 if (norm == 0.0 || prev_norm == 0.0
572 || (j > 0 && norm <= tolerance*prev_norm)) {
577 for(
size_t ii=0;ii<M;ii++) {
578 O2SCL_IX2(A,ii,j)=0.0;
598 O2SCL_ERR(
"Jacobi iterations did not reach desired tolerance",
614 template<
class mat_t,
class vec_t>
617 for(
size_t j=0;j<N;j++) O2SCL_IX(D,j)=1.0;
619 for(
size_t j=0;j<N;j++) {
622 if (s==0.0 || !std::isfinite(s)) {
int bidiag_unpack2(size_t M, size_t N, mat_t &A, vec_t &tau_U, vec2_t &tau_V, mat2_t &V)
Unpack a matrix A with the bidiagonal decomposition and create matrix V.
void SV_solve(size_t M, size_t N, mat_t &U, mat2_t &V, vec_t &S, vec2_t &b, vec3_t &x)
Solve the system A x = b using the SV decomposition.
void balance_columns(size_t M, size_t N, mat_t &A, vec_t &D)
Balance a general matrix A by scaling the columns by the diagonal matrix D.
exceeded max number of iterations
failed to reach the specified tolerance
requested feature not (yet) implemented
double householder_transform_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute the Householder transform of a vector formed with n rows of a column of a matrix...
void dgemv(const enum o2cblas_order order, const enum o2cblas_transpose TransA, const size_t M, const size_t N, const double alpha, const mat_t &A, const vec_t &X, const double beta, vec2_t &Y)
Compute .
void householder_hm_subcol(mat_t &M, const size_t ir, const size_t ic, const size_t nr, const size_t nc, const mat_t &M2, const size_t ir2, const size_t ic2, double tau)
Apply a Householder transformation to the lower-right part of M when the transformation is stored in ...
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
double dasum_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute for a subcolumn of a matrix.
The namespace for linear algebra classes and functions.
double dnrm2_subcol(const mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute the norm of a subcolumn of a matrix.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void qrstep_sub(size_t M, size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a)
A special form of qrstep() for SV_decomp()
void chop_small_elements(size_t N, vec_t &d, vec2_t &f)
Zero out small elements in f according to the scales set in d.
void SV_decomp(size_t M, size_t N, mat_t &A, mat2_t &V, vec_t &S, vec2_t &work)
Factorise a general matrix into its SV decomposition using the Golub-Reinsch algorithm.
void SV_decomp_mod(size_t M, size_t N, mat_t &A, mat2_t &X, mat3_t &V, vec_t &S, vec2_t &work)
SV decomposition by the modified Golub-Reinsch algorithm which is better for .
int bidiag_decomp(size_t M, size_t N, mat_t &A, vec_t &tau_U, vec2_t &tau_V)
Factor a matrix into bidiagonal form.
void SV_decomp_jacobi(size_t M, size_t N, mat_t &A, mat2_t &Q, vec_t &S)
SV decomposition using one-sided Jacobi orthogonalization.
void daxpy_subrow(const double alpha, const size_t N, const mat_t &X, const size_t ir, const size_t ic, vec_t &Y)
Compute for a subrow of a matrix.
void dscal_subvec(const double alpha, const size_t N, vec_t &X, const size_t ie)
Compute beginning with index ie and ending with index N-1.
void dscal_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M, const double alpha)
Compute for a subcolumn of a matrix.
void householder_hm1_sub(const size_t M, const size_t N, double tau, mat_t &A, size_t ir, size_t ic)
Apply a Householder transformation to a matrix being build up from the identity matrix, using the first column of A as a Householder vector.