Mercurial > octave-nkf

changeset 10822:23d2378512a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement rsf2csf
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 27 Jul 2010 11:26:43 +0200
parents 693e22af08ae
children 3d89d262f5d4
files libcruft/ChangeLog libcruft/lapack-xtra/crsf2csf.f libcruft/lapack-xtra/module.mk libcruft/lapack-xtra/zrsf2csf.f liboctave/ChangeLog liboctave/CmplxSCHUR.cc liboctave/CmplxSCHUR.h liboctave/dbleSCHUR.cc liboctave/dbleSCHUR.h liboctave/fCmplxSCHUR.cc liboctave/fCmplxSCHUR.h liboctave/floatSCHUR.cc liboctave/floatSCHUR.h src/ChangeLog src/DLD-FUNCTIONS/schur.cc
diffstat 15 files changed, 398 insertions(+), 1 deletions(-) [+]
line wrap: on
line diff
--- a/libcruft/ChangeLog	Mon Jul 26 21:25:36 2010 -0700
+++ b/libcruft/ChangeLog	Tue Jul 27 11:26:43 2010 +0200
@@ -1,3 +1,8 @@
+2010-07-27  Jaroslav Hajek  <highegg@gmail.com>
+
+	* lapack-xtra/zrsf2csf.f, lapack-xtra/crsf2csf.f: New sources.
+	* lapack-xtra/module.mk: Add them.
+
 2010-05-04  John W. Eaton  <jwe@octave.org>
 
 	* villad/dfopr.f, villad/dif.f, villad/intrp.f, villad/jcobi.f,
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/lapack-xtra/crsf2csf.f	Tue Jul 27 11:26:43 2010 +0200
@@ -0,0 +1,101 @@
+c Copyright (C) 2010  VZLU Prague, a.s., Czech Republic
+c
+c Author: Jaroslav Hajek <highegg@gmail.com>
+c
+c This file is part of Octave.
+c
+c Octave is free software; you can redistribute it and/or modify
+c it under the terms of the GNU General Public License as published by
+c the Free Software Foundation; either version 3 of the License, or
+c (at your option) any later version.
+c
+c This program is distributed in the hope that it will be useful,
+c but WITHOUT ANY WARRANTY; without even the implied warranty of
+c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+c GNU General Public License for more details.
+c
+c You should have received a copy of the GNU General Public License
+c along with this software; see the file COPYING.  If not, see
+c <http://www.gnu.org/licenses/>.
+c
+
+       subroutine crsf2csf(n,t,u,c,s)
+       integer n
+       complex t(n,n),u(n,n)
+       real c(n-1),s(n-1)
+       real x,y,z
+       integer j
+       j = 1
+       do while (j < n)
+c apply previous rotations to rows
+         call crcrot1(j,t(1,j),c,s)
+
+         y = t(j+1,j)
+         if (y /= 0) then
+c 2x2 block, form Givens rotation [c, i*s; i*s, c] 
+           x = t(j,j)
+           z = t(j,j+1)
+           c(j) = sqrt(abs(z/(y-z)))
+           s(j) = sign(sqrt(abs(y/(y-z))),z)
+c apply new rotation to t(j:j+1,j)
+           call crcrot1(2,t(j,j),c(j),s(j))
+c apply all rotations to t(1:j+1,j+1)
+           call crcrot1(j+1,t(1,j+1),c,s)
+c apply new rotation to columns j,j+1
+           call crcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))           
+c zero subdiagonal entry, skip next row
+           t(j+1,j) = 0
+           c(j+1) = 1
+           j = j + 2
+         else
+           c(j) = 1
+           j = j + 1
+         end if
+       end do
+
+c apply rotations to last column if needed
+       if (j == n) then
+         call crcrot1(j,t(1,j),c,s)
+       end if
+
+c apply stored rotations to all columns of u
+       do j = 1,n-1
+         if (c(j) /= 1) then
+           call crcrot2(n,u(1,j),u(1,j+1),c(j),s(j))
+         end if
+       end do
+
+       end subroutine
+
+       subroutine crcrot1(n,x,c,s)
+c apply rotations to a column from the left
+       integer n
+       complex x(n), t
+       real c(n-1),s(n-1)
+       integer i
+       do i = 1,n-1
+         if (c(i) /= 1) then
+           t = x(i)*c(i) - x(i+1)*cmplx(0,s(i))
+           x(i+1) = x(i+1)*c(i) - x(i)*cmplx(0,s(i))
+           x(i) = t
+         endif
+       end do
+       end subroutine
+
+       subroutine crcrot2(n,x,y,c,s)
+c apply a single rotation from the right to a pair of columns
+       integer n
+       complex x(n),y(n),t
+       real c, s
+       integer i
+       do i = 1,n
+         t = x(i)*c + y(i)*cmplx(0,s)
+         y(i) = y(i)*c + x(i)*cmplx(0,s)
+         x(i) = t
+       end do
+       end subroutine
+
+
+
+
+
--- a/libcruft/lapack-xtra/module.mk	Mon Jul 26 21:25:36 2010 -0700
+++ b/libcruft/lapack-xtra/module.mk	Tue Jul 27 11:26:43 2010 +0200
@@ -7,4 +7,6 @@
   lapack-xtra/xilaenv.f \
   lapack-xtra/xslamch.f \
   lapack-xtra/xslange.f \
-  lapack-xtra/xzlange.f
+  lapack-xtra/xzlange.f \
+  lapack-xtra/zrsf2csf.f \
+  lapack-xtra/crsf2csf.f
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/lapack-xtra/zrsf2csf.f	Tue Jul 27 11:26:43 2010 +0200
@@ -0,0 +1,101 @@
+c Copyright (C) 2010  VZLU Prague, a.s., Czech Republic
+c
+c Author: Jaroslav Hajek <highegg@gmail.com>
+c
+c This file is part of Octave.
+c
+c Octave is free software; you can redistribute it and/or modify
+c it under the terms of the GNU General Public License as published by
+c the Free Software Foundation; either version 3 of the License, or
+c (at your option) any later version.
+c
+c This program is distributed in the hope that it will be useful,
+c but WITHOUT ANY WARRANTY; without even the implied warranty of
+c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+c GNU General Public License for more details.
+c
+c You should have received a copy of the GNU General Public License
+c along with this software; see the file COPYING.  If not, see
+c <http://www.gnu.org/licenses/>.
+c
+
+       subroutine zrsf2csf(n,t,u,c,s)
+       integer n
+       double complex t(n,n),u(n,n)
+       double precision c(n-1),s(n-1)
+       double precision x,y,z
+       integer j
+       j = 1
+       do while (j < n)
+c apply previous rotations to rows
+         call zrcrot1(j,t(1,j),c,s)
+
+         y = t(j+1,j)
+         if (y /= 0) then
+c 2x2 block, form Givens rotation [c, i*s; i*s, c] 
+           x = t(j,j)
+           z = t(j,j+1)
+           c(j) = sqrt(abs(z/(y-z)))
+           s(j) = sign(sqrt(abs(y/(y-z))),z)
+c apply new rotation to t(j:j+1,j)
+           call zrcrot1(2,t(j,j),c(j),s(j))
+c apply all rotations to t(1:j+1,j+1)
+           call zrcrot1(j+1,t(1,j+1),c,s)
+c apply new rotation to columns j,j+1
+           call zrcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))           
+c zero subdiagonal entry, skip next row
+           t(j+1,j) = 0
+           c(j+1) = 1
+           j = j + 2
+         else
+           c(j) = 1
+           j = j + 1
+         end if
+       end do
+
+c apply rotations to last column if needed
+       if (j == n) then
+         call zrcrot1(j,t(1,j),c,s)
+       end if
+
+c apply stored rotations to all columns of u
+       do j = 1,n-1
+         if (c(j) /= 1) then
+           call zrcrot2(n,u(1,j),u(1,j+1),c(j),s(j))
+         end if
+       end do
+
+       end subroutine
+
+       subroutine zrcrot1(n,x,c,s)
+c apply rotations to a column from the left
+       integer n
+       double complex x(n), t
+       double precision c(n-1),s(n-1)
+       integer i
+       do i = 1,n-1
+         if (c(i) /= 1) then
+           t = x(i)*c(i) - x(i+1)*dcmplx(0,s(i))
+           x(i+1) = x(i+1)*c(i) - x(i)*dcmplx(0,s(i))
+           x(i) = t
+         endif
+       end do
+       end subroutine
+
+       subroutine zrcrot2(n,x,y,c,s)
+c apply a single rotation from the right to a pair of columns
+       integer n
+       double complex x(n),y(n),t
+       double precision c, s
+       integer i
+       do i = 1,n
+         t = x(i)*c + y(i)*dcmplx(0,s)
+         y(i) = y(i)*c + x(i)*dcmplx(0,s)
+         x(i) = t
+       end do
+       end subroutine
+
+
+
+
+
--- a/liboctave/ChangeLog	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/ChangeLog	Tue Jul 27 11:26:43 2010 +0200
@@ -1,3 +1,20 @@
+2010-07-27  Jaroslav Hajek  <highegg@gmail.com>
+
+	* dbleSCHUR.cc (SCHUR::SCHUR (const Matrix&, const Matrix&)): 
+	New ctor.
+	* dbleSCHUR.h: Declare it.
+	* floatSCHUR.cc (FloatSCHUR::FloatSCHUR (const FloatMatrix&, const
+	FloatMatrix&)): New ctor.
+	* floatSCHUR.h: Declare it.
+	* CmplxSCHUR.cc (ComplexSCHUR::ComplexSCHUR (const ComplexMatrix&,
+	const ComplexMatrix&),
+	ComplexSCHUR::ComplexSCHUR (const SCHUR&)): New ctors.
+	* CmplxSCHUR.h: Declare them.
+	* fCmplxSCHUR.cc (FloatComplexSCHUR::FloatComplexSCHUR 
+	(const FloatComplexMatrix&, const FloatComplexMatrix&),
+	FloatComplexSCHUR::FloatComplexSCHUR (const FloatSCHUR&)): New ctors.
+	* fCmplxSCHUR.h: Declare them.
+
 2010-07-22  Jaroslav Hajek  <highegg@gmail.com>
 
 	* dMatrix.cc (Matrix::lssolve): Fix decision test for workaround.
