Mercurial > jwe > octave

changeset 22436:09005ac7d56c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gsvd<T>: add class template support for FloatMatrix and FloatComplexMatrix. * liboctave/numeric/gsvd.h: template class more so that it works with FloatMatrix and ComplexFloatMatrix. * liboctave/numeric/gsvd.cc: add methods FloatMatrix and FloatComplexMatrix. * liboctave/numeric/lo-lapack-proto.h: declare interface to LAPACK's CGGSVD and SGGSVD.
author Carnë Draug <carandraug@octave.org>
date Sun, 04 Sep 2016 17:12:03 +0100
parents de24ca103c21
children 0aee8b620864
files liboctave/numeric/gsvd.cc liboctave/numeric/gsvd.h liboctave/numeric/lo-lapack-proto.h
diffstat 3 files changed, 141 insertions(+), 21 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/numeric/gsvd.cc	Sun Sep 04 16:09:02 2016 +0100
+++ b/liboctave/numeric/gsvd.cc	Sun Sep 04 17:12:03 2016 +0100
@@ -19,14 +19,19 @@
 #  include <config.h>
 #endif
 
+#include <vector>
+
 #include "gsvd.h"
+
 #include "lo-error.h"
 #include "lo-lapack-proto.h"
+#include "dMatrix.h"
+#include "fMatrix.h"
 #include "CMatrix.h"
+#include "fCMatrix.h"
 #include "dDiagMatrix.h"
-#include "dMatrix.h"
+#include "fDiagMatrix.h"
 
-#include <vector>
 
 template <>
 void
@@ -36,7 +41,8 @@
                      double *tmp_dataB, octave_idx_type p1, Matrix& alpha,
                      Matrix& beta, double *u, octave_idx_type nrow_u, double *v,
                      octave_idx_type nrow_v, double *q, octave_idx_type nrow_q,
-                     Matrix& work, octave_idx_type* iwork, octave_idx_type& info)
+                     Matrix& work, octave_idx_type* iwork,
+                     octave_idx_type& info)
 {
   F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
                              F77_CONST_CHAR_ARG2 (&jobv, 1),
@@ -51,6 +57,31 @@
                              F77_CHAR_ARG_LEN (1)));
 }
 
+template <>
+void
+gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, octave_idx_type m,
+                          octave_idx_type n, octave_idx_type p,
+                          octave_idx_type& k, octave_idx_type& l,
+                          float *tmp_dataA, octave_idx_type m1,
+                          float *tmp_dataB, octave_idx_type p1,
+                          FloatMatrix& alpha, FloatMatrix& beta, float *u,
+                          octave_idx_type nrow_u, float *v,
+                          octave_idx_type nrow_v, float *q,
+                          octave_idx_type nrow_q, FloatMatrix& work,
+                          octave_idx_type* iwork, octave_idx_type& info)
+{
+  F77_XFCN (sggsvd, SGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
+                             F77_CONST_CHAR_ARG2 (&jobv, 1),
+                             F77_CONST_CHAR_ARG2 (&jobq, 1),
+                             m, n, p, k, l, tmp_dataA, m1,
+                             tmp_dataB, p1, alpha.fortran_vec (),
+                             beta.fortran_vec (), u, nrow_u,
+                             v, nrow_v, q, nrow_q, work.fortran_vec (),
+                             iwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
+}
 
 template <>
 void
@@ -82,6 +113,38 @@
                              F77_CHAR_ARG_LEN (1)));
 }
 
+template <>
+void
+gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
+                                 octave_idx_type m, octave_idx_type n,
+                                 octave_idx_type p, octave_idx_type& k,
+                                 octave_idx_type& l, FloatComplex *tmp_dataA,
+                                 octave_idx_type m1, FloatComplex *tmp_dataB,
+                                 octave_idx_type p1, FloatMatrix& alpha,
+                                 FloatMatrix& beta, FloatComplex *u,
+                                 octave_idx_type nrow_u, FloatComplex *v,
+                                 octave_idx_type nrow_v, FloatComplex *q,
+                                 octave_idx_type nrow_q,
+                                 FloatComplexMatrix& work,
+                                 octave_idx_type* iwork, octave_idx_type& info)
+{
+  FloatMatrix rwork(2*n, 1);
+  F77_XFCN (cggsvd, CGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
+                             F77_CONST_CHAR_ARG2 (&jobv, 1),
+                             F77_CONST_CHAR_ARG2 (&jobq, 1),
+                             m, n, p, k, l, F77_CMPLX_ARG (tmp_dataA), m1,
+                             F77_CMPLX_ARG (tmp_dataB), p1,
+                             alpha.fortran_vec (), beta.fortran_vec (),
+                             F77_CMPLX_ARG (u), nrow_u,
+                             F77_CMPLX_ARG (v), nrow_v,
+                             F77_CMPLX_ARG (q), nrow_q,
+                             F77_CMPLX_ARG (work.fortran_vec ()),
+                             rwork.fortran_vec (), iwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
+}
+
 template <typename T>
 T
 gsvd<T>::left_singular_matrix_A (void) const
@@ -207,8 +270,8 @@
   lwork = (lwork > p ? lwork : p) + n;
 
   T work (lwork, 1);
-  Matrix alpha (n, 1);
-  Matrix beta (n, 1);
+  real_matrix alpha (n, 1);
+  real_matrix beta (n, 1);
 
   std::vector<octave_idx_type> iwork (n);
 
@@ -288,8 +351,7 @@
 }
 
 // Instantiations we need.
