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 / dbleGEPBAL.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  5KB  |  187 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 <math.h>
  29.  
  30. #include "dbleGEPBAL.h"
  31. #include "f77-uscore.h"
  32.  
  33. extern "C"
  34. {
  35.   int F77_FCN (dgebak) (const char*, const char*, const int*, const int*,
  36.             const int*, double*, const int*, double*, const int*,
  37.             int*, long, long);
  38.  
  39.   int F77_FCN (reduce) (const int*, const int*, double*,
  40.                    const int*, double*,
  41.             int*, int*, double*, double*);
  42.  
  43.   int F77_FCN (scaleg) (const int*, const int*, double*,
  44.                    const int*, double*,
  45.             const int*, const int*, double*, double*, double*);
  46.  
  47.   int F77_FCN (gradeq) (const int*, const int*, double*,
  48.                    const int*, double*,
  49.             int*, int*, double*, double*);
  50. }
  51.  
  52. int
  53. GEPBALANCE::init (const Matrix& a, const Matrix& b, const char *balance_job)
  54. {
  55.   int a_nr = a.rows ();
  56.   int a_nc = a.cols ();
  57.   int b_nr = b.rows ();
  58.   if (a_nr != a_nc || a_nr != b_nr || b_nr != b.cols ())
  59.     {
  60.       (*current_liboctave_error_handler)
  61.     ("GEPBALANCE requires square matrices of the same size");
  62.       return -1;
  63.     }
  64.  
  65.   int n = a_nc;
  66.  
  67. // Parameters for balance call.
  68.  
  69.   int info;
  70.   int ilo;
  71.   int ihi;
  72.   double *cscale = new double [n];
  73.   double *cperm = new double [n];
  74.   Matrix wk (n, 6, 0.0);
  75.  
  76. // Back out the permutations:
  77. //
  78. // cscale contains the exponents of the column scaling factors in its 
  79. // ilo through ihi locations and the reducing column permutations in 
  80. // its first ilo-1 and its ihi+1 through n locations.
  81. //
  82. // cperm contains the column permutations applied in grading the a and b 
  83. // submatrices in its ilo through ihi locations.
  84. //
  85. // wk contains the exponents of the row scaling factors in its ilo 
  86. // through ihi locations, the reducing row permutations in its first 
  87. // ilo-1 and its ihi+1 through n locations, and the row permutations
  88. // applied in grading the a and b submatrices in its n+ilo through 
  89. // n+ihi locations.
  90.   
  91. // Copy matrices into local structure.
  92.  
  93.   balanced_a_mat = a;
  94.   balanced_b_mat = b;
  95.  
  96. // Initialize balancing matrices to identity.
  97.  
  98.   left_balancing_mat = Matrix (n, n, 0.0);
  99.   for (int i = 0; i < n; i++)
  100.     left_balancing_mat (i, i) = 1.0;
  101.  
  102.   right_balancing_mat = left_balancing_mat;
  103.  
  104. // Check for permutation option.
  105.  
  106.   if (*balance_job == 'P' || *balance_job == 'B')
  107.     {
  108.       F77_FCN (reduce) (&n, &n, balanced_a_mat.fortran_vec (),
  109.             &n, balanced_b_mat.fortran_vec (), &ilo, &ihi,
  110.             cscale, wk.fortran_vec ());
  111.     }
  112.   else
  113.     {
  114.  
  115. // Set up for scaling later.
  116.  
  117.       ilo = 1;
  118.       ihi = n;
  119.     }
  120.  
  121. // Check for scaling option.
  122.  
  123.   if ((*balance_job == 'S' || *balance_job == 'B') && ilo != ihi)
  124.     {
  125.       F77_FCN (scaleg) (&n, &n, balanced_a_mat.fortran_vec (), 
  126.             &n, balanced_b_mat.fortran_vec (), &ilo, &ihi,
  127.             cscale, cperm, wk.fortran_vec ());
  128.     }
  129.   else
  130.     {
  131.  
  132. // Set scaling data to 0's.
  133.  
  134.       for (int tmp = ilo-1; tmp < ihi; tmp++)
  135.     {
  136.       cscale[tmp] = 0.0;
  137.       wk.elem (tmp, 0) = 0.0;
  138.     }
  139.     }
  140.  
  141. // Scaleg returns exponents, not values, so...
  142.  
  143.   for (int tmp = ilo-1; tmp < ihi; tmp++)
  144.     {
  145.       cscale[tmp] = pow (2.0, cscale[tmp]);
  146.       wk.elem (tmp, 0) = pow (2.0, -wk.elem (tmp, 0));
  147.     }
  148.  
  149. // Column permutations/scaling.
  150.  
  151.   F77_FCN (dgebak) (balance_job, "R", &n, &ilo, &ihi, cscale, &n, 
  152.             right_balancing_mat.fortran_vec (), &n, &info, 1L,
  153.             1L);
  154.     
  155. // Row permutations/scaling.
  156.  
  157.   F77_FCN (dgebak) (balance_job, "L", &n, &ilo, &ihi, &wk.elem (0, 0), &n, 
  158.             left_balancing_mat.fortran_vec (), &n, &info, 1L, 1L);
  159.  
  160. // XXX FIXME XXX --- these four lines need to be added and debugged.
  161. // GEPBALANCE::init will work without them, though, so here they are.
  162.  
  163. #if 0
  164.   if ((*balance_job == 'P' || *balance_job == 'B') && ilo != ihi)
  165.     {
  166.       F77_FCN (gradeq) (&n, &n, balanced_a_mat.fortran_vec (),
  167.             &n, balanced_b_mat.fortran_vec (), &ilo, &ihi,
  168.             cperm, &wk.elem (0, 1));
  169.     }
  170. #endif
  171.  
  172. // Transpose for aa = cc*a*dd convention...
  173.   left_balancing_mat = left_balancing_mat.transpose ();
  174.  
  175.   delete [] cscale;
  176.   delete [] cperm;
  177.  
  178.   return info;
  179. }
  180.  
  181. /*
  182. ;;; Local Variables: ***
  183. ;;; mode: C++ ***
  184. ;;; page-delimiter: "^/\\*" ***
  185. ;;; End: ***
  186. */
  187.