Mercurial > octave

changeset 29418:9245ae55b6bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use Octave coding conventions in patch for GEJSV SVD driver (bug #55727). * libinterp/corefcn/svd.cc, liboctave/numeric/svd.cc: Eliminate trailing spaces. Use indent of 2 spaces in #define. Use space between function name and opening parenthesis.
author Rik <rik@octave.org>
date Tue, 09 Mar 2021 14:23:47 -0800
parents a6ab7069a87c
children 2cfdcae08e84
files libinterp/corefcn/svd.cc liboctave/numeric/svd.cc
diffstat 2 files changed, 46 insertions(+), 48 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/svd.cc	Tue Mar 09 21:05:54 2021 +0100
+++ b/libinterp/corefcn/svd.cc	Tue Mar 09 14:23:47 2021 -0800
@@ -167,8 +167,8 @@
 singular matrices in addition to singular values) there is a choice of two
 routines in @sc{lapack}.  The default routine used by Octave is @code{gesvd}.
 The alternative is @code{gesdd} which is 5X faster, but may use more memory
-and may be inaccurate for some input matrices.  There is a third routine 
-@code{gejsv}, suitable for better accuracy at extreme scale.  See the 
+and may be inaccurate for some input matrices.  There is a third routine
+@code{gejsv}, suitable for better accuracy at extreme scale.  See the
 documentation for @code{svd_driver} for more information on choosing a driver.
 @seealso{svd_driver, svds, eig, lu, chol, hess, qr, qz}
 @end deftypefn */)
@@ -421,8 +421,8 @@
 
 The routine @code{gejsv} uses a preconditioned Jacobi SVD algorithm.  Unlike
 @code{gesvd} and @code{gesdd}, in @code{gejsv}, there is no bidiagonalization
-step that could contaminate accuracy in some extreme case. Also, @code{gejsv}
-is shown to be optimally accurate in some sense.  However, the speed is slower
+step that could contaminate accuracy in some extreme cases.  Also, @code{gejsv}
+is known to be optimally accurate in some sense.  However, the speed is slower
 (single threaded at its core) and uses more memory (O(min(M,N) ^ 2 + M + N)).
 
 Beyond speed and memory issues, there have been instances where some input
--- a/liboctave/numeric/svd.cc	Tue Mar 09 21:05:54 2021 +0100
+++ b/liboctave/numeric/svd.cc	Tue Mar 09 14:23:47 2021 -0800
@@ -62,7 +62,7 @@
 
 private:
   typedef typename T::element_type P;
-  
+
   // functions could be called from GEJSV
   static F77_INT geqp3_lwork (F77_INT m, F77_INT n,
                               P *a, F77_INT lda,
@@ -92,21 +92,19 @@
                               P *work, F77_INT lwork, F77_INT& info);
 };
 
-#define GEJSV_REAL_QP3_LWORK(f, F)                          \
-    F77_XFCN (f, F, (m, n, a, lda, jpvt,                    \
-                     tau, work, lwork, info))
+#define GEJSV_REAL_QP3_LWORK(f, F)                              \
+  F77_XFCN (f, F, (m, n, a, lda, jpvt, tau, work, lwork, info))
 
-#define GEJSV_REAL_QR_LWORK(f, F)                           \
-    F77_XFCN (f, F, (m, n, a, lda,                          \
-                     tau, work, lwork, info))
+#define GEJSV_REAL_QR_LWORK(f, F)                               \
+  F77_XFCN (f, F, (m, n, a, lda, tau, work, lwork, info))
 
-#define GEJSV_REAL_ORM_LWORK(f, F)                          \
-    F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&side,  1),       \
-                     F77_CONST_CHAR_ARG2 (&trans, 1),       \
-                     m, n, k, a, lda, tau,                  \
-                     c, ldc, work, lwork, info              \
-                     F77_CHAR_ARG_LEN (1)                   \
-                     F77_CHAR_ARG_LEN (1)))
+#define GEJSV_REAL_ORM_LWORK(f, F)                              \
+  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&side,  1),             \
+                   F77_CONST_CHAR_ARG2 (&trans, 1),             \
+                   m, n, k, a, lda, tau,                        \
+                   c, ldc, work, lwork, info                    \
+                   F77_CHAR_ARG_LEN (1)                         \
+                   F77_CHAR_ARG_LEN (1)))
 
 // For Matrix
 template<>
@@ -235,7 +233,7 @@
 {
   F77_INT lwork = -1;
   std::vector<P> work (2);  // dummy work space
-  
+
   // variables that mimic running environment of gejsv
   F77_INT lda  = std::max<F77_INT> (m, 1);
   F77_INT ierr = 0;
@@ -331,23 +329,23 @@
     // GESVD specializations
 
 #define GESVD_REAL_STEP(f, F)                                   \
-    F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1),            \
-                     F77_CONST_CHAR_ARG2 (&jobv, 1),            \
-                     m, n, tmp_data, m1, s_vec, u, m1, vt,      \
-                     nrow_vt1, work.data (), lwork, info        \
-                     F77_CHAR_ARG_LEN (1)                       \
-                     F77_CHAR_ARG_LEN (1)))
+  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1),              \
+                   F77_CONST_CHAR_ARG2 (&jobv, 1),              \
+                   m, n, tmp_data, m1, s_vec, u, m1, vt,        \
+                   nrow_vt1, work.data (), lwork, info          \
+                   F77_CHAR_ARG_LEN (1)                         \
+                   F77_CHAR_ARG_LEN (1)))
 
 #define GESVD_COMPLEX_STEP(f, F, CMPLX_ARG)             \
