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 / dbleLU.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  2KB  |  109 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 "dbleLU.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 (dgesv) (const int*, const int*, double*, const int*,
  36.                int*, double*, const int*, int*);
  37. }
  38.  
  39. LU::LU (const Matrix& a)
  40. {
  41.   int a_nr = a.rows ();
  42.   int a_nc = a.cols ();
  43.   if (a_nr == 0 || a_nc == 0 || a_nr != a_nc)
  44.     {
  45.       (*current_liboctave_error_handler) ("LU requires square matrix");
  46.       return;
  47.     }
  48.  
  49.   int n = a_nr;
  50.  
  51.   int *ipvt = new int [n];
  52.   int *pvt = new int [n];
  53.   double *tmp_data = dup (a.data (), a.length ());
  54.   int info = 0;
  55.   int zero = 0;
  56.   double b;
  57.  
  58.   F77_FCN (dgesv) (&n, &zero, tmp_data, &n, ipvt, &b, &n, &info);
  59.  
  60.   Matrix A_fact (tmp_data, n, n);
  61.  
  62.   int i;
  63.  
  64.   for (i = 0; i < n; i++)
  65.     {
  66.       ipvt[i] -= 1;
  67.       pvt[i] = i;
  68.     }
  69.  
  70.   for (i = 0; i < n - 1; i++)
  71.     {
  72.       int k = ipvt[i];
  73.       if (k != i)
  74.     {
  75.       int tmp = pvt[k];
  76.       pvt[k] = pvt[i];
  77.       pvt[i] = tmp;
  78.     }
  79.     }
  80.  
  81.   l.resize (n, n, 0.0);
  82.   u.resize (n, n, 0.0);
  83.   p.resize (n, n, 0.0);
  84.  
  85.   for (i = 0; i < n; i++)
  86.     {
  87.       p.elem (i, pvt[i]) = 1.0;
  88.  
  89.       int j;
  90.  
  91.       l.elem (i, i) = 1.0;
  92.       for (j = 0; j < i; j++)
  93.     l.elem (i, j) = A_fact.elem (i, j);
  94.  
  95.       for (j = i; j < n; j++)
  96.     u.elem (i, j) = A_fact.elem (i, j);
  97.     }
  98.  
  99.   delete [] ipvt;
  100.   delete [] pvt;
  101. }
  102.  
  103. /*
  104. ;;; Local Variables: ***
  105. ;;; mode: C++ ***
  106. ;;; page-delimiter: "^/\\*" ***
  107. ;;; End: ***
  108. */
  109.