--- a/liboctave/CmplxSCHUR.cc	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/CmplxSCHUR.cc	Tue Jul 27 11:26:43 2010 +0200
@@ -28,6 +28,7 @@
 #include "CmplxSCHUR.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
+#include "oct-locbuf.h"
 
 extern "C"
 {
@@ -43,6 +44,9 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (zrsf2csf, ZRSF2CSF) (const octave_idx_type&,
+                                 Complex *, Complex *, double *, double *);
 }
 
 static octave_idx_type
@@ -140,3 +144,27 @@
 
   return info;
 }
+
+ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& s, 
+                            const ComplexMatrix& u)
+: schur_mat (s), unitary_mat (u)
+{
+  octave_idx_type n = s.rows ();
+  if (s.columns () != n || u.rows () != n || u.columns () != n)
+    (*current_liboctave_error_handler)
+      ("schur: inconsistent matrix dimensions");
+}
+
+ComplexSCHUR::ComplexSCHUR (const SCHUR& s)
+: schur_mat (s.schur_matrix ()), unitary_mat (s.unitary_matrix ())
+{
+  octave_idx_type n = schur_mat.rows ();
+  if (n > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (double, c, n-1);
+      OCTAVE_LOCAL_BUFFER (double, sx, n-1);
+
+      F77_XFCN (zrsf2csf, ZRSF2CSF, (n, schur_mat.fortran_vec (),
+                                     unitary_mat.fortran_vec (), c, sx));
+    }
+}
--- a/liboctave/CmplxSCHUR.h	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/CmplxSCHUR.h	Tue Jul 27 11:26:43 2010 +0200
@@ -28,6 +28,7 @@
 #include <string>
 
 #include "CMatrix.h"
+#include "dbleSCHUR.h"
 
 class
 OCTAVE_API
@@ -49,6 +50,10 @@
   ComplexSCHUR (const ComplexSCHUR& a)
     : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat) { }
 
