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 / CmplxSVD.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  2KB  |  97 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 "CmplxSVD.h"
  29. #include "mx-inlines.cc"
  30. #include "f77-uscore.h"
  31.  
  32. extern "C"
  33. {
  34.   int F77_FCN (zgesvd) (const char*, const char*, const int*,
  35.             const int*, Complex*, const int*, double*,
  36.             Complex*, const int*, Complex*, const int*,
  37.             Complex*, const int*, double*, int*, long, long);
  38. }
  39.  
  40. int
  41. ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type)
  42. {
  43.   int info;
  44.  
  45.   int m = a.rows ();
  46.   int n = a.cols ();
  47.  
  48.   Complex *tmp_data = dup (a.data (), a.length ());
  49.  
  50.   int min_mn = m < n ? m : n;
  51.   int max_mn = m > n ? m : n;
  52.  
  53.   char jobu = 'A';
  54.   char jobv = 'A';
  55.  
  56.   int ncol_u = m;
  57.   int nrow_vt = n;
  58.   int nrow_s = m;
  59.   int ncol_s = n;
  60.  
  61.   if (svd_type == SVD::economy)
  62.     {
  63.       jobu = jobv = 'S';
  64.       ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
  65.     }
  66.  
  67.   Complex *u = new Complex[m * ncol_u];
  68.   double *s_vec  = new double[min_mn];
  69.   Complex *vt = new Complex[nrow_vt * n];
  70.  
  71.   int lwork = 2*min_mn + max_mn;
  72.   Complex *work = new Complex[lwork];
  73.  
  74.   int lrwork = 5*max_mn;
  75.   double *rwork = new double[lrwork];
  76.  
  77.   F77_FCN (zgesvd) (&jobu, &jobv, &m, &n, tmp_data, &m, s_vec, u, &m,
  78.             vt, &nrow_vt, work, &lwork, rwork, &info, 1L, 1L);
  79.  
  80.   left_sm = ComplexMatrix (u, m, ncol_u);
  81.   sigma = DiagMatrix (s_vec, nrow_s, ncol_s);
  82.   ComplexMatrix vt_m (vt, nrow_vt, n);
  83.   right_sm = vt_m.hermitian ();
  84.  
  85.   delete [] tmp_data;
  86.   delete [] work;
  87.  
  88.   return info;
  89. }
  90.  
  91. /*
  92. ;;; Local Variables: ***
  93. ;;; mode: C++ ***
  94. ;;; page-delimiter: "^/\\*" ***
  95. ;;; End: ***
  96. */
  97.