svd_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2010-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* linalg/svd.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Gerard
26  * Jungman, Brian Gough
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 /** \file svd_base.h
44  \brief File for SVD composition
45 */
46 
47 #ifdef DOXYGEN
48 namespace o2scl_linalg {
49 #endif
50 
51  /** \brief Factorise a general matrix into its SV
52  decomposition using the Golub-Reinsch algorithm
53 
54  This factors matrix \c A of size <tt>(M,N)</tt> into
55  \f[
56  A = U~D~V^T
57  \f]
58  where \c U is a column-orthogonal matrix of size <tt>(M,N)</tt>
59  (stored in \c A), \c D is a diagonal matrix of size
60  <tt>(N,N)</tt> (stored in the vector \c S of size \c N), and \c
61  V is a orthogonal matrix of size <tt>(N,N)</tt>. The vector \c
62  work is a workspace vector of size \c N. The matrices \c U and
63  \c V are constructed so that
64  \f[
65  U^T~U = I \qquad \mathrm{and} \qquad V^T~V = V~V^T = I
66  \f]
67 
68  This algorithm requres \f$ M \geq N \f$.
69 
70  \todo Test N=1 case, N=2 case, and non-square matrices.
71  */
72  template<class mat_t, class mat2_t, class vec_t, class vec2_t>
73  void SV_decomp(size_t M, size_t N,
74  mat_t &A, mat2_t &V, vec_t &S, vec2_t &work) {
75 
76  if (M<N) {
77  O2SCL_ERR2("SVD decomposition of MxN matrix with M<N, ",
78  "is not implemented in SV_decomp().",o2scl::exc_eunimpl);
79  }
80 
81  // K should be the smaller of M and N, which is always N
82  const size_t K=N;
83 
84  // Handle the case of N = 1 (SVD of a column vector)
85 
86  if (N==1) {
87  double norm=O2SCL_CBLAS_NAMESPACE::dnrm2_subcol(A,0,0,M);
88  O2SCL_IX(S,0)=norm;
89  O2SCL_IX2(V,0,0)=1.0;
90 
91  if (norm != 0.0) {
92  O2SCL_CBLAS_NAMESPACE::dscal_subcol(A,0,0,M,1.0/norm);
93  }
94 
95  return;
96  }
97 
98  // [GSL] bidiagonalize matrix A, unpack A into U S V */
99  bidiag_decomp(M,N,A,S,work);
100 
101  //std::cout << "A: " << A(0,0) << " " << A(M-1,N-1) << std::endl;
102  //std::cout << "S: " << S[0] << " " << S[S.size()-1] << std::endl;
103 
104  bidiag_unpack2(M,N,A,S,work,V);
105 
106  //std::cout << "S2: " << S[0] << " " << S[S.size()-1] << std::endl;
107 
108  // [GSL] apply reduction steps to B=(S,Sd)
109 
110  chop_small_elements(N,S,work);
111 
112  //std::cout << "S3: " << S[0] << " " << S[S.size()-1] << std::endl;
113 
114  // [GSL] Progressively reduce the matrix until it is diagonal
115 
116  size_t b=N-1;
117  size_t iter=0;
118 
119  while (b>0) {
120 
121  double fbm1=O2SCL_IX(work,b-1);
122  if (fbm1==0.0 || gsl_isnan(fbm1)) {
123  b--;
124  continue;
125  }
126  //std::cout << "b,fbm1: " << b << " " << fbm1 << std::endl;
127 
128  // [GSL] Find the largest unreduced block (a,b) starting from b
129  // and working backwards
130 
131  size_t a=b-1;
132 
133  while (a>0) {
134 
135  double fam1=O2SCL_IX(work,a-1);
136  if (fam1==0.0 || gsl_isnan(fam1)) {
137  break;
138  }
139  a--;
140 
141  //std::cout << "a,fam1: " << a << " " << fam1 << std::endl;
142  }
143 
144  iter++;
145  if (iter>100*N) {
146  O2SCL_ERR("SV decomposition failed to converge in SV_decomp().",
148  }
149 
150  int rescale=0;
151  double scale=1.0;
152  double norm=0.0;
153  size_t n_block=b-a+1;
154 
155  // [GSL] Find the maximum absolute values of the diagonal
156  // and subdiagonal
157 
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;
161  //std::cout << "aa: " << aa << std::endl;
162  }
163 
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;
167  //std::cout << "aa2: " << aa << std::endl;
168  }
169 
170  // [GSL] Temporarily scale the submatrix if necessary
171 
172  if (norm>GSL_SQRT_DBL_MAX) {
173  scale=norm/GSL_SQRT_DBL_MAX;
174  rescale=1;
175  } else if (norm<GSL_SQRT_DBL_MIN && norm>0.0) {
176  scale=norm/GSL_SQRT_DBL_MIN;
177  rescale=1;
178  }
179 
180  //std::cout << "rescale: " << rescale << std::endl;
181 
182  if (rescale) {
183  O2SCL_CBLAS_NAMESPACE::dscal_subvec(1.0/scale,N,S,a);
184  O2SCL_CBLAS_NAMESPACE::dscal_subvec(1.0/scale,N,work,a);
185  }
186 
187  // [GSL] Perform the implicit QR step
188 
189  /*
190  for(size_t ii=0;ii<M;ii++) {
191  for(size_t jj=0;jj<N;jj++) {
192  std::cout << ii << "." << jj << "." << A(ii,jj) << std::endl;
193  }
194  }
195  for(size_t ii=0;ii<N;ii++) {
196  for(size_t jj=0;jj<N;jj++) {
197  std::cout << "V: " << ii << "." << jj << "." << V(ii,jj) << std::endl;
198  }
199  }
200  */
201 
202  // At this point, we want to perform qrstep on the arguments
203  // assuming S is a vector of size n_block,
204  // f is a vector of size n_block, U is matrix of size (M,n_block),
205  // and V is a matrix of size (N,n_block).
206  qrstep_sub(M,N,n_block+a,S,work,A,V,a);
207 
208  /*
209  for(size_t ii=0;ii<M;ii++) {
210  for(size_t jj=0;jj<N;jj++) {
211  std::cout << ii << " " << jj << " " << A(ii,jj) << std::endl;
212  }
213  }
214  for(size_t ii=0;ii<N;ii++) {
215  for(size_t jj=0;jj<N;jj++) {
216  std::cout << "V: " << ii << " " << jj << " " << V(ii,jj) << std::endl;
217  }
218  }
219  */
220 
221  // [GSL] Remove any small off-diagonal elements
222 
223  chop_small_elements(N,S,work);
224 
225  // [GSL] Undo the scaling if needed
226  if (rescale) {
228  O2SCL_CBLAS_NAMESPACE::dscal_subvec(scale,N,work,a);
229  }
230  }
231 
232  // [GSL] Make singular values positive by reflections if necessary
233 
234  for (size_t j=0;j<K;j++) {
235  double Sj=O2SCL_IX(S,j);
236  if (Sj<0.0) {
237  for (size_t i=0;i<N;i++) {
238  O2SCL_IX2(V,i,j)=-O2SCL_IX2(V,i,j);
239  }
240  O2SCL_IX(S,j)=-Sj;
241  }
242  }
243  //std::cout << "Here8" << std::endl;
244 
245  // [GSL] Sort singular values into decreasing order
246 
247  for (size_t i=0;i<K;i++) {
248 
249  double S_max=O2SCL_IX(S,i);
250  size_t i_max=i;
251 
252  for (size_t j=i+1;j<K;j++) {
253  double Sj=O2SCL_IX(S,j);
254  if (Sj>S_max) {
255  S_max=Sj;
256  i_max=j;
257  }
258  }
259 
260  if (i_max != i) {
261  // [GSL] swap eigenvalues
262  o2scl::vector_swap<vec_t,double>(S,i,i_max);
263 
264  // [GSL] swap eigenvectors
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);
267  }
268  }
269  //std::cout << "Here9" << std::endl;
270 
271  return;
272  }
273 
274  /** \brief SV decomposition by the modified Golub-Reinsch
275  algorithm which is better for \f$ M \gg N \f$
276 
277  This factors matrix \c A of size <tt>(M,N)</tt> into
278  \f[
279  A = U~D~V^T
280  \f]
281  where \c U is a column-orthogonal matrix of size <tt>(M,N)</tt>
282  (stored in \c A), \c D is a diagonal matrix of size
283  <tt>(N,N)</tt> (stored in the vector \c S of size \c N), and \c
284  V is a orthogonal matrix of size <tt>(N,N)</tt>. The vector \c
285  work is a workspace vector of size \c N and the matrix
286  \c X is a workspace of size <tt>(N,N)</tt>.
287  */
288  template<class mat_t, class mat2_t, class mat3_t, class vec_t,
289  class vec2_t> void SV_decomp_mod
290  (size_t M, size_t N, mat_t &A, mat2_t &X, mat3_t &V, vec_t &S,
291  vec2_t &work) {
292 
293  if (N == 1) {
294  double norm=O2SCL_CBLAS_NAMESPACE::dnrm2_subcol(A,0,0,M);
295  O2SCL_IX(S,0)=norm;
296  O2SCL_IX2(V,0,0)=1.0;
297 
298  if (norm != 0.0) {
299  O2SCL_CBLAS_NAMESPACE::dscal_subcol(A,0,0,M,1.0/norm);
300  }
301 
302  return;
303  }
304 
305  // AWS: This next loop is just a QR decomposition. Replace
306  // with QR_decomp()?
307 
308  // [GSL] Convert A into an upper triangular matrix R
309  for (size_t i=0;i<N;i++) {
310  double tau_i=householder_transform_subcol(A,i,i,M);
311  // [GSL] Apply the transformation to the remaining columns
312  if (i+1<N) {
313  householder_hm_subcol(A,i,i+1,M,N,A,i,i,tau_i);
314  }
315  O2SCL_IX(S,i)=tau_i;
316  }
317 
318  // [GSL] Copy the upper triangular part of A into X
319 
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;
323  }
324 
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);
328  }
329  }
330 
331  // [GSL] Convert A into an orthogonal matrix L
332 
333  for (size_t j=N;j-- > 0;) {
334  // [GSL] Householder column transformation to accumulate L
335  double tj=O2SCL_IX(S,j);
336  householder_hm1_sub(M,N,tj,A,j,j);
337  }
338 
339  // [GSL] unpack R into X V S
340 
341  SV_decomp(N,N,X,V,S,work);
342 
343  // [GSL] Multiply L by X, to obtain U = L X, stored in U
344 
345  for (size_t i=0;i<M;i++) {
346 
347  for(size_t j=0;j<N;j++) {
348  O2SCL_IX(work,j)=0.0;
349  }
350 
351  for (size_t j=0;j<N;j++) {
352  double Lij=O2SCL_IX2(A,i,j);
353  O2SCL_CBLAS_NAMESPACE::daxpy_subrow(Lij,N,X,j,0,work);
354  }
355 
356  for (size_t j=0;j<N;j++) {
357  O2SCL_IX2(A,i,j)=O2SCL_IX(work,j);
358  }
359 
360  }
361 
362  return;
363  }
364 
365  /** \brief Solve the system A x = b using the SV decomposition
366 
367  Solves a linear system using the output of \ref SV_decomp().
368  Only non-zero singular values are used in computing solution. In
369  the over-determined case, \f$ M>N \f$, the system is solved in
370  the least-squares sense.
371  */
372  template<class mat_t, class mat2_t, class vec_t, class vec2_t, class vec3_t>
373  void SV_solve(size_t M, size_t N,
374  mat_t &U, mat2_t &V, vec_t &S, vec2_t &b, vec3_t &x) {
375 
376  //if (U->size1 != b->size)
377  //else if (U->size2 != S->size)
378  //else if (V->size1 != V->size2)
379  //else if (S->size != V->size1)
380  //else if (V->size2 != x->size)
381 
382  double *w=new double[N];
383 
385  (O2SCL_CBLAS_NAMESPACE::o2cblas_RowMajor,
386  O2SCL_CBLAS_NAMESPACE::o2cblas_Trans,M,N,1.0,U,b,0.0,w);
387 
388  for (size_t i=0;i<N;i++) {
389  double alpha=O2SCL_IX(S,i);
390  if (alpha != 0) alpha=1.0/alpha;
391  w[i]*=alpha;
392  }
393 
395  (O2SCL_CBLAS_NAMESPACE::o2cblas_RowMajor,
396  O2SCL_CBLAS_NAMESPACE::o2cblas_NoTrans,N,N,1.0,V,w,0.0,x);
397 
398  delete[] w;
399 
400  return;
401  }
402 
403  /** \brief SV decomposition using one-sided Jacobi orthogonalization
404 
405  This factors matrix \c A of size <tt>(M,N)</tt> into
406  \f[
407  A = U~D~V^T
408  \f]
409  where \c U is a column-orthogonal matrix of size <tt>(M,N)</tt>
410  (stored in \c A), \c D is a diagonal matrix of size
411  <tt>(N,N)</tt> (stored in the vector \c S of size \c N), and \c
412  V is a orthogonal matrix of size <tt>(N,N)</tt>.
413 
414  This function computes singular values to higher relative
415  accuracy than \ref SV_decomp() and \ref SV_decomp_mod().
416 
417  \comment
418  Algorithm due to J.C. Nash, Compact Numerical Methods for
419  Computers (New York: Wiley and Sons, 1979), chapter 3.
420  See also Algorithm 4.1 in
421  James Demmel, Kresimir Veselic, "Jacobi's Method is more
422  accurate than QR", Lapack Working Note 15 (LAWN15), October 1989.
423  Available from netlib.
424 
425  Based on code by Arthur Kosowsky, Rutgers University
426  kosowsky@physics.rutgers.edu
427 
428  Another relevant paper is, P.P.M. De Rijk, "A One-Sided Jacobi
429  Algorithm for computing the singular value decomposition on a
430  vector computer", SIAM Journal of Scientific and Statistical
431  Computing, Vol 10, No 2, pp 359-371, March 1989.
432  \endcomment
433 
434  \future There were originally a few GSL_COERCE_DBL calls which
435  have been temporarily removed and could be restored.
436  */
437  template<class mat_t, class mat2_t, class vec_t>
438  void SV_decomp_jacobi(size_t M, size_t N, mat_t &A, mat2_t &Q, vec_t &S) {
439 
440  //if (A->size1<A->size2)
441  //else if (Q->size1 != A->size2)
442  //else if (Q->size1 != Q->size2)
443  //else if (S->size != A->size2)
444  //const size_t M=A->size1;
445  //const size_t N=A->size2;
446 
447  // [GSL] Initialize the rotation counter and the sweep counter.
448  int count=1;
449  int sweep=0;
450  int sweepmax=5*N;
451 
452  double dbl_eps=std::numeric_limits<double>::epsilon();
453 
454  double tolerance=10*M*dbl_eps;
455 
456  // [GSL] Always do at least 12 sweeps.
457  if (sweepmax<12) sweepmax=12;
458 
459  // [GSL] Set Q to the identity matrix.
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;
464  }
465  }
466 
467  // [GSL] Store the column error estimates in S, pfor use during the
468  // orthogonalization
469 
470  for (size_t j=0;j<N;j++) {
471  double sj=O2SCL_CBLAS_NAMESPACE::dnrm2_subcol(A,0,j,M);
472  O2SCL_IX(S,j)=dbl_eps*sj;
473  }
474 
475  // [GSL] Orthogonalize A by plane rotations.
476 
477  while (count > 0 && sweep <= sweepmax) {
478 
479  // [GSL] Initialize rotation counter.
480  count=N*(N-1)/2;
481 
482  for (size_t j=0;j<N-1;j++) {
483  for (size_t k=j+1;k<N;k++) {
484 
485  double a=0.0;
486  double b=0.0;
487  double p=0.0;
488  double q=0.0;
489  double cosine, sine;
490  double v;
491  double abserr_a, abserr_b;
492  int sorted, orthog, noisya, noisyb;
493 
494  //gsl_blas_ddot(&cj.vector,&ck.vector,&p);
495  for(size_t ii=0;ii<M;ii++) {
496  p+=O2SCL_IX2(A,ii,j)*O2SCL_IX2(A,ii,k);
497  }
498  // [GSL] equation 9a: p = 2 x.y
499  p*=2.0;
500 
503 
504  q=a*a-b*b;
505  v=hypot(p,q);
506 
507  // [GSL] test for columns j,k orthogonal, or dominant errors
508 
509  abserr_a=O2SCL_IX(S,j);
510  abserr_b=O2SCL_IX(S,k);
511 
512  //sorted=(GSL_COERCE_DBL(a) >= GSL_COERCE_DBL(b));
513  //orthog=(fabs (p) <= tolerance*GSL_COERCE_DBL(a*b));
514  sorted=(a>=b);
515  orthog=(fabs (p) <= tolerance*a*b);
516  noisya=(a<abserr_a);
517  noisyb=(b<abserr_b);
518 
519  if (sorted && (orthog || noisya || noisyb)) {
520  count--;
521  continue;
522  }
523 
524  // [GSL] calculate rotation angles
525  if (v == 0 || !sorted) {
526  cosine=0.0;
527  sine=1.0;
528  } else {
529  cosine=sqrt((v+q)/(2.0*v));
530  sine=p/(2.0*v*cosine);
531  }
532 
533  // [GSL] apply rotation to A
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;
539  }
540 
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;
543 
544  // [GSL] apply rotation to Q
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;
550  }
551  }
552  }
553 
554  // [GSL] Sweep completed.
555  sweep++;
556  }
557 
558  // [GSL] Orthogonalization complete. Compute singular values.
559 
560  {
561  double prev_norm=-1.0;
562 
563  for (size_t j=0;j<N;j++) {
564 
565  double norm=O2SCL_CBLAS_NAMESPACE::dnrm2_subcol(A,0,j,M);
566 
567  // [GSL] Determine if singular value is zero, according to the
568  // criteria used in the main loop above (i.e. comparison
569  // with norm of previous column).
570 
571  if (norm == 0.0 || prev_norm == 0.0
572  || (j > 0 && norm <= tolerance*prev_norm)) {
573 
574  // [GSL] singular
575  O2SCL_IX(S,j)=0.0;
576  // [GSL] annihilate column
577  for(size_t ii=0;ii<M;ii++) {
578  O2SCL_IX2(A,ii,j)=0.0;
579  }
580 
581  prev_norm=0.0;
582 
583  } else {
584 
585  // [GSL] non-singular
586  O2SCL_IX(S,j)=norm;
587  // [GSL] normalize column
588  O2SCL_CBLAS_NAMESPACE::dscal_subcol(A,0,j,M,1.0/norm);
589 
590  prev_norm=norm;
591 
592  }
593  }
594  }
595 
596  if (count > 0) {
597  // [GSL] reached sweep limit
598  O2SCL_ERR("Jacobi iterations did not reach desired tolerance",
600  }
601 
602  return;
603  }
604 
605  /** \brief Balance a general matrix A by scaling the columns
606  by the diagonal matrix D
607 
608  The function computes
609  \f$
610  A_{\mathrm{new}} = A D^{-1}
611  \f$
612  where \f$ D \f$ is a diagonal matrix.
613  */
614  template<class mat_t, class vec_t>
615  void balance_columns(size_t M, size_t N, mat_t &A, vec_t &D) {
616 
617  for(size_t j=0;j<N;j++) O2SCL_IX(D,j)=1.0;
618 
619  for(size_t j=0;j<N;j++) {
620  double s=O2SCL_CBLAS_NAMESPACE::dasum_subcol(A,0,j,M);
621  double f=1.0;
622  if (s==0.0 || !std::isfinite(s)) {
623  O2SCL_IX(D,j)=f;
624  continue;
625  }
626  while (s>1.0) {
627  s/=2.0;
628  f*=2.0;
629  }
630  while (s<0.5) {
631  s*=2.0;
632  f/=2.0;
633  }
634  O2SCL_IX(D,j)=f;
635  if (f!=1.0) {
637  }
638  }
639 
640  return;
641  }
642 
643 #ifdef DOXYGEN
644 }
645 #endif
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.
Definition: bidiag_base.h:170
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.
Definition: svd_base.h:373
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.
Definition: svd_base.h:615
exceeded max number of iterations
Definition: err_hnd.h:73
failed to reach the specified tolerance
Definition: err_hnd.h:79
requested feature not (yet) implemented
Definition: err_hnd.h:99
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 .
Definition: cblas_base.h:218
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.
Definition: err_hnd.h:281
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.
Definition: cblas_base.h:1131
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
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.
Definition: cblas_base.h:1044
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
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()
Definition: svdstep_base.h:762
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.
Definition: svdstep_base.h:58
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.
Definition: svd_base.h:73
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 .
Definition: svd_base.h:290
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.
Definition: bidiag_base.h:65
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.
Definition: svd_base.h:438
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.
Definition: cblas_base.h:1177
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.
Definition: cblas_base.h:920
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.
Definition: cblas_base.h:1093
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.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).