+  ComplexSCHUR (const ComplexMatrix& s, const ComplexMatrix& u);
+
+  ComplexSCHUR (const SCHUR& s);
+
   ComplexSCHUR& operator = (const ComplexSCHUR& a)
     {
       if (this != &a)
--- a/liboctave/dbleSCHUR.cc	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/dbleSCHUR.cc	Tue Jul 27 11:26:43 2010 +0200
@@ -154,3 +154,13 @@
 
   return os;
 }
+
+SCHUR::SCHUR (const Matrix& s, const Matrix& u)
+: schur_mat (s), unitary_mat (u)
+{
+  octave_idx_type n = s.rows ();
+  if (s.columns () != n || u.rows () != n || u.columns () != n)
+    (*current_liboctave_error_handler)
+      ("schur: inconsistent matrix dimensions");
+}
+
--- a/liboctave/dbleSCHUR.h	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/dbleSCHUR.h	Tue Jul 27 11:26:43 2010 +0200
@@ -48,6 +48,8 @@
   SCHUR (const SCHUR& a)
     : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat) { }
 
+  SCHUR (const Matrix& s, const Matrix& u);
+
   SCHUR& operator = (const SCHUR& a)
     {
       if (this != &a)
--- a/liboctave/fCmplxSCHUR.cc	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/fCmplxSCHUR.cc	Tue Jul 27 11:26:43 2010 +0200
@@ -28,6 +28,7 @@
 #include "fCmplxSCHUR.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
+#include "oct-locbuf.h"
 
 extern "C"
 {
@@ -43,6 +44,9 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (crsf2csf, CRSF2CSF) (const octave_idx_type&,
+                                 FloatComplex *, FloatComplex *, float *, float *);
 }
 
 static octave_idx_type
@@ -140,3 +144,27 @@
 
   return info;
 }
+
+FloatComplexSCHUR::FloatComplexSCHUR (const FloatComplexMatrix& s, 
+                                      const FloatComplexMatrix& u)
+: schur_mat (s), unitary_mat (u)
+{
+  octave_idx_type n = s.rows ();
+  if (s.columns () != n || u.rows () != n || u.columns () != n)
+    (*current_liboctave_error_handler)
+      ("schur: inconsistent matrix dimensions");
+}
+
+FloatComplexSCHUR::FloatComplexSCHUR (const FloatSCHUR& s)
+: schur_mat (s.schur_matrix ()), unitary_mat (s.unitary_matrix ())
+{
+  octave_idx_type n = schur_mat.rows ();
+  if (n > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (float, c, n-1);
+      OCTAVE_LOCAL_BUFFER (float, sx, n-1);
+
+      F77_XFCN (crsf2csf, CRSF2CSF, (n, schur_mat.fortran_vec (),
+                                     unitary_mat.fortran_vec (), c, sx));
+    }
+}
--- a/liboctave/fCmplxSCHUR.h	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/fCmplxSCHUR.h	Tue Jul 27 11:26:43 2010 +0200
@@ -28,6 +28,7 @@
 #include <string>
 
 #include "fCMatrix.h"
+#include "floatSCHUR.h"
 
 class
 OCTAVE_API
@@ -49,6 +50,10 @@
   FloatComplexSCHUR (const FloatComplexSCHUR& a)
     : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat) { }
 
+  FloatComplexSCHUR (const FloatComplexMatrix& s, const FloatComplexMatrix& u);
+
+  FloatComplexSCHUR (const FloatSCHUR& s);
+
   FloatComplexSCHUR& operator = (const FloatComplexSCHUR& a)
     {
       if (this != &a)
--- a/liboctave/floatSCHUR.cc	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/floatSCHUR.cc	Tue Jul 27 11:26:43 2010 +0200
@@ -146,6 +146,15 @@
   return info;
 }
 