-    F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1),    \
-                     F77_CONST_CHAR_ARG2 (&jobv, 1),    \
-                     m, n, CMPLX_ARG (tmp_data),        \
-                     m1, s_vec, CMPLX_ARG (u), m1,      \
-                     CMPLX_ARG (vt), nrow_vt1,          \
-                     CMPLX_ARG (work.data ()),          \
-                     lwork, rwork.data (), info         \
-                     F77_CHAR_ARG_LEN (1)               \
-                     F77_CHAR_ARG_LEN (1)))
+  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1),      \
+                   F77_CONST_CHAR_ARG2 (&jobv, 1),      \
+                   m, n, CMPLX_ARG (tmp_data),          \
+                   m1, s_vec, CMPLX_ARG (u), m1,        \
+                   CMPLX_ARG (vt), nrow_vt1,            \
+                   CMPLX_ARG (work.data ()),            \
+                   lwork, rwork.data (), info           \
+                   F77_CHAR_ARG_LEN (1)                 \
+                   F77_CHAR_ARG_LEN (1)))
 
     // DGESVD
     template<>
@@ -583,9 +581,9 @@
                         F77_INT& lwork, std::vector<F77_INT>& iwork,
                         F77_INT& info)
     {
-      lwork = gejsv_lwork<Matrix>::optimal(joba, jobu, jobv, m, n);
+      lwork = gejsv_lwork<Matrix>::optimal (joba, jobu, jobv, m, n);
       work.reserve (lwork);
-      
+
       GEJSV_REAL_STEP (dgejsv, DGEJSV);
     }
 
@@ -600,7 +598,7 @@
                              F77_INT& lwork, std::vector<F77_INT>& iwork,
                              F77_INT& info)
     {
-      lwork = gejsv_lwork<FloatMatrix>::optimal(joba, jobu, jobv, m, n);
+      lwork = gejsv_lwork<FloatMatrix>::optimal (joba, jobu, jobv, m, n);
       work.reserve (lwork);
 
       GEJSV_REAL_STEP (sgejsv, SGEJSV);
@@ -622,13 +620,13 @@
       work.reserve (2);
 
       GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
-      
+
       lwork = static_cast<F77_INT> (work[0].real ());
       work.reserve (lwork);
-      
+
       lrwork = static_cast<F77_INT> (rwork[0]);
       rwork.reserve (lrwork);
-      
+
       F77_INT liwork = static_cast<F77_INT> (iwork[0]);
       iwork.reserve (liwork);
 
@@ -651,13 +649,13 @@
       work.reserve (2);
 
       GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
-      
+
       lwork = static_cast<F77_INT> (work[0].real ());
       work.reserve (lwork);
-      
+
       lrwork = static_cast<F77_INT> (rwork[0]);
       rwork.reserve (lrwork);
-      
+
       F77_INT liwork = static_cast<F77_INT> (iwork[0]);
       iwork.reserve (liwork);
 
@@ -805,7 +803,7 @@
               u  = right_sm.fortran_vec ();
               vt = left_sm.fortran_vec ();
             }
-          
+
           // translate jobu and jobv from gesvd to gejsv.
           assert ('A' <= jobu && jobu <= 'S' && 'A' <= jobv && jobv <= 'S');
           char job_svd2jsv[1 + 'S' - 'A'][2] = {0};
@@ -819,22 +817,22 @@
           job_svd2jsv['N' - 'A'][1] = 'N';
           jobu = job_svd2jsv[jobu - 'A'][0];
           jobv = job_svd2jsv[jobv - 'A'][1];
-          
+
           char joba = 'F';  // 'F': most conservative
           char jobr = 'R';  // 'R' is recommended.
           char jobt = 'N';  // or 'T', but that requires U and V appear together
           char jobp = 'N';  // use 'P' if denormal is poorly implemented.
-          
+
           std::vector<F77_INT> iwork (std::max<F77_INT> (m + 3*n, 1));
 
           gejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, tmp_data, m1,
                  s_vec, u, vt, nrow_vt1, work, lwork, iwork, info);
-          
+
           if (iwork[2] == 1)
             (*current_liboctave_warning_with_id_handler)
             ("Octave:convergence", "svd: (driver: GEJSV) "
              "Denormal occured, possible loss of accuracy.");
-          
+
           if (info < 0)
             (*current_liboctave_error_handler)
             ("svd: (driver: GEJSV) Illegal argument at #%d",