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 / EIG.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  4KB  |  166 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 "EIG.h"
  29. #include "mx-inlines.cc"
  30. #include "lo-error.h"
  31. #include "f77-uscore.h"
  32.  
  33. extern "C"
  34. {
  35.   int F77_FCN (dgeev) (const char*, const char*, const int*, double*,
  36.                const int*, double*, double*, double*,
  37.                const int*, double*, const int*, double*,
  38.                const int*, int*, long, long);
  39.  
  40.   int F77_FCN (zgeev) (const char*, const char*, const int*, Complex*,
  41.                const int*, Complex*, Complex*, const int*,
  42.                Complex*, const int*, Complex*, const int*,
  43.                double*, int*, long, long);
  44. }
  45.  
  46. int
  47. EIG::init (const Matrix& a)
  48. {
  49.   int a_nr = a.rows ();
  50.   if (a_nr != a.cols ())
  51.     {
  52.       (*current_liboctave_error_handler) ("EIG requires square matrix");
  53.       return -1;
  54.     }
  55.  
  56.   int n = a_nr;
  57.  
  58.   int info;
  59.  
  60.   char jobvl = 'N';
  61.   char jobvr = 'V';
  62.  
  63.   double *tmp_data = dup (a.data (), a.length ());
  64.   double *wr = new double[n];
  65.   double *wi = new double[n];
  66.   Matrix vr (n, n);
  67.   double *pvr = vr.fortran_vec ();
  68.   int lwork = 8*n;
  69.   double *work = new double[lwork];
  70.  
  71.   double dummy;
  72.   int idummy = 1;
  73.  
  74.   F77_FCN (dgeev) (&jobvl, &jobvr, &n, tmp_data, &n, wr, wi, &dummy,
  75.            &idummy, pvr, &n, work, &lwork, &info, 1L, 1L);
  76.  
  77.   lambda.resize (n);
  78.   v.resize (n, n);
  79.  
  80.   for (int j = 0; j < n; j++)
  81.     {
  82.       if (wi[j] == 0.0)
  83.     {
  84.       lambda.elem (j) = Complex (wr[j]);
  85.       for (int i = 0; i < n; i++)
  86.         v.elem (i, j) = vr.elem (i, j);
  87.     }
  88.       else
  89.     {
  90.       if (j+1 >= n)
  91.         {
  92.           (*current_liboctave_error_handler) ("EIG: internal error");
  93.           return -1;
  94.         }
  95.  
  96.       for (int i = 0; i < n; i++)
  97.         {
  98.           lambda.elem (j) = Complex (wr[j], wi[j]);
  99.           lambda.elem (j+1) = Complex (wr[j+1], wi[j+1]);
  100.           double real_part = vr.elem (i, j);
  101.           double imag_part = vr.elem (i, j+1);
  102.           v.elem (i, j) = Complex (real_part, imag_part);
  103.           v.elem (i, j+1) = Complex (real_part, -imag_part);
  104.         }
  105.       j++;
  106.     }
  107.     }
  108.  
  109.   delete [] tmp_data;
  110.   delete [] wr;
  111.   delete [] wi;
  112.   delete [] work;
  113.  
  114.   return info;
  115. }
  116.  
  117. int
  118. EIG::init (const ComplexMatrix& a)
  119. {
  120.   int a_nr = a.rows ();
  121.   if (a_nr != a.cols ())
  122.     {
  123.       (*current_liboctave_error_handler) ("EIG requires square matrix");
  124.       return -1;
  125.     }
  126.  
  127.   int n = a_nr;
  128.  
  129.   int info;
  130.  
  131.   char jobvl = 'N';
  132.   char jobvr = 'V';
  133.  
  134.   lambda.resize (n);
  135.   v.resize (n, n);
  136.  
  137.   Complex *pw = lambda.fortran_vec ();
  138.   Complex *pvr = v.fortran_vec ();
  139.  
  140.   Complex *tmp_data = dup (a.data (), a.length ());
  141.  
  142.   int lwork = 8*n;
  143.   Complex *work = new Complex[lwork];
  144.   double *rwork = new double[4*n];
  145.  
  146.   Complex dummy;
  147.   int idummy = 1;
  148.  
  149.   F77_FCN (zgeev) (&jobvl, &jobvr, &n, tmp_data, &n, pw, &dummy,
  150.            &idummy, pvr, &n, work, &lwork, rwork, &info, 1L,
  151.            1L);
  152.  
  153.   delete [] tmp_data;
  154.   delete [] work;
  155.   delete [] rwork;
  156.  
  157.   return info;
  158. }
  159.  
  160. /*
  161. ;;; Local Variables: ***
  162. ;;; mode: C++ ***
  163. ;;; page-delimiter: "^/\\*" ***
  164. ;;; End: ***
  165. */
  166.