CDiagMatrix.cc

Go to the documentation of this file.
00001 // DiagMatrix manipulations.
00002 /*
00003 
00004 Copyright (C) 1994-2012 John W. Eaton
00005 Copyright (C) 2009 VZLU Prague
00006 
00007 This file is part of Octave.
00008 
00009 Octave is free software; you can redistribute it and/or modify it
00010 under the terms of the GNU General Public License as published by the
00011 Free Software Foundation; either version 3 of the License, or (at your
00012 option) any later version.
00013 
00014 Octave is distributed in the hope that it will be useful, but WITHOUT
00015 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00016 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00017 for more details.
00018 
00019 You should have received a copy of the GNU General Public License
00020 along with Octave; see the file COPYING.  If not, see
00021 <http://www.gnu.org/licenses/>.
00022 
00023 */
00024 
00025 #ifdef HAVE_CONFIG_H
00026 #include <config.h>
00027 #endif
00028 
00029 #include <iostream>
00030 
00031 #include "Array-util.h"
00032 #include "lo-error.h"
00033 #include "lo-ieee.h"
00034 #include "mx-base.h"
00035 #include "mx-inlines.cc"
00036 #include "oct-cmplx.h"
00037 
00038 // Complex Diagonal Matrix class
00039 
00040 ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a)
00041   : MDiagArray2<Complex> (a.rows (), a.cols ())
00042 {
00043   for (octave_idx_type i = 0; i < length (); i++)
00044     elem (i, i) = a.elem (i, i);
00045 }
00046 
00047 bool
00048 ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const
00049 {
00050   if (rows () != a.rows () || cols () != a.cols ())
00051     return 0;
00052 
00053   return mx_inline_equal (length (), data (), a.data ());
00054 }
00055 
00056 bool
00057 ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const
00058 {
00059   return !(*this == a);
00060 }
00061 
00062 ComplexDiagMatrix&
00063 ComplexDiagMatrix::fill (double val)
00064 {
00065   for (octave_idx_type i = 0; i < length (); i++)
00066     elem (i, i) = val;
00067   return *this;
00068 }
00069 
00070 ComplexDiagMatrix&
00071 ComplexDiagMatrix::fill (const Complex& val)
00072 {
00073   for (octave_idx_type i = 0; i < length (); i++)
00074     elem (i, i) = val;
00075   return *this;
00076 }
00077 
00078 ComplexDiagMatrix&
00079 ComplexDiagMatrix::fill (double val, octave_idx_type beg, octave_idx_type end)
00080 {
00081   if (beg < 0 || end >= length () || end < beg)
00082     {
00083       (*current_liboctave_error_handler) ("range error for fill");
00084       return *this;
00085     }
00086 
00087   for (octave_idx_type i = beg; i <= end; i++)
00088     elem (i, i) = val;
00089 
00090   return *this;
00091 }
00092 
00093 ComplexDiagMatrix&
00094 ComplexDiagMatrix::fill (const Complex& val, octave_idx_type beg, octave_idx_type end)
00095 {
00096   if (beg < 0 || end >= length () || end < beg)
00097     {
00098       (*current_liboctave_error_handler) ("range error for fill");
00099       return *this;
00100     }
00101 
00102   for (octave_idx_type i = beg; i <= end; i++)
00103     elem (i, i) = val;
00104 
00105   return *this;
00106 }
00107 
00108 ComplexDiagMatrix&
00109 ComplexDiagMatrix::fill (const ColumnVector& a)
00110 {
00111   octave_idx_type len = length ();
00112   if (a.length () != len)
00113     {
00114       (*current_liboctave_error_handler) ("range error for fill");
00115       return *this;
00116     }
00117 
00118   for (octave_idx_type i = 0; i < len; i++)
00119     elem (i, i) = a.elem (i);
00120 
00121   return *this;
00122 }
00123 
00124 ComplexDiagMatrix&
00125 ComplexDiagMatrix::fill (const ComplexColumnVector& a)
00126 {
00127   octave_idx_type len = length ();
00128   if (a.length () != len)
00129     {
00130       (*current_liboctave_error_handler) ("range error for fill");
00131       return *this;
00132     }
00133 
00134   for (octave_idx_type i = 0; i < len; i++)
00135     elem (i, i) = a.elem (i);
00136 
00137   return *this;
00138 }
00139 
00140 ComplexDiagMatrix&
00141 ComplexDiagMatrix::fill (const RowVector& a)
00142 {
00143   octave_idx_type len = length ();
00144   if (a.length () != len)
00145     {
00146       (*current_liboctave_error_handler) ("range error for fill");
00147       return *this;
00148     }
00149 
00150   for (octave_idx_type i = 0; i < len; i++)
00151     elem (i, i) = a.elem (i);
00152 
00153   return *this;
00154 }
00155 
00156 ComplexDiagMatrix&
00157 ComplexDiagMatrix::fill (const ComplexRowVector& a)
00158 {
00159   octave_idx_type len = length ();
00160   if (a.length () != len)
00161     {
00162       (*current_liboctave_error_handler) ("range error for fill");
00163       return *this;
00164     }
00165 
00166   for (octave_idx_type i = 0; i < len; i++)
00167     elem (i, i) = a.elem (i);
00168 
00169   return *this;
00170 }
00171 
00172 ComplexDiagMatrix&
00173 ComplexDiagMatrix::fill (const ColumnVector& a, octave_idx_type beg)
00174 {
00175   octave_idx_type a_len = a.length ();
00176   if (beg < 0 || beg + a_len >= length ())
00177     {
00178       (*current_liboctave_error_handler) ("range error for fill");
00179       return *this;
00180     }
00181 
00182   for (octave_idx_type i = 0; i < a_len; i++)
00183     elem (i+beg, i+beg) = a.elem (i);
00184 
00185   return *this;
00186 }
00187 
00188 ComplexDiagMatrix&
00189 ComplexDiagMatrix::fill (const ComplexColumnVector& a, octave_idx_type beg)
00190 {
00191   octave_idx_type a_len = a.length ();
00192   if (beg < 0 || beg + a_len >= length ())
00193     {
00194       (*current_liboctave_error_handler) ("range error for fill");
00195       return *this;
00196     }
00197 
00198   for (octave_idx_type i = 0; i < a_len; i++)
00199     elem (i+beg, i+beg) = a.elem (i);
00200 
00201   return *this;
00202 }
00203 
00204 ComplexDiagMatrix&
00205 ComplexDiagMatrix::fill (const RowVector& a, octave_idx_type beg)
00206 {
00207   octave_idx_type a_len = a.length ();
00208   if (beg < 0 || beg + a_len >= length ())
00209     {
00210       (*current_liboctave_error_handler) ("range error for fill");
00211       return *this;
00212     }
00213 
00214   for (octave_idx_type i = 0; i < a_len; i++)
00215     elem (i+beg, i+beg) = a.elem (i);
00216 
00217   return *this;
00218 }
00219 
00220 ComplexDiagMatrix&
00221 ComplexDiagMatrix::fill (const ComplexRowVector& a, octave_idx_type beg)
00222 {
00223   octave_idx_type a_len = a.length ();
00224   if (beg < 0 || beg + a_len >= length ())
00225     {
00226       (*current_liboctave_error_handler) ("range error for fill");
00227       return *this;
00228     }
00229 
00230   for (octave_idx_type i = 0; i < a_len; i++)
00231     elem (i+beg, i+beg) = a.elem (i);
00232 
00233   return *this;
00234 }
00235 
00236 DiagMatrix
00237 ComplexDiagMatrix::abs (void) const
00238 {
00239   return DiagMatrix (diag ().abs (), rows (), columns ());
00240 }
00241 
00242 ComplexDiagMatrix
00243 conj (const ComplexDiagMatrix& a)
00244 {
00245   return ComplexDiagMatrix (conj (a.diag ()), a.rows (), a.columns ());
00246 }
00247 
00248 // resize is the destructive analog for this one
00249 
00250 ComplexMatrix
00251 ComplexDiagMatrix::extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
00252 {
00253   if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
00254   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00255 
00256   octave_idx_type new_r = r2 - r1 + 1;
00257   octave_idx_type new_c = c2 - c1 + 1;
00258 
00259   ComplexMatrix result (new_r, new_c);
00260 
00261   for (octave_idx_type j = 0; j < new_c; j++)
00262     for (octave_idx_type i = 0; i < new_r; i++)
00263       result.elem (i, j) = elem (r1+i, c1+j);
00264 
00265   return result;
00266 }
00267 
00268 // extract row or column i.
00269 
00270 ComplexRowVector
00271 ComplexDiagMatrix::row (octave_idx_type i) const
00272 {
00273   octave_idx_type r = rows ();
00274   octave_idx_type c = cols ();
00275   if (i < 0 || i >= r)
00276     {
00277       (*current_liboctave_error_handler) ("invalid row selection");
00278       return ComplexRowVector ();
00279     }
00280 
00281   ComplexRowVector retval (c, 0.0);
00282   if (r <= c || (r > c && i < c))
00283     retval.elem (i) = elem (i, i);
00284 
00285   return retval;
00286 }
00287 
00288 ComplexRowVector
00289 ComplexDiagMatrix::row (char *s) const
00290 {
00291   if (! s)
00292     {
00293       (*current_liboctave_error_handler) ("invalid row selection");
00294       return ComplexRowVector ();
00295     }
00296 
00297   char c = *s;
00298   if (c == 'f' || c == 'F')
00299     return row (static_cast<octave_idx_type>(0));
00300   else if (c == 'l' || c == 'L')
00301     return row (rows () - 1);
00302   else
00303     {
00304       (*current_liboctave_error_handler) ("invalid row selection");
00305       return ComplexRowVector ();
00306     }
00307 }
00308 
00309 ComplexColumnVector
00310 ComplexDiagMatrix::column (octave_idx_type i) const
00311 {
00312   octave_idx_type r = rows ();
00313   octave_idx_type c = cols ();
00314   if (i < 0 || i >= c)
00315     {
00316       (*current_liboctave_error_handler) ("invalid column selection");
00317       return ComplexColumnVector ();
00318     }
00319 
00320   ComplexColumnVector retval (r, 0.0);
00321   if (r >= c || (r < c && i < r))
00322     retval.elem (i) = elem (i, i);
00323 
00324   return retval;
00325 }
00326 
00327 ComplexColumnVector
00328 ComplexDiagMatrix::column (char *s) const
00329 {
00330   if (! s)
00331     {
00332       (*current_liboctave_error_handler) ("invalid column selection");
00333       return ComplexColumnVector ();
00334     }
00335 
00336   char c = *s;
00337   if (c == 'f' || c == 'F')
00338     return column (static_cast<octave_idx_type>(0));
00339   else if (c == 'l' || c == 'L')
00340     return column (cols () - 1);
00341   else
00342     {
00343       (*current_liboctave_error_handler) ("invalid column selection");
00344       return ComplexColumnVector ();
00345     }
00346 }
00347 
00348 ComplexDiagMatrix
00349 ComplexDiagMatrix::inverse (void) const
00350 {
00351   octave_idx_type info;
00352   return inverse (info);
00353 }
00354 
00355 ComplexDiagMatrix
00356 ComplexDiagMatrix::inverse (octave_idx_type& info) const
00357 {
00358   octave_idx_type r = rows ();
00359   octave_idx_type c = cols ();
00360   if (r != c)
00361     {
00362       (*current_liboctave_error_handler) ("inverse requires square matrix");
00363       return ComplexDiagMatrix ();
00364     }
00365 
00366   ComplexDiagMatrix retval (r, c);
00367 
00368   info = 0;
00369   for (octave_idx_type i = 0; i < length (); i++)
00370     {
00371       if (elem (i, i) == 0.0)
00372         {
00373           info = -1;
00374           return *this;
00375         }
00376       else
00377         retval.elem (i, i) = 1.0 / elem (i, i);
00378     }
00379 
00380   return retval;
00381 }
00382 
00383 ComplexDiagMatrix
00384 ComplexDiagMatrix::pseudo_inverse (void) const
00385 {
00386   octave_idx_type r = rows ();
00387   octave_idx_type c = cols ();
00388   octave_idx_type len = length ();
00389 
00390   ComplexDiagMatrix retval (c, r);
00391 
00392   for (octave_idx_type i = 0; i < len; i++)
00393     {
00394       if (elem (i, i) != 0.0)
00395         retval.elem (i, i) = 1.0 / elem (i, i);
00396       else
00397         retval.elem (i, i) = 0.0;
00398     }
00399 
00400   return retval;
00401 }
00402 
00403 bool
00404 ComplexDiagMatrix::all_elements_are_real (void) const
00405 {
00406   return mx_inline_all_real (length (), data ());
00407 }
00408 
00409 // diagonal matrix by diagonal matrix -> diagonal matrix operations
00410 
00411 ComplexDiagMatrix&
00412 ComplexDiagMatrix::operator += (const DiagMatrix& a)
00413 {
00414   octave_idx_type r = rows ();
00415   octave_idx_type c = cols ();
00416 
00417   octave_idx_type a_nr = a.rows ();
00418   octave_idx_type a_nc = a.cols ();
00419 
00420   if (r != a_nr || c != a_nc)
00421     {
00422       gripe_nonconformant ("operator +=", r, c, a_nr, a_nc);
00423       return *this;
00424     }
00425 
00426   if (r == 0 || c == 0)
00427     return *this;
00428 
00429   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
00430 
00431   mx_inline_add2 (length (), d, a.data ());
00432   return *this;
00433 }
00434 
00435 ComplexDiagMatrix
00436 operator * (const ComplexDiagMatrix& a, const DiagMatrix& b)
00437 {
00438   octave_idx_type a_nr = a.rows ();
00439   octave_idx_type a_nc = a.cols ();
00440 
00441   octave_idx_type b_nr = b.rows ();
00442   octave_idx_type b_nc = b.cols ();
00443 
00444   if (a_nc != b_nr)
00445     gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
00446 
00447   ComplexDiagMatrix c (a_nr, b_nc);
00448 
00449   octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
00450 
00451   for (octave_idx_type i = 0; i < lenm; i++)
00452     c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
00453   for (octave_idx_type i = lenm; i < len; i++)
00454     c.dgxelem (i) = 0.0;
00455 
00456   return c;
00457 }
00458 
00459 ComplexDiagMatrix
00460 operator * (const DiagMatrix& a, const ComplexDiagMatrix& b)
00461 {
00462   octave_idx_type a_nr = a.rows ();
00463   octave_idx_type a_nc = a.cols ();
00464 
00465   octave_idx_type b_nr = b.rows ();
00466   octave_idx_type b_nc = b.cols ();
00467 
00468   if (a_nc != b_nr)
00469     {
00470       gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
00471       return ComplexDiagMatrix ();
00472     }
00473 
00474   if (a_nr == 0 || a_nc == 0 || b_nc == 0)
00475     return ComplexDiagMatrix (a_nr, a_nc, 0.0);
00476 
00477   ComplexDiagMatrix c (a_nr, b_nc);
00478 
00479   octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
00480 
00481   for (octave_idx_type i = 0; i < len; i++)
00482     {
00483       double a_element = a.elem (i, i);
00484       Complex b_element = b.elem (i, i);
00485 
00486       c.elem (i, i) = a_element * b_element;
00487     }
00488 
00489   return c;
00490 }
00491 
00492 ComplexDiagMatrix
00493 operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
00494 {
00495   octave_idx_type a_nr = a.rows ();
00496   octave_idx_type a_nc = a.cols ();
00497 
00498   octave_idx_type b_nr = b.rows ();
00499   octave_idx_type b_nc = b.cols ();
00500 
00501   if (a_nc != b_nr)
00502     {
00503       gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
00504       return ComplexDiagMatrix ();
00505     }
00506 
00507   if (a_nr == 0 || a_nc == 0 || b_nc == 0)
00508     return ComplexDiagMatrix (a_nr, a_nc, 0.0);
00509 
00510   ComplexDiagMatrix c (a_nr, b_nc);
00511 
00512   octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
00513 
00514   for (octave_idx_type i = 0; i < len; i++)
00515     {
00516       Complex a_element = a.elem (i, i);
00517       Complex b_element = b.elem (i, i);
00518 
00519       c.elem (i, i) = a_element * b_element;
00520     }
00521 
00522   return c;
00523 }
00524 
00525 // other operations
00526 
00527 ComplexDET
00528 ComplexDiagMatrix::determinant (void) const
00529 {
00530   ComplexDET det (1.0);
00531   if (rows () != cols ())
00532     {
00533       (*current_liboctave_error_handler) ("determinant requires square matrix");
00534       det = ComplexDET (0.0);
00535     }
00536   else
00537     {
00538       octave_idx_type len = length ();
00539       for (octave_idx_type i = 0; i < len; i++)
00540         det *= elem (i, i);
00541     }
00542 
00543   return det;
00544 }
00545 
00546 double
00547 ComplexDiagMatrix::rcond (void) const
00548 {
00549   ColumnVector av = diag (0).map<double> (std::abs);
00550   double amx = av.max (), amn = av.min ();
00551   return amx == 0 ? 0.0 : amn / amx;
00552 }
00553 
00554 // i/o
00555 
00556 std::ostream&
00557 operator << (std::ostream& os, const ComplexDiagMatrix& a)
00558 {
00559   Complex ZERO (0.0);
00560 //  int field_width = os.precision () + 7;
00561   for (octave_idx_type i = 0; i < a.rows (); i++)
00562     {
00563       for (octave_idx_type j = 0; j < a.cols (); j++)
00564         {
00565           if (i == j)
00566             os << " " /* setw (field_width) */ << a.elem (i, i);
00567           else
00568             os << " " /* setw (field_width) */ << ZERO;
00569         }
00570       os << "\n";
00571     }
00572   return os;
00573 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines