home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / liboctave / dbleSCHUR.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  3KB  |  151 lines
  1. //                                        -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <iostream.h>
  29.  
  30. #include "dbleSCHUR.h"
  31. #include "mx-inlines.cc"
  32. #include "lo-error.h"
  33. #include "f77-uscore.h"
  34.  
  35. extern "C"
  36. {
  37.   int F77_FCN (dgeesx) (const char*, const char*,
  38.             int (*)(double*, double*), const char*,
  39.             const int*, double*, const int*, int*, double*,
  40.             double*, double*, const int*, double*, double*, 
  41.             double*, const int*, int*, const int*, int*,
  42.             int*, long, long);
  43. }
  44.  
  45. static int
  46. select_ana (double *a, double *b)
  47. {
  48.    return (*a < 0.0);
  49. }
  50.  
  51. static int
  52. select_dig (double *a, double *b)
  53. {
  54.   return (hypot (*a, *b) < 1.0);
  55. }
  56.  
  57. int
  58. SCHUR::init (const Matrix& a, const char *ord)
  59. {
  60.   int a_nr = a.rows ();
  61.   int a_nc = a.cols ();
  62.   if (a_nr != a_nc)
  63.     {
  64.       (*current_liboctave_error_handler) ("SCHUR requires square matrix");
  65.       return -1;
  66.     }
  67.  
  68.   char jobvs = 'V';
  69.   char sort;
  70.  
  71.   if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd')
  72.     sort = 'S';
  73.   else
  74.     sort = 'N';
  75.  
  76.   char sense = 'N';
  77.  
  78.   int n = a_nc;
  79.   int lwork = 8 * n;
  80.   int liwork = 1;
  81.   int info;
  82.   int sdim;
  83.   double rconde;
  84.   double rcondv;
  85.  
  86.   double *s = dup (a.data (), a.length ());
  87.  
  88.   double *wr = new double [n];
  89.   double *wi = new double [n];
  90.   double *q = new double [n*n];
  91.   double *work = new double [lwork];
  92.  
  93. // These are not referenced for the non-ordered Schur routine.
  94.  
  95.   int *iwork = 0;
  96.   int *bwork = 0;
  97.   if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd')
  98.     {
  99.       iwork = new int [liwork];
  100.       bwork = new int [n];
  101.     }
  102.  
  103.   if (*ord == 'A' || *ord == 'a')
  104.     {
  105.       F77_FCN (dgeesx) (&jobvs, &sort, select_ana, &sense, &n, s, &n,
  106.             &sdim, wr, wi, q, &n, &rconde, &rcondv, work,
  107.             &lwork, iwork, &liwork, bwork, &info, 1L, 1L);
  108.     }
  109.   else if (*ord == 'D' || *ord == 'd')
  110.     {
  111.       F77_FCN (dgeesx) (&jobvs, &sort, select_dig, &sense, &n, s, &n,
  112.             &sdim, wr, wi, q, &n, &rconde, &rcondv, work,
  113.             &lwork, iwork, &liwork, bwork, &info, 1L, 1L);
  114.       
  115.     }
  116.   else
  117.     {
  118.       F77_FCN (dgeesx) (&jobvs, &sort, (void *) 0, &sense, &n, s,
  119.             &n, &sdim, wr, wi, q, &n, &rconde, &rcondv,
  120.             work, &lwork, iwork, &liwork, bwork, &info,
  121.             1L, 1L);
  122.     }
  123.  
  124.   schur_mat = Matrix (s, n, n);
  125.   unitary_mat = Matrix (q, n, n);
  126.  
  127.   delete [] wr;
  128.   delete [] wi;
  129.   delete [] work;
  130.   delete [] iwork;
  131.   delete [] bwork;
  132.  
  133.   return info;
  134. }
  135.  
  136. ostream&
  137. operator << (ostream& os, const SCHUR& a)
  138. {
  139.   os << a.schur_matrix () << "\n";
  140.   os << a.unitary_matrix () << "\n";
  141.  
  142.   return os;
  143. }
  144.  
  145. /*
  146. ;;; Local Variables: ***
  147. ;;; mode: C++ ***
  148. ;;; page-delimiter: "^/\\*" ***
  149. ;;; End: ***
  150. */
  151.