-
 template class gsvd<Matrix>;
-
+template class gsvd<FloatMatrix>;
 template class gsvd<ComplexMatrix>;
-
+template class gsvd<FloatComplexMatrix>;
--- a/liboctave/numeric/gsvd.h	Sun Sep 04 16:09:02 2016 +0100
+++ b/liboctave/numeric/gsvd.h	Sun Sep 04 17:12:03 2016 +0100
@@ -20,9 +20,6 @@
 
 #include "octave-config.h"
 
-#include "dDiagMatrix.h"
-#include "dMatrix.h"
-
 template <typename T>
 class
 gsvd
@@ -64,8 +61,11 @@
 
   ~gsvd (void) { }
 
-  DiagMatrix singular_values_A (void) const { return sigmaA; }
-  DiagMatrix singular_values_B (void) const { return sigmaB; }
+  typename T::real_diag_matrix_type
+  singular_values_A (void) const { return sigmaA; }
+
+  typename T::real_diag_matrix_type
+  singular_values_B (void) const { return sigmaB; }
 
   T left_singular_matrix_A (void) const;
   T left_singular_matrix_B (void) const;
@@ -74,22 +74,21 @@
   T R_matrix (void) const;
 
 private:
+  typedef typename T::value_type P;
+  typedef typename T::real_matrix_type real_matrix;
 
   gsvd::Type type;
-
-  typedef typename T::element_type P;
-
-  DiagMatrix sigmaA, sigmaB;
+  typename T::real_diag_matrix_type sigmaA, sigmaB;
   T left_smA, left_smB;
   T right_sm, R;
 
   void ggsvd (char& jobu, char& jobv, char& jobq, octave_idx_type m,
               octave_idx_type n, octave_idx_type p, octave_idx_type& k,
               octave_idx_type& l, P *tmp_dataA, octave_idx_type m1,
-              P *tmp_dataB, octave_idx_type p1, Matrix& alpha, Matrix& beta,
-              P *u, octave_idx_type nrow_u, P *v, octave_idx_type nrow_v, P *q,
-              octave_idx_type nrow_q, T& work, octave_idx_type* iwork,
-              octave_idx_type& info);
+              P *tmp_dataB, octave_idx_type p1, real_matrix& alpha,
+              real_matrix& beta, P *u, octave_idx_type nrow_u, P *v,
+              octave_idx_type nrow_v, P *q, octave_idx_type nrow_q, T& work,
+              octave_idx_type* iwork, octave_idx_type& info);
 };
 
 #endif
--- a/liboctave/numeric/lo-lapack-proto.h	Sun Sep 04 16:09:02 2016 +0100
+++ b/liboctave/numeric/lo-lapack-proto.h	Sun Sep 04 17:12:03 2016 +0100
@@ -839,6 +839,35 @@
      F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
+  F77_FUNC (sggsvd, SGGSVD)
+   (F77_CONST_CHAR_ARG_DECL,    // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     F77_REAL*,                 // A
+     const F77_INT&,            // LDA
+     F77_REAL*,                 // B
+     const F77_INT&,            // LDB
+     F77_REAL*,                 // ALPHA
+     F77_REAL*,                 // BETA
+     F77_REAL*,                 // U
+     const F77_INT&,            // LDU
+     F77_REAL*,                 // V
+     const F77_INT&,            // LDV
+     F77_REAL*,                 // Q
+     const F77_INT&,            // LDQ
+     F77_REAL*,                 // WORK
+     int*,                      // IWORK
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
   F77_FUNC (zggsvd, ZGGSVD)
     (F77_CONST_CHAR_ARG_DECL,   // JOBU
      F77_CONST_CHAR_ARG_DECL,   // JOBV
@@ -868,6 +897,36 @@
      F77_CHAR_ARG_LEN_DECL
      F77_CHAR_ARG_LEN_DECL);
 
+  F77_RET_T
+  F77_FUNC (cggsvd, CGGSVD)
+   (F77_CONST_CHAR_ARG_DECL,    // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     F77_CMPLX*,                // A
+     const F77_INT&,            // LDA
+     F77_CMPLX*,                // B
+     const F77_INT&,            // LDB
+     F77_REAL*,                 // ALPHA
+     F77_REAL*,                 // BETA
+     F77_CMPLX*,                // U
+     const F77_INT&,            // LDU
+     F77_CMPLX*,                // V
+     const F77_INT&,            // LDV
+     F77_CMPLX*,                // Q
+     const F77_INT&,            // LDQ
+     F77_CMPLX*,                // WORK
+     F77_REAL*,                 // RWORK
+     int*,                      // IWORK
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+
   // GTSV
 
   F77_RET_T