+FloatSCHUR::FloatSCHUR (const FloatMatrix& s, const FloatMatrix& u)
+: schur_mat (s), unitary_mat (u)
+{
+  octave_idx_type n = s.rows ();
+  if (s.columns () != n || u.rows () != n || u.columns () != n)
+    (*current_liboctave_error_handler)
+      ("schur: inconsistent matrix dimensions");
+}
+
 std::ostream&
 operator << (std::ostream& os, const FloatSCHUR& a)
 {
--- a/liboctave/floatSCHUR.h	Mon Jul 26 21:25:36 2010 -0700
+++ b/liboctave/floatSCHUR.h	Tue Jul 27 11:26:43 2010 +0200
@@ -48,6 +48,8 @@
   FloatSCHUR (const FloatSCHUR& a)
     : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat) { }
 
+  FloatSCHUR (const FloatMatrix& s, const FloatMatrix& u);
+
   FloatSCHUR& operator = (const FloatSCHUR& a)
     {
       if (this != &a)
--- a/src/ChangeLog	Mon Jul 26 21:25:36 2010 -0700
+++ b/src/ChangeLog	Tue Jul 27 11:26:43 2010 +0200
@@ -1,3 +1,7 @@
+2010-07-27  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/schur.cc (Frsf2csf): New DEFUN.
+
 2010-07-23  Jaroslav Hajek  <highegg@gmail.com>
 
 	* ov-base-scalar.cc (octave_base_scalar::diag): Implement here. Fix.
--- a/src/DLD-FUNCTIONS/schur.cc	Mon Jul 26 21:25:36 2010 -0700
+++ b/src/DLD-FUNCTIONS/schur.cc	Tue Jul 27 11:26:43 2010 +0200
@@ -390,3 +390,81 @@
 %!error schur ([1, 2, 3; 4, 5, 6]);
 
 */
+
+DEFUN_DLD (rsf2csf, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Function File} {[@var{U}, @var{T}] =} rsf2csf (@var{UR}, @var{TR})\n\
+Converts a real, upper quasi-triangular Schur form @var{TR} to a complex,\n\
+upper triangular Schur form @var{T}.\n\
+\n\
+Note that the following relations hold: \n\
+\n\
+@code{@var{UR} * @var{TR} * @var{UR}' = @var{U} * @var{T} * @var{U}'} and\n\
+@code{@var{U}' * @var{U}} is identity.\n\
+\n\
+Note also that U and T are not unique.\n\
+\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  if (args.length () == 2 && nargout <= 2)
+    {
+      if (! args(0).is_numeric_type ())
+        gripe_wrong_type_arg ("rsf2csf", args(0));
+      else if (! args(1).is_numeric_type ())
+        gripe_wrong_type_arg ("rsf2csf", args(1));
+      else if (args(0).is_complex_type () || args(1).is_complex_type ())
+        error ("rsf2csf: both matrices should be real");
+      else
+        {
+
+          if (args(0).is_single_type () || args(1).is_single_type ())
+            {
+              FloatMatrix u = args(0).float_matrix_value ();
+              FloatMatrix t = args(1).float_matrix_value ();
+              if (! error_state)
+                {
+                  FloatComplexSCHUR cs (FloatSCHUR (t, u));
+
+                  retval(1) = cs.schur_matrix ();
+                  retval(0) = cs.unitary_matrix ();
+                }
+            }
+          else
+            {
+              Matrix u = args(0).matrix_value ();
+              Matrix t = args(1).matrix_value ();
+              if (! error_state)
+                {
+                  ComplexSCHUR cs (SCHUR (t, u));
+
+                  retval(1) = cs.schur_matrix ();
+                  retval(0) = cs.unitary_matrix ();
+                }
+            }
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+
+%!test
+%! A = [1, 1, 1, 2; 1, 2, 1, 1; 1, 1, 3, 1; -2, 1, 1, 1];
+%! [u, t] = schur (A);
+%! [U, T] = rsf2csf (u, t);
+%! assert (norm (u * t * u' - U * T * U'), 0, 1e-12)
+%! assert (norm (A - U * T * U'), 0, 1e-12)
+
+%!test
+%! A = rand (10);
+%! [u, t] = schur (A);
+%! [U, T] = rsf2csf (u, t);
+%! assert (norm (tril (T, -1)), 0)
+%! assert (norm (U * U'), 1, 1e-14)
+
+*/