fCMatrix.cc

Go to the documentation of this file.
00001 // Matrix manipulations.
00002 /*
00003 
00004 Copyright (C) 1994-2012 John W. Eaton
00005 Copyright (C) 2008-2009 Jaroslav Hajek
00006 Copyright (C) 2009 VZLU Prague, a.s.
00007 
00008 This file is part of Octave.
00009 
00010 Octave is free software; you can redistribute it and/or modify it
00011 under the terms of the GNU General Public License as published by the
00012 Free Software Foundation; either version 3 of the License, or (at your
00013 option) any later version.
00014 
00015 Octave is distributed in the hope that it will be useful, but WITHOUT
00016 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00017 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00018 for more details.
00019 
00020 You should have received a copy of the GNU General Public License
00021 along with Octave; see the file COPYING.  If not, see
00022 <http://www.gnu.org/licenses/>.
00023 
00024 */
00025 
00026 #ifdef HAVE_CONFIG_H
00027 #include <config.h>
00028 #endif
00029 
00030 #include <cfloat>
00031 
00032 #include <iostream>
00033 #include <vector>
00034 
00035 // FIXME
00036 #include <sys/types.h>
00037 
00038 #include "Array-util.h"
00039 #include "DET.h"
00040 #include "f77-fcn.h"
00041 #include "fCMatrix.h"
00042 #include "fCmplxCHOL.h"
00043 #include "fCmplxSCHUR.h"
00044 #include "fCmplxSVD.h"
00045 #include "functor.h"
00046 #include "lo-error.h"
00047 #include "lo-ieee.h"
00048 #include "lo-mappers.h"
00049 #include "lo-utils.h"
00050 #include "mx-base.h"
00051 #include "mx-fcm-fdm.h"
00052 #include "mx-fcm-fs.h"
00053 #include "mx-fdm-fcm.h"
00054 #include "mx-inlines.cc"
00055 #include "mx-op-defs.h"
00056 #include "oct-cmplx.h"
00057 #include "oct-fftw.h"
00058 #include "oct-locbuf.h"
00059 #include "oct-norm.h"
00060 
00061 // Fortran functions we call.
00062 
00063 extern "C"
00064 {
00065   F77_RET_T
00066   F77_FUNC (xilaenv, XILAENV) (const octave_idx_type&,
00067                                F77_CONST_CHAR_ARG_DECL,
00068                                F77_CONST_CHAR_ARG_DECL,
00069                                const octave_idx_type&, const octave_idx_type&,
00070                                const octave_idx_type&, const octave_idx_type&,
00071                                octave_idx_type&
00072                                F77_CHAR_ARG_LEN_DECL
00073                                F77_CHAR_ARG_LEN_DECL);
00074 
00075   F77_RET_T
00076   F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL,
00077                              const octave_idx_type&, FloatComplex*,
00078                              const octave_idx_type&, octave_idx_type&,
00079                              octave_idx_type&, float*, octave_idx_type&
00080                              F77_CHAR_ARG_LEN_DECL);
00081 
00082   F77_RET_T
00083   F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL,
00084                              F77_CONST_CHAR_ARG_DECL,
00085                              const octave_idx_type&, const octave_idx_type&,
00086                              const octave_idx_type&, float*,
00087                              const octave_idx_type&, float*,
00088                              const octave_idx_type&, octave_idx_type&
00089                              F77_CHAR_ARG_LEN_DECL
00090                              F77_CHAR_ARG_LEN_DECL);
00091 
00092   F77_RET_T
00093   F77_FUNC (cgemm, CGEMM) (F77_CONST_CHAR_ARG_DECL,
00094                            F77_CONST_CHAR_ARG_DECL,
00095                            const octave_idx_type&, const octave_idx_type&,
00096                            const octave_idx_type&, const FloatComplex&,
00097                            const FloatComplex*, const octave_idx_type&,
00098                            const FloatComplex*, const octave_idx_type&,
00099                            const FloatComplex&, FloatComplex*,
00100                            const octave_idx_type&
00101                            F77_CHAR_ARG_LEN_DECL
00102                            F77_CHAR_ARG_LEN_DECL);
00103 
00104   F77_RET_T
00105   F77_FUNC (cgemv, CGEMV) (F77_CONST_CHAR_ARG_DECL,
00106                            const octave_idx_type&, const octave_idx_type&,
00107                            const FloatComplex&, const FloatComplex*,
00108                            const octave_idx_type&, const FloatComplex*,
00109                            const octave_idx_type&, const FloatComplex&,
00110                            FloatComplex*, const octave_idx_type&
00111                            F77_CHAR_ARG_LEN_DECL);
00112 
00113   F77_RET_T
00114   F77_FUNC (xcdotu, XCDOTU) (const octave_idx_type&, const FloatComplex*,
00115                              const octave_idx_type&, const FloatComplex*,
00116                              const octave_idx_type&, FloatComplex&);
00117 
00118   F77_RET_T
00119   F77_FUNC (xcdotc, XCDOTC) (const octave_idx_type&, const FloatComplex*,
00120                              const octave_idx_type&, const FloatComplex*,
00121                              const octave_idx_type&, FloatComplex&);
00122 
00123   F77_RET_T
00124   F77_FUNC (csyrk, CSYRK) (F77_CONST_CHAR_ARG_DECL,
00125                            F77_CONST_CHAR_ARG_DECL,
00126                            const octave_idx_type&, const octave_idx_type&,
00127                            const FloatComplex&, const FloatComplex*,
00128                            const octave_idx_type&, const FloatComplex&,
00129                            FloatComplex*, const octave_idx_type&
00130                            F77_CHAR_ARG_LEN_DECL
00131                            F77_CHAR_ARG_LEN_DECL);
00132 
00133   F77_RET_T
00134   F77_FUNC (cherk, CHERK) (F77_CONST_CHAR_ARG_DECL,
00135                            F77_CONST_CHAR_ARG_DECL,
00136                            const octave_idx_type&, const octave_idx_type&,
00137                            const float&, const FloatComplex*,
00138                            const octave_idx_type&, const float&,
00139                            FloatComplex*, const octave_idx_type&
00140                            F77_CHAR_ARG_LEN_DECL
00141                            F77_CHAR_ARG_LEN_DECL);
00142 
00143   F77_RET_T
00144   F77_FUNC (cgetrf, CGETRF) (const octave_idx_type&, const octave_idx_type&,
00145                              FloatComplex*, const octave_idx_type&,
00146                              octave_idx_type*, octave_idx_type&);
00147 
00148   F77_RET_T
00149   F77_FUNC (cgetrs, CGETRS) (F77_CONST_CHAR_ARG_DECL,
00150                              const octave_idx_type&, const octave_idx_type&,
00151                              FloatComplex*, const octave_idx_type&,
00152                              const octave_idx_type*, FloatComplex*,
00153                              const octave_idx_type&, octave_idx_type&
00154                              F77_CHAR_ARG_LEN_DECL);
00155 
00156   F77_RET_T
00157   F77_FUNC (cgetri, CGETRI) (const octave_idx_type&, FloatComplex*,
00158                              const octave_idx_type&, const octave_idx_type*,
00159                              FloatComplex*, const octave_idx_type&,
00160                              octave_idx_type&);
00161 
00162   F77_RET_T
00163   F77_FUNC (cgecon, CGECON) (F77_CONST_CHAR_ARG_DECL,
00164                              const octave_idx_type&, FloatComplex*,
00165                              const octave_idx_type&, const float&, float&,
00166                              FloatComplex*, float*, octave_idx_type&
00167                              F77_CHAR_ARG_LEN_DECL);
00168 
00169   F77_RET_T
00170   F77_FUNC (cgelsy, CGELSY) (const octave_idx_type&, const octave_idx_type&,
00171                              const octave_idx_type&, FloatComplex*,
00172                              const octave_idx_type&, FloatComplex*,
00173                              const octave_idx_type&, octave_idx_type*,
00174                              float&, octave_idx_type&, FloatComplex*,
00175                              const octave_idx_type&, float*, octave_idx_type&);
00176 
00177   F77_RET_T
00178   F77_FUNC (cgelsd, CGELSD) (const octave_idx_type&, const octave_idx_type&,
00179                              const octave_idx_type&, FloatComplex*,
00180                              const octave_idx_type&, FloatComplex*,
00181                              const octave_idx_type&, float*, float&,
00182                              octave_idx_type&, FloatComplex*,
00183                              const octave_idx_type&, float*,
00184                              octave_idx_type*, octave_idx_type&);
00185 
00186   F77_RET_T
00187   F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL,
00188                              const octave_idx_type&, FloatComplex*,
00189                              const octave_idx_type&, octave_idx_type&
00190                              F77_CHAR_ARG_LEN_DECL);
00191 
00192   F77_RET_T
00193   F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL,
00194                              const octave_idx_type&, FloatComplex*,
00195                              const octave_idx_type&, const float&, float&,
00196                              FloatComplex*, float*, octave_idx_type&
00197                              F77_CHAR_ARG_LEN_DECL);
00198 
00199   F77_RET_T
00200   F77_FUNC (cpotrs, CPOTRS) (F77_CONST_CHAR_ARG_DECL,
00201                              const octave_idx_type&,
00202                              const octave_idx_type&, const FloatComplex*,
00203                              const octave_idx_type&, FloatComplex*,
00204                              const octave_idx_type&, octave_idx_type&
00205                              F77_CHAR_ARG_LEN_DECL);
00206 
00207   F77_RET_T
00208   F77_FUNC (ctrtri, CTRTRI) (F77_CONST_CHAR_ARG_DECL,
00209                              F77_CONST_CHAR_ARG_DECL,
00210                              const octave_idx_type&, const FloatComplex*,
00211                              const octave_idx_type&, octave_idx_type&
00212                              F77_CHAR_ARG_LEN_DECL
00213                              F77_CHAR_ARG_LEN_DECL);
00214 
00215   F77_RET_T
00216   F77_FUNC (ctrcon, CTRCON) (F77_CONST_CHAR_ARG_DECL,
00217                              F77_CONST_CHAR_ARG_DECL,
00218                              F77_CONST_CHAR_ARG_DECL,
00219                              const octave_idx_type&, const FloatComplex*,
00220                              const octave_idx_type&, float&, FloatComplex*,
00221                              float*, octave_idx_type&
00222                              F77_CHAR_ARG_LEN_DECL
00223                              F77_CHAR_ARG_LEN_DECL
00224                              F77_CHAR_ARG_LEN_DECL);
00225 
00226   F77_RET_T
00227   F77_FUNC (ctrtrs, CTRTRS) (F77_CONST_CHAR_ARG_DECL,
00228                              F77_CONST_CHAR_ARG_DECL,
00229                              F77_CONST_CHAR_ARG_DECL,
00230                              const octave_idx_type&, const octave_idx_type&,
00231                              const FloatComplex*, const octave_idx_type&,
00232                              FloatComplex*, const octave_idx_type&,
00233                              octave_idx_type&
00234                              F77_CHAR_ARG_LEN_DECL
00235                              F77_CHAR_ARG_LEN_DECL
00236                              F77_CHAR_ARG_LEN_DECL);
00237 
00238   F77_RET_T
00239   F77_FUNC (clartg, CLARTG) (const FloatComplex&, const FloatComplex&,
00240                              float&, FloatComplex&, FloatComplex&);
00241 
00242   F77_RET_T
00243   F77_FUNC (ctrsyl, CTRSYL) (F77_CONST_CHAR_ARG_DECL,
00244                              F77_CONST_CHAR_ARG_DECL,
00245                              const octave_idx_type&, const octave_idx_type&,
00246                              const octave_idx_type&, const FloatComplex*,
00247                              const octave_idx_type&, const FloatComplex*,
00248                              const octave_idx_type&, const FloatComplex*,
00249                              const octave_idx_type&, float&, octave_idx_type&
00250                              F77_CHAR_ARG_LEN_DECL
00251                              F77_CHAR_ARG_LEN_DECL);
00252 
00253   F77_RET_T
00254   F77_FUNC (xclange, XCLANGE) (F77_CONST_CHAR_ARG_DECL,
00255                                const octave_idx_type&, const octave_idx_type&,
00256                                const FloatComplex*, const octave_idx_type&,
00257                                float*, float&
00258                                F77_CHAR_ARG_LEN_DECL);
00259 }
00260 
00261 static const FloatComplex FloatComplex_NaN_result (octave_Float_NaN, octave_Float_NaN);
00262 
00263 // FloatComplex Matrix class
00264 
00265 FloatComplexMatrix::FloatComplexMatrix (const FloatMatrix& a)
00266   : MArray<FloatComplex> (a)
00267 {
00268 }
00269 
00270 FloatComplexMatrix::FloatComplexMatrix (const FloatRowVector& rv)
00271   : MArray<FloatComplex> (rv)
00272 {
00273 }
00274 
00275 FloatComplexMatrix::FloatComplexMatrix (const FloatColumnVector& cv)
00276   : MArray<FloatComplex> (cv)
00277 {
00278 }
00279 
00280 FloatComplexMatrix::FloatComplexMatrix (const FloatDiagMatrix& a)
00281   : MArray<FloatComplex> (a.dims (), 0.0)
00282 {
00283   for (octave_idx_type i = 0; i < a.length (); i++)
00284     elem (i, i) = a.elem (i, i);
00285 }
00286 
00287 FloatComplexMatrix::FloatComplexMatrix (const FloatComplexRowVector& rv)
00288   : MArray<FloatComplex> (rv)
00289 {
00290 }
00291 
00292 FloatComplexMatrix::FloatComplexMatrix (const FloatComplexColumnVector& cv)
00293   : MArray<FloatComplex> (cv)
00294 {
00295 }
00296 
00297 FloatComplexMatrix::FloatComplexMatrix (const FloatComplexDiagMatrix& a)
00298   : MArray<FloatComplex> (a.dims (), 0.0)
00299 {
00300   for (octave_idx_type i = 0; i < a.length (); i++)
00301     elem (i, i) = a.elem (i, i);
00302 }
00303 
00304 // FIXME -- could we use a templated mixed-type copy function
00305 // here?
00306 
00307 FloatComplexMatrix::FloatComplexMatrix (const boolMatrix& a)
00308   : MArray<FloatComplex> (a)
00309 {
00310 }
00311 
00312 FloatComplexMatrix::FloatComplexMatrix (const charMatrix& a)
00313   : MArray<FloatComplex> (a.dims (), 0.0)
00314 {
00315   for (octave_idx_type i = 0; i < a.rows (); i++)
00316     for (octave_idx_type j = 0; j < a.cols (); j++)
00317       elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
00318 }
00319 
00320 FloatComplexMatrix::FloatComplexMatrix (const FloatMatrix& re,
00321                                         const FloatMatrix& im)
00322   : MArray<FloatComplex> (re.dims ())
00323 {
00324   if (im.rows () != rows () || im.cols () != cols ())
00325     (*current_liboctave_error_handler) ("complex: internal error");
00326 
00327   octave_idx_type nel = numel ();
00328   for (octave_idx_type i = 0; i < nel; i++)
00329     xelem (i) = FloatComplex (re(i), im(i));
00330 }
00331 
00332 bool
00333 FloatComplexMatrix::operator == (const FloatComplexMatrix& a) const
00334 {
00335   if (rows () != a.rows () || cols () != a.cols ())
00336     return false;
00337 
00338   return mx_inline_equal (length (), data (), a.data ());
00339 }
00340 
00341 bool
00342 FloatComplexMatrix::operator != (const FloatComplexMatrix& a) const
00343 {
00344   return !(*this == a);
00345 }
00346 
00347 bool
00348 FloatComplexMatrix::is_hermitian (void) const
00349 {
00350   octave_idx_type nr = rows ();
00351   octave_idx_type nc = cols ();
00352 
00353   if (is_square () && nr > 0)
00354     {
00355       for (octave_idx_type i = 0; i < nr; i++)
00356         for (octave_idx_type j = i; j < nc; j++)
00357           if (elem (i, j) != conj (elem (j, i)))
00358             return false;
00359 
00360       return true;
00361     }
00362 
00363   return false;
00364 }
00365 
00366 // destructive insert/delete/reorder operations
00367 
00368 FloatComplexMatrix&
00369 FloatComplexMatrix::insert (const FloatMatrix& a, octave_idx_type r, octave_idx_type c)
00370 {
00371   octave_idx_type a_nr = a.rows ();
00372   octave_idx_type a_nc = a.cols ();
00373 
00374   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00375     {
00376       (*current_liboctave_error_handler) ("range error for insert");
00377       return *this;
00378     }
00379 
00380   if (a_nr >0 && a_nc > 0)
00381     {
00382       make_unique ();
00383 
00384       for (octave_idx_type j = 0; j < a_nc; j++)
00385         for (octave_idx_type i = 0; i < a_nr; i++)
00386           xelem (r+i, c+j) = a.elem (i, j);
00387     }
00388 
00389   return *this;
00390 }
00391 
00392 FloatComplexMatrix&
00393 FloatComplexMatrix::insert (const FloatRowVector& a, octave_idx_type r, octave_idx_type c)
00394 {
00395   octave_idx_type a_len = a.length ();
00396 
00397   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
00398     {
00399       (*current_liboctave_error_handler) ("range error for insert");
00400       return *this;
00401     }
00402 
00403   if (a_len > 0)
00404     {
00405       make_unique ();
00406 
00407       for (octave_idx_type i = 0; i < a_len; i++)
00408         xelem (r, c+i) = a.elem (i);
00409     }
00410 
00411   return *this;
00412 }
00413 
00414 FloatComplexMatrix&
00415 FloatComplexMatrix::insert (const FloatColumnVector& a, octave_idx_type r, octave_idx_type c)
00416 {
00417   octave_idx_type a_len = a.length ();
00418 
00419   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
00420     {
00421       (*current_liboctave_error_handler) ("range error for insert");
00422       return *this;
00423     }
00424 
00425   if (a_len > 0)
00426     {
00427       make_unique ();
00428 
00429       for (octave_idx_type i = 0; i < a_len; i++)
00430         xelem (r+i, c) = a.elem (i);
00431     }
00432 
00433   return *this;
00434 }
00435 
00436 FloatComplexMatrix&
00437 FloatComplexMatrix::insert (const FloatDiagMatrix& a, octave_idx_type r, octave_idx_type c)
00438 {
00439   octave_idx_type a_nr = a.rows ();
00440   octave_idx_type a_nc = a.cols ();
00441 
00442   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00443     {
00444       (*current_liboctave_error_handler) ("range error for insert");
00445       return *this;
00446     }
00447 
00448   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
00449 
00450   octave_idx_type a_len = a.length ();
00451 
00452   if (a_len > 0)
00453     {
00454       make_unique ();
00455 
00456       for (octave_idx_type i = 0; i < a_len; i++)
00457         xelem (r+i, c+i) = a.elem (i, i);
00458     }
00459 
00460   return *this;
00461 }
00462 
00463 FloatComplexMatrix&
00464 FloatComplexMatrix::insert (const FloatComplexMatrix& a, octave_idx_type r, octave_idx_type c)
00465 {
00466   Array<FloatComplex>::insert (a, r, c);
00467   return *this;
00468 }
00469 
00470 FloatComplexMatrix&
00471 FloatComplexMatrix::insert (const FloatComplexRowVector& a, octave_idx_type r, octave_idx_type c)
00472 {
00473   octave_idx_type a_len = a.length ();
00474   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
00475     {
00476       (*current_liboctave_error_handler) ("range error for insert");
00477       return *this;
00478     }
00479 
00480   for (octave_idx_type i = 0; i < a_len; i++)
00481     elem (r, c+i) = a.elem (i);
00482 
00483   return *this;
00484 }
00485 
00486 FloatComplexMatrix&
00487 FloatComplexMatrix::insert (const FloatComplexColumnVector& a, octave_idx_type r, octave_idx_type c)
00488 {
00489   octave_idx_type a_len = a.length ();
00490 
00491   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
00492     {
00493       (*current_liboctave_error_handler) ("range error for insert");
00494       return *this;
00495     }
00496 
00497   if (a_len > 0)
00498     {
00499       make_unique ();
00500 
00501       for (octave_idx_type i = 0; i < a_len; i++)
00502         xelem (r+i, c) = a.elem (i);
00503     }
00504 
00505   return *this;
00506 }
00507 
00508 FloatComplexMatrix&
00509 FloatComplexMatrix::insert (const FloatComplexDiagMatrix& a, octave_idx_type r, octave_idx_type c)
00510 {
00511   octave_idx_type a_nr = a.rows ();
00512   octave_idx_type a_nc = a.cols ();
00513 
00514   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00515     {
00516       (*current_liboctave_error_handler) ("range error for insert");
00517       return *this;
00518     }
00519 
00520   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
00521 
00522   octave_idx_type a_len = a.length ();
00523 
00524   if (a_len > 0)
00525     {
00526       make_unique ();
00527 
00528       for (octave_idx_type i = 0; i < a_len; i++)
00529         xelem (r+i, c+i) = a.elem (i, i);
00530     }
00531 
00532   return *this;
00533 }
00534 
00535 FloatComplexMatrix&
00536 FloatComplexMatrix::fill (float val)
00537 {
00538   octave_idx_type nr = rows ();
00539   octave_idx_type nc = cols ();
00540 
00541   if (nr > 0 && nc > 0)
00542     {
00543       make_unique ();
00544 
00545       for (octave_idx_type j = 0; j < nc; j++)
00546         for (octave_idx_type i = 0; i < nr; i++)
00547           xelem (i, j) = val;
00548     }
00549 
00550   return *this;
00551 }
00552 
00553 FloatComplexMatrix&
00554 FloatComplexMatrix::fill (const FloatComplex& val)
00555 {
00556   octave_idx_type nr = rows ();
00557   octave_idx_type nc = cols ();
00558 
00559   if (nr > 0 && nc > 0)
00560     {
00561       make_unique ();
00562 
00563       for (octave_idx_type j = 0; j < nc; j++)
00564         for (octave_idx_type i = 0; i < nr; i++)
00565           xelem (i, j) = val;
00566     }
00567 
00568   return *this;
00569 }
00570 
00571 FloatComplexMatrix&
00572 FloatComplexMatrix::fill (float val, octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2)
00573 {
00574   octave_idx_type nr = rows ();
00575   octave_idx_type nc = cols ();
00576 
00577   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
00578       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
00579     {
00580       (*current_liboctave_error_handler) ("range error for fill");
00581       return *this;
00582     }
00583 
00584   if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
00585   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00586 
00587   if (r2 >= r1 && c2 >= c1)
00588     {
00589       make_unique ();
00590 
00591       for (octave_idx_type j = c1; j <= c2; j++)
00592         for (octave_idx_type i = r1; i <= r2; i++)
00593           xelem (i, j) = val;
00594     }
00595 
00596   return *this;
00597 }
00598 
00599 FloatComplexMatrix&
00600 FloatComplexMatrix::fill (const FloatComplex& val, octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2)
00601 {
00602   octave_idx_type nr = rows ();
00603   octave_idx_type nc = cols ();
00604 
00605   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
00606       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
00607     {
00608       (*current_liboctave_error_handler) ("range error for fill");
00609       return *this;
00610     }
00611 
00612   if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
00613   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00614 
00615   if (r2 >= r1 && c2 >=c1)
00616     {
00617       make_unique ();
00618 
00619       for (octave_idx_type j = c1; j <= c2; j++)
00620         for (octave_idx_type i = r1; i <= r2; i++)
00621           xelem (i, j) = val;
00622     }
00623 
00624   return *this;
00625 }
00626 
00627 FloatComplexMatrix
00628 FloatComplexMatrix::append (const FloatMatrix& a) const
00629 {
00630   octave_idx_type nr = rows ();
00631   octave_idx_type nc = cols ();
00632   if (nr != a.rows ())
00633     {
00634       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00635       return *this;
00636     }
00637 
00638   octave_idx_type nc_insert = nc;
00639   FloatComplexMatrix retval (nr, nc + a.cols ());
00640   retval.insert (*this, 0, 0);
00641   retval.insert (a, 0, nc_insert);
00642   return retval;
00643 }
00644 
00645 FloatComplexMatrix
00646 FloatComplexMatrix::append (const FloatRowVector& a) const
00647 {
00648   octave_idx_type nr = rows ();
00649   octave_idx_type nc = cols ();
00650   if (nr != 1)
00651     {
00652       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00653       return *this;
00654     }
00655 
00656   octave_idx_type nc_insert = nc;
00657   FloatComplexMatrix retval (nr, nc + a.length ());
00658   retval.insert (*this, 0, 0);
00659   retval.insert (a, 0, nc_insert);
00660   return retval;
00661 }
00662 
00663 FloatComplexMatrix
00664 FloatComplexMatrix::append (const FloatColumnVector& a) const
00665 {
00666   octave_idx_type nr = rows ();
00667   octave_idx_type nc = cols ();
00668   if (nr != a.length ())
00669     {
00670       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00671       return *this;
00672     }
00673 
00674   octave_idx_type nc_insert = nc;
00675   FloatComplexMatrix retval (nr, nc + 1);
00676   retval.insert (*this, 0, 0);
00677   retval.insert (a, 0, nc_insert);
00678   return retval;
00679 }
00680 
00681 FloatComplexMatrix
00682 FloatComplexMatrix::append (const FloatDiagMatrix& a) const
00683 {
00684   octave_idx_type nr = rows ();
00685   octave_idx_type nc = cols ();
00686   if (nr != a.rows ())
00687     {
00688       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00689       return *this;
00690     }
00691 
00692   octave_idx_type nc_insert = nc;
00693   FloatComplexMatrix retval (nr, nc + a.cols ());
00694   retval.insert (*this, 0, 0);
00695   retval.insert (a, 0, nc_insert);
00696   return retval;
00697 }
00698 
00699 FloatComplexMatrix
00700 FloatComplexMatrix::append (const FloatComplexMatrix& a) const
00701 {
00702   octave_idx_type nr = rows ();
00703   octave_idx_type nc = cols ();
00704   if (nr != a.rows ())
00705     {
00706       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00707       return *this;
00708     }
00709 
00710   octave_idx_type nc_insert = nc;
00711   FloatComplexMatrix retval (nr, nc + a.cols ());
00712   retval.insert (*this, 0, 0);
00713   retval.insert (a, 0, nc_insert);
00714   return retval;
00715 }
00716 
00717 FloatComplexMatrix
00718 FloatComplexMatrix::append (const FloatComplexRowVector& a) const
00719 {
00720   octave_idx_type nr = rows ();
00721   octave_idx_type nc = cols ();
00722   if (nr != 1)
00723     {
00724       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00725       return *this;
00726     }
00727 
00728   octave_idx_type nc_insert = nc;
00729   FloatComplexMatrix retval (nr, nc + a.length ());
00730   retval.insert (*this, 0, 0);
00731   retval.insert (a, 0, nc_insert);
00732   return retval;
00733 }
00734 
00735 FloatComplexMatrix
00736 FloatComplexMatrix::append (const FloatComplexColumnVector& a) const
00737 {
00738   octave_idx_type nr = rows ();
00739   octave_idx_type nc = cols ();
00740   if (nr != a.length ())
00741     {
00742       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00743       return *this;
00744     }
00745 
00746   octave_idx_type nc_insert = nc;
00747   FloatComplexMatrix retval (nr, nc + 1);
00748   retval.insert (*this, 0, 0);
00749   retval.insert (a, 0, nc_insert);
00750   return retval;
00751 }
00752 
00753 FloatComplexMatrix
00754 FloatComplexMatrix::append (const FloatComplexDiagMatrix& a) const
00755 {
00756   octave_idx_type nr = rows ();
00757   octave_idx_type nc = cols ();
00758   if (nr != a.rows ())
00759     {
00760       (*current_liboctave_error_handler) ("row dimension mismatch for append");
00761       return *this;
00762     }
00763 
00764   octave_idx_type nc_insert = nc;
00765   FloatComplexMatrix retval (nr, nc + a.cols ());
00766   retval.insert (*this, 0, 0);
00767   retval.insert (a, 0, nc_insert);
00768   return retval;
00769 }
00770 
00771 FloatComplexMatrix
00772 FloatComplexMatrix::stack (const FloatMatrix& a) const
00773 {
00774   octave_idx_type nr = rows ();
00775   octave_idx_type nc = cols ();
00776   if (nc != a.cols ())
00777     {
00778       (*current_liboctave_error_handler)
00779         ("column dimension mismatch for stack");
00780       return *this;
00781     }
00782 
00783   octave_idx_type nr_insert = nr;
00784   FloatComplexMatrix retval (nr + a.rows (), nc);
00785   retval.insert (*this, 0, 0);
00786   retval.insert (a, nr_insert, 0);
00787   return retval;
00788 }
00789 
00790 FloatComplexMatrix
00791 FloatComplexMatrix::stack (const FloatRowVector& a) const
00792 {
00793   octave_idx_type nr = rows ();
00794   octave_idx_type nc = cols ();
00795   if (nc != a.length ())
00796     {
00797       (*current_liboctave_error_handler)
00798         ("column dimension mismatch for stack");
00799       return *this;
00800     }
00801 
00802   octave_idx_type nr_insert = nr;
00803   FloatComplexMatrix retval (nr + 1, nc);
00804   retval.insert (*this, 0, 0);
00805   retval.insert (a, nr_insert, 0);
00806   return retval;
00807 }
00808 
00809 FloatComplexMatrix
00810 FloatComplexMatrix::stack (const FloatColumnVector& a) const
00811 {
00812   octave_idx_type nr = rows ();
00813   octave_idx_type nc = cols ();
00814   if (nc != 1)
00815     {
00816       (*current_liboctave_error_handler)
00817         ("column dimension mismatch for stack");
00818       return *this;
00819     }
00820 
00821   octave_idx_type nr_insert = nr;
00822   FloatComplexMatrix retval (nr + a.length (), nc);
00823   retval.insert (*this, 0, 0);
00824   retval.insert (a, nr_insert, 0);
00825   return retval;
00826 }
00827 
00828 FloatComplexMatrix
00829 FloatComplexMatrix::stack (const FloatDiagMatrix& a) const
00830 {
00831   octave_idx_type nr = rows ();
00832   octave_idx_type nc = cols ();
00833   if (nc != a.cols ())
00834     {
00835       (*current_liboctave_error_handler)
00836         ("column dimension mismatch for stack");
00837       return *this;
00838     }
00839 
00840   octave_idx_type nr_insert = nr;
00841   FloatComplexMatrix retval (nr + a.rows (), nc);
00842   retval.insert (*this, 0, 0);
00843   retval.insert (a, nr_insert, 0);
00844   return retval;
00845 }
00846 
00847 FloatComplexMatrix
00848 FloatComplexMatrix::stack (const FloatComplexMatrix& a) const
00849 {
00850   octave_idx_type nr = rows ();
00851   octave_idx_type nc = cols ();
00852   if (nc != a.cols ())
00853     {
00854       (*current_liboctave_error_handler)
00855         ("column dimension mismatch for stack");
00856       return *this;
00857     }
00858 
00859   octave_idx_type nr_insert = nr;
00860   FloatComplexMatrix retval (nr + a.rows (), nc);
00861   retval.insert (*this, 0, 0);
00862   retval.insert (a, nr_insert, 0);
00863   return retval;
00864 }
00865 
00866 FloatComplexMatrix
00867 FloatComplexMatrix::stack (const FloatComplexRowVector& a) const
00868 {
00869   octave_idx_type nr = rows ();
00870   octave_idx_type nc = cols ();
00871   if (nc != a.length ())
00872     {
00873       (*current_liboctave_error_handler)
00874         ("column dimension mismatch for stack");
00875       return *this;
00876     }
00877 
00878   octave_idx_type nr_insert = nr;
00879   FloatComplexMatrix retval (nr + 1, nc);
00880   retval.insert (*this, 0, 0);
00881   retval.insert (a, nr_insert, 0);
00882   return retval;
00883 }
00884 
00885 FloatComplexMatrix
00886 FloatComplexMatrix::stack (const FloatComplexColumnVector& a) const
00887 {
00888   octave_idx_type nr = rows ();
00889   octave_idx_type nc = cols ();
00890   if (nc != 1)
00891     {
00892       (*current_liboctave_error_handler)
00893         ("column dimension mismatch for stack");
00894       return *this;
00895     }
00896 
00897   octave_idx_type nr_insert = nr;
00898   FloatComplexMatrix retval (nr + a.length (), nc);
00899   retval.insert (*this, 0, 0);
00900   retval.insert (a, nr_insert, 0);
00901   return retval;
00902 }
00903 
00904 FloatComplexMatrix
00905 FloatComplexMatrix::stack (const FloatComplexDiagMatrix& a) const
00906 {
00907   octave_idx_type nr = rows ();
00908   octave_idx_type nc = cols ();
00909   if (nc != a.cols ())
00910     {
00911       (*current_liboctave_error_handler)
00912         ("column dimension mismatch for stack");
00913       return *this;
00914     }
00915 
00916   octave_idx_type nr_insert = nr;
00917   FloatComplexMatrix retval (nr + a.rows (), nc);
00918   retval.insert (*this, 0, 0);
00919   retval.insert (a, nr_insert, 0);
00920   return retval;
00921 }
00922 
00923 FloatComplexMatrix
00924 conj (const FloatComplexMatrix& a)
00925 {
00926   return do_mx_unary_map<FloatComplex, FloatComplex, std::conj<float> > (a);
00927 }
00928 
00929 // resize is the destructive equivalent for this one
00930 
00931 FloatComplexMatrix
00932 FloatComplexMatrix::extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
00933 {
00934   if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
00935   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00936 
00937   return index (idx_vector (r1, r2+1), idx_vector (c1, c2+1));
00938 }
00939 
00940 FloatComplexMatrix
00941 FloatComplexMatrix::extract_n (octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const
00942 {
00943   return index (idx_vector (r1, r1 + nr), idx_vector (c1, c1 + nc));
00944 }
00945 
00946 // extract row or column i.
00947 
00948 FloatComplexRowVector
00949 FloatComplexMatrix::row (octave_idx_type i) const
00950 {
00951   return index (idx_vector (i), idx_vector::colon);
00952 }
00953 
00954 FloatComplexColumnVector
00955 FloatComplexMatrix::column (octave_idx_type i) const
00956 {
00957   return index (idx_vector::colon, idx_vector (i));
00958 }
00959 
00960 FloatComplexMatrix
00961 FloatComplexMatrix::inverse (void) const
00962 {
00963   octave_idx_type info;
00964   float rcon;
00965   MatrixType mattype (*this);
00966   return inverse (mattype, info, rcon, 0, 0);
00967 }
00968 
00969 FloatComplexMatrix
00970 FloatComplexMatrix::inverse (octave_idx_type& info) const
00971 {
00972   float rcon;
00973   MatrixType mattype (*this);
00974   return inverse (mattype, info, rcon, 0, 0);
00975 }
00976 
00977 FloatComplexMatrix
00978 FloatComplexMatrix::inverse (octave_idx_type& info, float& rcon, int force,
00979                         int calc_cond) const
00980 {
00981   MatrixType mattype (*this);
00982   return inverse (mattype, info, rcon, force, calc_cond);
00983 }
00984 
00985 FloatComplexMatrix
00986 FloatComplexMatrix::inverse (MatrixType &mattype) const
00987 {
00988   octave_idx_type info;
00989   float rcon;
00990   return inverse (mattype, info, rcon, 0, 0);
00991 }
00992 
00993 FloatComplexMatrix
00994 FloatComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const
00995 {
00996   float rcon;
00997   return inverse (mattype, info, rcon, 0, 0);
00998 }
00999 
01000 FloatComplexMatrix
01001 FloatComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info,
01002                          float& rcon, int force, int calc_cond) const
01003 {
01004   FloatComplexMatrix retval;
01005 
01006   octave_idx_type nr = rows ();
01007   octave_idx_type nc = cols ();
01008 
01009   if (nr != nc || nr == 0 || nc == 0)
01010     (*current_liboctave_error_handler) ("inverse requires square matrix");
01011   else
01012     {
01013       int typ = mattype.type ();
01014       char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
01015       char udiag = 'N';
01016       retval = *this;
01017       FloatComplex *tmp_data = retval.fortran_vec ();
01018 
01019       F77_XFCN (ctrtri, CTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
01020                                  F77_CONST_CHAR_ARG2 (&udiag, 1),
01021                                  nr, tmp_data, nr, info
01022                                  F77_CHAR_ARG_LEN (1)
01023                                  F77_CHAR_ARG_LEN (1)));
01024 
01025       // Throw-away extra info LAPACK gives so as to not change output.
01026       rcon = 0.0;
01027       if (info != 0)
01028         info = -1;
01029       else if (calc_cond)
01030         {
01031           octave_idx_type ztrcon_info = 0;
01032           char job = '1';
01033 
01034           OCTAVE_LOCAL_BUFFER (FloatComplex, cwork, 2*nr);
01035           OCTAVE_LOCAL_BUFFER (float, rwork, nr);
01036 
01037           F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
01038                                      F77_CONST_CHAR_ARG2 (&uplo, 1),
01039                                      F77_CONST_CHAR_ARG2 (&udiag, 1),
01040                                      nr, tmp_data, nr, rcon,
01041                                      cwork, rwork, ztrcon_info
01042                                      F77_CHAR_ARG_LEN (1)
01043                                      F77_CHAR_ARG_LEN (1)
01044                                      F77_CHAR_ARG_LEN (1)));
01045 
01046           if (ztrcon_info != 0)
01047             info = -1;
01048         }
01049 
01050       if (info == -1 && ! force)
01051         retval = *this; // Restore matrix contents.
01052     }
01053 
01054   return retval;
01055 }
01056 
01057 FloatComplexMatrix
01058 FloatComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info,
01059                          float& rcon, int force, int calc_cond) const
01060 {
01061   FloatComplexMatrix retval;
01062 
01063   octave_idx_type nr = rows ();
01064   octave_idx_type nc = cols ();
01065 
01066   if (nr != nc)
01067     (*current_liboctave_error_handler) ("inverse requires square matrix");
01068   else
01069     {
01070       Array<octave_idx_type> ipvt (dim_vector (nr, 1));
01071       octave_idx_type *pipvt = ipvt.fortran_vec ();
01072 
01073       retval = *this;
01074       FloatComplex *tmp_data = retval.fortran_vec ();
01075 
01076       Array<FloatComplex> z (dim_vector (1, 1));
01077       octave_idx_type lwork = -1;
01078 
01079       // Query the optimum work array size.
01080 
01081       F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt,
01082                                  z.fortran_vec (), lwork, info));
01083 
01084       lwork = static_cast<octave_idx_type> (std::real(z(0)));
01085       lwork = (lwork <  2 *nc ? 2*nc : lwork);
01086       z.resize (dim_vector (lwork, 1));
01087       FloatComplex *pz = z.fortran_vec ();
01088 
01089       info = 0;
01090 
01091       // Calculate the norm of the matrix, for later use.
01092       float anorm;
01093       if (calc_cond)
01094         anorm  = retval.abs().sum().row(static_cast<octave_idx_type>(0)).max();
01095 
01096       F77_XFCN (cgetrf, CGETRF, (nc, nc, tmp_data, nr, pipvt, info));
01097 
01098       // Throw-away extra info LAPACK gives so as to not change output.
01099       rcon = 0.0;
01100       if (info != 0)
01101         info = -1;
01102       else if (calc_cond)
01103         {
01104           // Now calculate the condition number for non-singular matrix.
01105           octave_idx_type zgecon_info = 0;
01106           char job = '1';
01107           Array<float> rz (dim_vector (2 * nc, 1));
01108           float *prz = rz.fortran_vec ();
01109           F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01110                                      nc, tmp_data, nr, anorm,
01111                                      rcon, pz, prz, zgecon_info
01112                                      F77_CHAR_ARG_LEN (1)));
01113 
01114           if (zgecon_info != 0)
01115             info = -1;
01116         }
01117 
01118       if (info == -1 && ! force)
01119         retval = *this;  // Restore contents.
01120       else
01121         {
01122           octave_idx_type zgetri_info = 0;
01123 
01124           F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt,
01125                                      pz, lwork, zgetri_info));
01126 
01127           if (zgetri_info != 0)
01128             info = -1;
01129         }
01130 
01131       if (info != 0)
01132         mattype.mark_as_rectangular();
01133     }
01134 
01135   return retval;
01136 }
01137 
01138 FloatComplexMatrix
01139 FloatComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info,
01140                         float& rcon, int force, int calc_cond) const
01141 {
01142   int typ = mattype.type (false);
01143   FloatComplexMatrix ret;
01144 
01145   if (typ == MatrixType::Unknown)
01146     typ = mattype.type (*this);
01147 
01148   if (typ == MatrixType::Upper || typ == MatrixType::Lower)
01149     ret = tinverse (mattype, info, rcon, force, calc_cond);
01150   else
01151     {
01152       if (mattype.is_hermitian ())
01153         {
01154           FloatComplexCHOL chol (*this, info, calc_cond);
01155           if (info == 0)
01156             {
01157               if (calc_cond)
01158                 rcon = chol.rcond();
01159               else
01160                 rcon = 1.0;
01161               ret = chol.inverse ();
01162             }
01163           else
01164             mattype.mark_as_unsymmetric ();
01165         }
01166 
01167       if (!mattype.is_hermitian ())
01168         ret = finverse(mattype, info, rcon, force, calc_cond);
01169 
01170       if ((mattype.is_hermitian () || calc_cond) && rcon == 0.)
01171         ret = FloatComplexMatrix (rows (), columns (), FloatComplex (octave_Float_Inf, 0.));
01172     }
01173 
01174   return ret;
01175 }
01176 
01177 FloatComplexMatrix
01178 FloatComplexMatrix::pseudo_inverse (float tol) const
01179 {
01180   FloatComplexMatrix retval;
01181 
01182   FloatComplexSVD result (*this, SVD::economy);
01183 
01184   FloatDiagMatrix S = result.singular_values ();
01185   FloatComplexMatrix U = result.left_singular_matrix ();
01186   FloatComplexMatrix V = result.right_singular_matrix ();
01187 
01188   FloatColumnVector sigma = S.diag ();
01189 
01190   octave_idx_type r = sigma.length () - 1;
01191   octave_idx_type nr = rows ();
01192   octave_idx_type nc = cols ();
01193 
01194   if (tol <= 0.0)
01195     {
01196       if (nr > nc)
01197         tol = nr * sigma.elem (0) * DBL_EPSILON;
01198       else
01199         tol = nc * sigma.elem (0) * DBL_EPSILON;
01200     }
01201 
01202   while (r >= 0 && sigma.elem (r) < tol)
01203     r--;
01204 
01205   if (r < 0)
01206     retval = FloatComplexMatrix (nc, nr, 0.0);
01207   else
01208     {
01209       FloatComplexMatrix Ur = U.extract (0, 0, nr-1, r);
01210       FloatDiagMatrix D = FloatDiagMatrix (sigma.extract (0, r)) . inverse ();
01211       FloatComplexMatrix Vr = V.extract (0, 0, nc-1, r);
01212       retval = Vr * D * Ur.hermitian ();
01213     }
01214 
01215   return retval;
01216 }
01217 
01218 #if defined (HAVE_FFTW)
01219 
01220 FloatComplexMatrix
01221 FloatComplexMatrix::fourier (void) const
01222 {
01223   size_t nr = rows ();
01224   size_t nc = cols ();
01225 
01226   FloatComplexMatrix retval (nr, nc);
01227 
01228   size_t npts, nsamples;
01229 
01230   if (nr == 1 || nc == 1)
01231     {
01232       npts = nr > nc ? nr : nc;
01233       nsamples = 1;
01234     }
01235   else
01236     {
01237       npts = nr;
01238       nsamples = nc;
01239     }
01240 
01241   const FloatComplex *in (data ());
01242   FloatComplex *out (retval.fortran_vec ());
01243 
01244   octave_fftw::fft (in, out, npts, nsamples);
01245 
01246   return retval;
01247 }
01248 
01249 FloatComplexMatrix
01250 FloatComplexMatrix::ifourier (void) const
01251 {
01252   size_t nr = rows ();
01253   size_t nc = cols ();
01254 
01255   FloatComplexMatrix retval (nr, nc);
01256 
01257   size_t npts, nsamples;
01258 
01259   if (nr == 1 || nc == 1)
01260     {
01261       npts = nr > nc ? nr : nc;
01262       nsamples = 1;
01263     }
01264   else
01265     {
01266       npts = nr;
01267       nsamples = nc;
01268     }
01269 
01270   const FloatComplex *in (data ());
01271   FloatComplex *out (retval.fortran_vec ());
01272 
01273   octave_fftw::ifft (in, out, npts, nsamples);
01274 
01275   return retval;
01276 }
01277 
01278 FloatComplexMatrix
01279 FloatComplexMatrix::fourier2d (void) const
01280 {
01281   dim_vector dv(rows (), cols ());
01282 
01283   FloatComplexMatrix retval (rows (), cols ());
01284   const FloatComplex *in (data ());
01285   FloatComplex *out (retval.fortran_vec ());
01286 
01287   octave_fftw::fftNd (in, out, 2, dv);
01288 
01289   return retval;
01290 }
01291 
01292 FloatComplexMatrix
01293 FloatComplexMatrix::ifourier2d (void) const
01294 {
01295   dim_vector dv(rows (), cols ());
01296 
01297   FloatComplexMatrix retval (rows (), cols ());
01298   const FloatComplex *in (data ());
01299   FloatComplex *out (retval.fortran_vec ());
01300 
01301   octave_fftw::ifftNd (in, out, 2, dv);
01302 
01303   return retval;
01304 }
01305 
01306 #else
01307 
01308 extern "C"
01309 {
01310   F77_RET_T
01311   F77_FUNC (cffti, CFFTI) (const octave_idx_type&, FloatComplex*);
01312 
01313   F77_RET_T
01314   F77_FUNC (cfftf, CFFTF) (const octave_idx_type&, FloatComplex*, FloatComplex*);
01315 
01316   F77_RET_T
01317   F77_FUNC (cfftb, CFFTB) (const octave_idx_type&, FloatComplex*, FloatComplex*);
01318 }
01319 
01320 FloatComplexMatrix
01321 FloatComplexMatrix::fourier (void) const
01322 {
01323   FloatComplexMatrix retval;
01324 
01325   octave_idx_type nr = rows ();
01326   octave_idx_type nc = cols ();
01327 
01328   octave_idx_type npts, nsamples;
01329 
01330   if (nr == 1 || nc == 1)
01331     {
01332       npts = nr > nc ? nr : nc;
01333       nsamples = 1;
01334     }
01335   else
01336     {
01337       npts = nr;
01338       nsamples = nc;
01339     }
01340 
01341   octave_idx_type nn = 4*npts+15;
01342 
01343   Array<FloatComplex> wsave (dim_vector (nn, 1));
01344   FloatComplex *pwsave = wsave.fortran_vec ();
01345 
01346   retval = *this;
01347   FloatComplex *tmp_data = retval.fortran_vec ();
01348 
01349   F77_FUNC (cffti, CFFTI) (npts, pwsave);
01350 
01351   for (octave_idx_type j = 0; j < nsamples; j++)
01352     {
01353       octave_quit ();
01354 
01355       F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
01356     }
01357 
01358   return retval;
01359 }
01360 
01361 FloatComplexMatrix
01362 FloatComplexMatrix::ifourier (void) const
01363 {
01364   FloatComplexMatrix retval;
01365 
01366   octave_idx_type nr = rows ();
01367   octave_idx_type nc = cols ();
01368 
01369   octave_idx_type npts, nsamples;
01370 
01371   if (nr == 1 || nc == 1)
01372     {
01373       npts = nr > nc ? nr : nc;
01374       nsamples = 1;
01375     }
01376   else
01377     {
01378       npts = nr;
01379       nsamples = nc;
01380     }
01381 
01382   octave_idx_type nn = 4*npts+15;
01383 
01384   Array<FloatComplex> wsave (dim_vector (nn, 1));
01385   FloatComplex *pwsave = wsave.fortran_vec ();
01386 
01387   retval = *this;
01388   FloatComplex *tmp_data = retval.fortran_vec ();
01389 
01390   F77_FUNC (cffti, CFFTI) (npts, pwsave);
01391 
01392   for (octave_idx_type j = 0; j < nsamples; j++)
01393     {
01394       octave_quit ();
01395 
01396       F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
01397     }
01398 
01399   for (octave_idx_type j = 0; j < npts*nsamples; j++)
01400     tmp_data[j] = tmp_data[j] / static_cast<float> (npts);
01401 
01402   return retval;
01403 }
01404 
01405 FloatComplexMatrix
01406 FloatComplexMatrix::fourier2d (void) const
01407 {
01408   FloatComplexMatrix retval;
01409 
01410   octave_idx_type nr = rows ();
01411   octave_idx_type nc = cols ();
01412 
01413   octave_idx_type npts, nsamples;
01414 
01415   if (nr == 1 || nc == 1)
01416     {
01417       npts = nr > nc ? nr : nc;
01418       nsamples = 1;
01419     }
01420   else
01421     {
01422       npts = nr;
01423       nsamples = nc;
01424     }
01425 
01426   octave_idx_type nn = 4*npts+15;
01427 
01428   Array<FloatComplex> wsave (dim_vector (nn, 1));
01429   FloatComplex *pwsave = wsave.fortran_vec ();
01430 
01431   retval = *this;
01432   FloatComplex *tmp_data = retval.fortran_vec ();
01433 
01434   F77_FUNC (cffti, CFFTI) (npts, pwsave);
01435 
01436   for (octave_idx_type j = 0; j < nsamples; j++)
01437     {
01438       octave_quit ();
01439 
01440       F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
01441     }
01442 
01443   npts = nc;
01444   nsamples = nr;
01445   nn = 4*npts+15;
01446 
01447   wsave.resize (dim_vector (nn, 1));
01448   pwsave = wsave.fortran_vec ();
01449 
01450   Array<FloatComplex> tmp (dim_vector (npts, 1));
01451   FloatComplex *prow = tmp.fortran_vec ();
01452 
01453   F77_FUNC (cffti, CFFTI) (npts, pwsave);
01454 
01455   for (octave_idx_type j = 0; j < nsamples; j++)
01456     {
01457       octave_quit ();
01458 
01459       for (octave_idx_type i = 0; i < npts; i++)
01460         prow[i] = tmp_data[i*nr + j];
01461 
01462       F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
01463 
01464       for (octave_idx_type i = 0; i < npts; i++)
01465         tmp_data[i*nr + j] = prow[i];
01466     }
01467 
01468   return retval;
01469 }
01470 
01471 FloatComplexMatrix
01472 FloatComplexMatrix::ifourier2d (void) const
01473 {
01474   FloatComplexMatrix retval;
01475 
01476   octave_idx_type nr = rows ();
01477   octave_idx_type nc = cols ();
01478 
01479   octave_idx_type npts, nsamples;
01480 
01481   if (nr == 1 || nc == 1)
01482     {
01483       npts = nr > nc ? nr : nc;
01484       nsamples = 1;
01485     }
01486   else
01487     {
01488       npts = nr;
01489       nsamples = nc;
01490     }
01491 
01492   octave_idx_type nn = 4*npts+15;
01493 
01494   Array<FloatComplex> wsave (dim_vector (nn, 1));
01495   FloatComplex *pwsave = wsave.fortran_vec ();
01496 
01497   retval = *this;
01498   FloatComplex *tmp_data = retval.fortran_vec ();
01499 
01500   F77_FUNC (cffti, CFFTI) (npts, pwsave);
01501 
01502   for (octave_idx_type j = 0; j < nsamples; j++)
01503     {
01504       octave_quit ();
01505 
01506       F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
01507     }
01508 
01509   for (octave_idx_type j = 0; j < npts*nsamples; j++)
01510     tmp_data[j] = tmp_data[j] / static_cast<float> (npts);
01511 
01512   npts = nc;
01513   nsamples = nr;
01514   nn = 4*npts+15;
01515 
01516   wsave.resize (dim_vector (nn, 1));
01517   pwsave = wsave.fortran_vec ();
01518 
01519   Array<FloatComplex> tmp (dim_vector (npts, 1));
01520   FloatComplex *prow = tmp.fortran_vec ();
01521 
01522   F77_FUNC (cffti, CFFTI) (npts, pwsave);
01523 
01524   for (octave_idx_type j = 0; j < nsamples; j++)
01525     {
01526       octave_quit ();
01527 
01528       for (octave_idx_type i = 0; i < npts; i++)
01529         prow[i] = tmp_data[i*nr + j];
01530 
01531       F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
01532 
01533       for (octave_idx_type i = 0; i < npts; i++)
01534         tmp_data[i*nr + j] = prow[i] / static_cast<float> (npts);
01535     }
01536 
01537   return retval;
01538 }
01539 
01540 #endif
01541 
01542 FloatComplexDET
01543 FloatComplexMatrix::determinant (void) const
01544 {
01545   octave_idx_type info;
01546   float rcon;
01547   return determinant (info, rcon, 0);
01548 }
01549 
01550 FloatComplexDET
01551 FloatComplexMatrix::determinant (octave_idx_type& info) const
01552 {
01553   float rcon;
01554   return determinant (info, rcon, 0);
01555 }
01556 
01557 FloatComplexDET
01558 FloatComplexMatrix::determinant (octave_idx_type& info, float& rcon, int calc_cond) const
01559 {
01560   MatrixType mattype (*this);
01561   return determinant (mattype, info, rcon, calc_cond);
01562 }
01563 
01564 FloatComplexDET
01565 FloatComplexMatrix::determinant (MatrixType& mattype,
01566                                  octave_idx_type& info, float& rcon, int calc_cond) const
01567 {
01568   FloatComplexDET retval (1.0);
01569 
01570   info = 0;
01571   rcon = 0.0;
01572 
01573   octave_idx_type nr = rows ();
01574   octave_idx_type nc = cols ();
01575 
01576   if (nr != nc)
01577     (*current_liboctave_error_handler) ("matrix must be square");
01578   else
01579     {
01580       volatile int typ = mattype.type ();
01581 
01582       // Even though the matrix is marked as singular (Rectangular), we may
01583       // still get a useful number from the LU factorization, because it always
01584       // completes.
01585 
01586       if (typ == MatrixType::Unknown)
01587         typ = mattype.type (*this);
01588       else if (typ == MatrixType::Rectangular)
01589         typ = MatrixType::Full;
01590 
01591       if (typ == MatrixType::Lower || typ == MatrixType::Upper)
01592         {
01593           for (octave_idx_type i = 0; i < nc; i++)
01594             retval *= elem (i,i);
01595         }
01596       else if (typ == MatrixType::Hermitian)
01597         {
01598           FloatComplexMatrix atmp = *this;
01599           FloatComplex *tmp_data = atmp.fortran_vec ();
01600 
01601           float anorm = 0;
01602           if (calc_cond) anorm = xnorm (*this, 1);
01603 
01604 
01605           char job = 'L';
01606           F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
01607                                      tmp_data, nr, info
01608                                      F77_CHAR_ARG_LEN (1)));
01609 
01610           if (info != 0)
01611             {
01612               rcon = 0.0;
01613               mattype.mark_as_unsymmetric ();
01614               typ = MatrixType::Full;
01615             }
01616           else
01617             {
01618               Array<FloatComplex> z (dim_vector (2 * nc, 1));
01619               FloatComplex *pz = z.fortran_vec ();
01620               Array<float> rz (dim_vector (nc, 1));
01621               float *prz = rz.fortran_vec ();
01622 
01623               F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
01624                                          nr, tmp_data, nr, anorm,
01625                                          rcon, pz, prz, info
01626                                          F77_CHAR_ARG_LEN (1)));
01627 
01628               if (info != 0)
01629                 rcon = 0.0;
01630 
01631               for (octave_idx_type i = 0; i < nc; i++)
01632                 retval *= atmp (i,i);
01633 
01634               retval = retval.square ();
01635             }
01636         }
01637       else if (typ != MatrixType::Full)
01638         (*current_liboctave_error_handler) ("det: invalid dense matrix type");
01639 
01640       if (typ == MatrixType::Full)
01641         {
01642           Array<octave_idx_type> ipvt (dim_vector (nr, 1));
01643           octave_idx_type *pipvt = ipvt.fortran_vec ();
01644 
01645           FloatComplexMatrix atmp = *this;
01646           FloatComplex *tmp_data = atmp.fortran_vec ();
01647 
01648           info = 0;
01649 
01650           // Calculate the norm of the matrix, for later use.
01651           float anorm = 0;
01652           if (calc_cond) anorm = xnorm (*this, 1);
01653 
01654           F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01655 
01656           // Throw-away extra info LAPACK gives so as to not change output.
01657           rcon = 0.0;
01658           if (info != 0)
01659             {
01660               info = -1;
01661               retval = FloatComplexDET ();
01662             }
01663           else
01664             {
01665               if (calc_cond)
01666                 {
01667                   // Now calc the condition number for non-singular matrix.
01668                   char job = '1';
01669                   Array<FloatComplex> z (dim_vector (2 * nc, 1));
01670                   FloatComplex *pz = z.fortran_vec ();
01671                   Array<float> rz (dim_vector (2 * nc, 1));
01672                   float *prz = rz.fortran_vec ();
01673 
01674                   F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01675                                              nc, tmp_data, nr, anorm,
01676                                              rcon, pz, prz, info
01677                                              F77_CHAR_ARG_LEN (1)));
01678                 }
01679 
01680               if (info != 0)
01681                 {
01682                   info = -1;
01683                   retval = FloatComplexDET ();
01684                 }
01685               else
01686                 {
01687                   for (octave_idx_type i = 0; i < nc; i++)
01688                     {
01689                       FloatComplex c = atmp(i,i);
01690                       retval *= (ipvt(i) != (i+1)) ? -c : c;
01691                     }
01692                 }
01693             }
01694         }
01695     }
01696 
01697   return retval;
01698 }
01699 
01700 float
01701 FloatComplexMatrix::rcond (void) const
01702 {
01703   MatrixType mattype (*this);
01704   return rcond (mattype);
01705 }
01706 
01707 float
01708 FloatComplexMatrix::rcond (MatrixType &mattype) const
01709 {
01710   float rcon;
01711   octave_idx_type nr = rows ();
01712   octave_idx_type nc = cols ();
01713 
01714   if (nr != nc)
01715     (*current_liboctave_error_handler) ("matrix must be square");
01716   else if (nr == 0 || nc == 0)
01717     rcon = octave_Inf;
01718   else
01719     {
01720       int typ = mattype.type ();
01721 
01722       if (typ == MatrixType::Unknown)
01723         typ = mattype.type (*this);
01724 
01725       // Only calculate the condition number for LU/Cholesky
01726       if (typ == MatrixType::Upper)
01727         {
01728           const FloatComplex *tmp_data = fortran_vec ();
01729           octave_idx_type info = 0;
01730           char norm = '1';
01731           char uplo = 'U';
01732           char dia = 'N';
01733 
01734           Array<FloatComplex> z (dim_vector (2 * nc, 1));
01735           FloatComplex *pz = z.fortran_vec ();
01736           Array<float> rz (dim_vector (nc, 1));
01737           float *prz = rz.fortran_vec ();
01738 
01739           F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
01740                                      F77_CONST_CHAR_ARG2 (&uplo, 1),
01741                                      F77_CONST_CHAR_ARG2 (&dia, 1),
01742                                      nr, tmp_data, nr, rcon,
01743                                      pz, prz, info
01744                                      F77_CHAR_ARG_LEN (1)
01745                                      F77_CHAR_ARG_LEN (1)
01746                                      F77_CHAR_ARG_LEN (1)));
01747 
01748           if (info != 0)
01749             rcon = 0;
01750         }
01751       else if  (typ == MatrixType::Permuted_Upper)
01752         (*current_liboctave_error_handler)
01753           ("permuted triangular matrix not implemented");
01754       else if (typ == MatrixType::Lower)
01755         {
01756           const FloatComplex *tmp_data = fortran_vec ();
01757           octave_idx_type info = 0;
01758           char norm = '1';
01759           char uplo = 'L';
01760           char dia = 'N';
01761 
01762           Array<FloatComplex> z (dim_vector (2 * nc, 1));
01763           FloatComplex *pz = z.fortran_vec ();
01764           Array<float> rz (dim_vector (nc, 1));
01765           float *prz = rz.fortran_vec ();
01766 
01767           F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
01768                                      F77_CONST_CHAR_ARG2 (&uplo, 1),
01769                                      F77_CONST_CHAR_ARG2 (&dia, 1),
01770                                      nr, tmp_data, nr, rcon,
01771                                      pz, prz, info
01772                                      F77_CHAR_ARG_LEN (1)
01773                                      F77_CHAR_ARG_LEN (1)
01774                                      F77_CHAR_ARG_LEN (1)));
01775 
01776           if (info != 0)
01777             rcon = 0.0;
01778         }
01779       else if (typ == MatrixType::Permuted_Lower)
01780         (*current_liboctave_error_handler)
01781           ("permuted triangular matrix not implemented");
01782       else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
01783         {
01784           float anorm = -1.0;
01785 
01786           if (typ == MatrixType::Hermitian)
01787             {
01788               octave_idx_type info = 0;
01789               char job = 'L';
01790 
01791               FloatComplexMatrix atmp = *this;
01792               FloatComplex *tmp_data = atmp.fortran_vec ();
01793 
01794               anorm = atmp.abs().sum().
01795                 row(static_cast<octave_idx_type>(0)).max();
01796 
01797               F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
01798                                          tmp_data, nr, info
01799                                          F77_CHAR_ARG_LEN (1)));
01800 
01801               if (info != 0)
01802                 {
01803                   rcon = 0.0;
01804 
01805                   mattype.mark_as_unsymmetric ();
01806                   typ = MatrixType::Full;
01807                 }
01808               else
01809                 {
01810                   Array<FloatComplex> z (dim_vector (2 * nc, 1));
01811                   FloatComplex *pz = z.fortran_vec ();
01812                   Array<float> rz (dim_vector (nc, 1));
01813                   float *prz = rz.fortran_vec ();
01814 
01815                   F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
01816                                              nr, tmp_data, nr, anorm,
01817                                              rcon, pz, prz, info
01818                                              F77_CHAR_ARG_LEN (1)));
01819 
01820                   if (info != 0)
01821                     rcon = 0.0;
01822                 }
01823             }
01824 
01825 
01826           if (typ == MatrixType::Full)
01827             {
01828               octave_idx_type info = 0;
01829 
01830               FloatComplexMatrix atmp = *this;
01831               FloatComplex *tmp_data = atmp.fortran_vec ();
01832 
01833               Array<octave_idx_type> ipvt (dim_vector (nr, 1));
01834               octave_idx_type *pipvt = ipvt.fortran_vec ();
01835 
01836               if(anorm < 0.)
01837                 anorm = atmp.abs().sum().
01838                   row(static_cast<octave_idx_type>(0)).max();
01839 
01840               Array<FloatComplex> z (dim_vector (2 * nc, 1));
01841               FloatComplex *pz = z.fortran_vec ();
01842               Array<float> rz (dim_vector (2 * nc, 1));
01843               float *prz = rz.fortran_vec ();
01844 
01845               F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01846 
01847               if (info != 0)
01848                 {
01849                   rcon = 0.0;
01850                   mattype.mark_as_rectangular ();
01851                 }
01852               else
01853                 {
01854                   char job = '1';
01855                   F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01856                                              nc, tmp_data, nr, anorm,
01857                                              rcon, pz, prz, info
01858                                              F77_CHAR_ARG_LEN (1)));
01859 
01860                   if (info != 0)
01861                     rcon = 0.0;
01862                 }
01863             }
01864         }
01865       else
01866         rcon = 0.0;
01867     }
01868 
01869   return rcon;
01870 }
01871 
01872 FloatComplexMatrix
01873 FloatComplexMatrix::utsolve (MatrixType &mattype, const FloatComplexMatrix& b,
01874                         octave_idx_type& info, float& rcon,
01875                         solve_singularity_handler sing_handler,
01876                         bool calc_cond, blas_trans_type transt) const
01877 {
01878   FloatComplexMatrix retval;
01879 
01880   octave_idx_type nr = rows ();
01881   octave_idx_type nc = cols ();
01882 
01883   if (nr != b.rows ())
01884     (*current_liboctave_error_handler)
01885       ("matrix dimension mismatch solution of linear equations");
01886   else if (nr == 0 || nc == 0 || b.cols () == 0)
01887     retval = FloatComplexMatrix (nc, b.cols (), FloatComplex (0.0, 0.0));
01888   else
01889     {
01890       volatile int typ = mattype.type ();
01891 
01892       if (typ == MatrixType::Permuted_Upper ||
01893           typ == MatrixType::Upper)
01894         {
01895           octave_idx_type b_nc = b.cols ();
01896           rcon = 1.;
01897           info = 0;
01898 
01899           if (typ == MatrixType::Permuted_Upper)
01900             {
01901               (*current_liboctave_error_handler)
01902                 ("permuted triangular matrix not implemented");
01903             }
01904           else
01905             {
01906               const FloatComplex *tmp_data = fortran_vec ();
01907 
01908               if (calc_cond)
01909                 {
01910                   char norm = '1';
01911                   char uplo = 'U';
01912                   char dia = 'N';
01913 
01914                   Array<FloatComplex> z (dim_vector (2 * nc, 1));
01915                   FloatComplex *pz = z.fortran_vec ();
01916                   Array<float> rz (dim_vector (nc, 1));
01917                   float *prz = rz.fortran_vec ();
01918 
01919                   F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
01920                                              F77_CONST_CHAR_ARG2 (&uplo, 1),
01921                                              F77_CONST_CHAR_ARG2 (&dia, 1),
01922                                              nr, tmp_data, nr, rcon,
01923                                              pz, prz, info
01924                                              F77_CHAR_ARG_LEN (1)
01925                                              F77_CHAR_ARG_LEN (1)
01926                                              F77_CHAR_ARG_LEN (1)));
01927 
01928                   if (info != 0)
01929                     info = -2;
01930 
01931                   volatile float rcond_plus_one = rcon + 1.0;
01932 
01933                   if (rcond_plus_one == 1.0 || xisnan (rcon))
01934                     {
01935                       info = -2;
01936 
01937                       if (sing_handler)
01938                         sing_handler (rcon);
01939                       else
01940                         (*current_liboctave_error_handler)
01941                           ("matrix singular to machine precision, rcond = %g",
01942                            rcon);
01943                     }
01944                 }
01945 
01946               if (info == 0)
01947                 {
01948                   retval = b;
01949                   FloatComplex *result = retval.fortran_vec ();
01950 
01951                   char uplo = 'U';
01952                   char trans = get_blas_char (transt);
01953                   char dia = 'N';
01954 
01955                   F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
01956                                              F77_CONST_CHAR_ARG2 (&trans, 1),
01957                                              F77_CONST_CHAR_ARG2 (&dia, 1),
01958                                              nr, b_nc, tmp_data, nr,
01959                                              result, nr, info
01960                                              F77_CHAR_ARG_LEN (1)
01961                                              F77_CHAR_ARG_LEN (1)
01962                                              F77_CHAR_ARG_LEN (1)));
01963                 }
01964             }
01965         }
01966       else
01967         (*current_liboctave_error_handler) ("incorrect matrix type");
01968     }
01969 
01970   return retval;
01971 }
01972 
01973 FloatComplexMatrix
01974 FloatComplexMatrix::ltsolve (MatrixType &mattype, const FloatComplexMatrix& b,
01975                         octave_idx_type& info, float& rcon,
01976                         solve_singularity_handler sing_handler,
01977                         bool calc_cond, blas_trans_type transt) const
01978 {
01979   FloatComplexMatrix retval;
01980 
01981   octave_idx_type nr = rows ();
01982   octave_idx_type nc = cols ();
01983 
01984   if (nr != b.rows ())
01985     (*current_liboctave_error_handler)
01986       ("matrix dimension mismatch solution of linear equations");
01987   else if (nr == 0 || nc == 0 || b.cols () == 0)
01988     retval = FloatComplexMatrix (nc, b.cols (), FloatComplex (0.0, 0.0));
01989   else
01990     {
01991       volatile int typ = mattype.type ();
01992 
01993       if (typ == MatrixType::Permuted_Lower ||
01994           typ == MatrixType::Lower)
01995         {
01996           octave_idx_type b_nc = b.cols ();
01997           rcon = 1.;
01998           info = 0;
01999 
02000           if (typ == MatrixType::Permuted_Lower)
02001             {
02002               (*current_liboctave_error_handler)
02003                 ("permuted triangular matrix not implemented");
02004             }
02005           else
02006             {
02007               const FloatComplex *tmp_data = fortran_vec ();
02008 
02009               if (calc_cond)
02010                 {
02011                   char norm = '1';
02012                   char uplo = 'L';
02013                   char dia = 'N';
02014 
02015                   Array<FloatComplex> z (dim_vector (2 * nc, 1));
02016                   FloatComplex *pz = z.fortran_vec ();
02017                   Array<float> rz (dim_vector (nc, 1));
02018                   float *prz = rz.fortran_vec ();
02019 
02020                   F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
02021                                              F77_CONST_CHAR_ARG2 (&uplo, 1),
02022                                              F77_CONST_CHAR_ARG2 (&dia, 1),
02023                                              nr, tmp_data, nr, rcon,
02024                                              pz, prz, info
02025                                              F77_CHAR_ARG_LEN (1)
02026                                              F77_CHAR_ARG_LEN (1)
02027                                              F77_CHAR_ARG_LEN (1)));
02028 
02029                   if (info != 0)
02030                     info = -2;
02031 
02032                   volatile float rcond_plus_one = rcon + 1.0;
02033 
02034                   if (rcond_plus_one == 1.0 || xisnan (rcon))
02035                     {
02036                       info = -2;
02037 
02038                       if (sing_handler)
02039                         sing_handler (rcon);
02040                       else
02041                         (*current_liboctave_error_handler)
02042                           ("matrix singular to machine precision, rcond = %g",
02043                            rcon);
02044                     }
02045                 }
02046 
02047               if (info == 0)
02048                 {
02049                   retval = b;
02050                   FloatComplex *result = retval.fortran_vec ();
02051 
02052                   char uplo = 'L';
02053                   char trans = get_blas_char (transt);
02054                   char dia = 'N';
02055 
02056                   F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
02057                                              F77_CONST_CHAR_ARG2 (&trans, 1),
02058                                              F77_CONST_CHAR_ARG2 (&dia, 1),
02059                                              nr, b_nc, tmp_data, nr,
02060                                              result, nr, info
02061                                              F77_CHAR_ARG_LEN (1)
02062                                              F77_CHAR_ARG_LEN (1)
02063                                              F77_CHAR_ARG_LEN (1)));
02064                 }
02065             }
02066         }
02067       else
02068         (*current_liboctave_error_handler) ("incorrect matrix type");
02069     }
02070 
02071   return retval;
02072 }
02073 
02074 FloatComplexMatrix
02075 FloatComplexMatrix::fsolve (MatrixType &mattype, const FloatComplexMatrix& b,
02076                        octave_idx_type& info, float& rcon,
02077                        solve_singularity_handler sing_handler,
02078                        bool calc_cond) const
02079 {
02080   FloatComplexMatrix retval;
02081 
02082   octave_idx_type nr = rows ();
02083   octave_idx_type nc = cols ();
02084 
02085 
02086   if (nr != nc || nr != b.rows ())
02087     (*current_liboctave_error_handler)
02088       ("matrix dimension mismatch solution of linear equations");
02089   else if (nr == 0 || b.cols () == 0)
02090     retval = FloatComplexMatrix (nc, b.cols (), FloatComplex (0.0, 0.0));
02091   else
02092     {
02093       volatile int typ = mattype.type ();
02094 
02095      // Calculate the norm of the matrix, for later use.
02096       float anorm = -1.;
02097 
02098       if (typ == MatrixType::Hermitian)
02099         {
02100           info = 0;
02101           char job = 'L';
02102 
02103           FloatComplexMatrix atmp = *this;
02104           FloatComplex *tmp_data = atmp.fortran_vec ();
02105 
02106           anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
02107 
02108           F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
02109                                      tmp_data, nr, info
02110                                      F77_CHAR_ARG_LEN (1)));
02111 
02112           // Throw-away extra info LAPACK gives so as to not change output.
02113           rcon = 0.0;
02114           if (info != 0)
02115             {
02116               info = -2;
02117 
02118               mattype.mark_as_unsymmetric ();
02119               typ = MatrixType::Full;
02120             }
02121           else
02122             {
02123               if (calc_cond)
02124                 {
02125                   Array<FloatComplex> z (dim_vector (2 * nc, 1));
02126                   FloatComplex *pz = z.fortran_vec ();
02127                   Array<float> rz (dim_vector (nc, 1));
02128                   float *prz = rz.fortran_vec ();
02129 
02130                   F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
02131                                              nr, tmp_data, nr, anorm,
02132                                              rcon, pz, prz, info
02133                                              F77_CHAR_ARG_LEN (1)));
02134 
02135                   if (info != 0)
02136                     info = -2;
02137 
02138                   volatile float rcond_plus_one = rcon + 1.0;
02139 
02140                   if (rcond_plus_one == 1.0 || xisnan (rcon))
02141                     {
02142                       info = -2;
02143 
02144                       if (sing_handler)
02145                         sing_handler (rcon);
02146                       else
02147                         (*current_liboctave_error_handler)
02148                           ("matrix singular to machine precision, rcond = %g",
02149                            rcon);
02150                     }
02151                 }
02152 
02153               if (info == 0)
02154                 {
02155                   retval = b;
02156                   FloatComplex *result = retval.fortran_vec ();
02157 
02158                   octave_idx_type b_nc = b.cols ();
02159 
02160                   F77_XFCN (cpotrs, CPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
02161                                              nr, b_nc, tmp_data, nr,
02162                                              result, b.rows(), info
02163                                              F77_CHAR_ARG_LEN (1)));
02164                 }
02165               else
02166                 {
02167                   mattype.mark_as_unsymmetric ();
02168                   typ = MatrixType::Full;
02169                 }
02170             }
02171         }
02172 
02173       if (typ == MatrixType::Full)
02174         {
02175           info = 0;
02176 
02177           Array<octave_idx_type> ipvt (dim_vector (nr, 1));
02178           octave_idx_type *pipvt = ipvt.fortran_vec ();
02179 
02180           FloatComplexMatrix atmp = *this;
02181           FloatComplex *tmp_data = atmp.fortran_vec ();
02182 
02183           Array<FloatComplex> z (dim_vector (2 * nc, 1));
02184           FloatComplex *pz = z.fortran_vec ();
02185           Array<float> rz (dim_vector (2 * nc, 1));
02186           float *prz = rz.fortran_vec ();
02187 
02188           // Calculate the norm of the matrix, for later use.
02189           if (anorm < 0.)
02190             anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
02191 
02192           F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
02193 
02194           // Throw-away extra info LAPACK gives so as to not change output.
02195           rcon = 0.0;
02196           if (info != 0)
02197             {
02198               info = -2;
02199 
02200               if (sing_handler)
02201                 sing_handler (rcon);
02202               else
02203                 (*current_liboctave_error_handler)
02204                   ("matrix singular to machine precision");
02205 
02206               mattype.mark_as_rectangular ();
02207             }
02208           else
02209             {
02210               if (calc_cond)
02211                 {
02212                   // Now calculate the condition number for
02213                   // non-singular matrix.
02214                   char job = '1';
02215                   F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
02216                                              nc, tmp_data, nr, anorm,
02217                                              rcon, pz, prz, info
02218                                              F77_CHAR_ARG_LEN (1)));
02219 
02220                   if (info != 0)
02221                     info = -2;
02222 
02223                   volatile float rcond_plus_one = rcon + 1.0;
02224 
02225                   if (rcond_plus_one == 1.0 || xisnan (rcon))
02226                     {
02227                       info = -2;
02228 
02229                       if (sing_handler)
02230                         sing_handler (rcon);
02231                       else
02232                         (*current_liboctave_error_handler)
02233                           ("matrix singular to machine precision, rcond = %g",
02234                            rcon);
02235                     }
02236                 }
02237 
02238               if (info == 0)
02239                 {
02240                   retval = b;
02241                   FloatComplex *result = retval.fortran_vec ();
02242 
02243                   octave_idx_type b_nc = b.cols ();
02244 
02245                   char job = 'N';
02246                   F77_XFCN (cgetrs, CGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
02247                                              nr, b_nc, tmp_data, nr,
02248                                              pipvt, result, b.rows(), info
02249                                              F77_CHAR_ARG_LEN (1)));
02250                 }
02251               else
02252                 mattype.mark_as_rectangular ();
02253             }
02254         }
02255     }
02256 
02257   return retval;
02258 }
02259 
02260 FloatComplexMatrix
02261 FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b) const
02262 {
02263   octave_idx_type info;
02264   float rcon;
02265   return solve (typ, b, info, rcon, 0);
02266 }
02267 
02268 FloatComplexMatrix
02269 FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b,
02270                       octave_idx_type& info) const
02271 {
02272   float rcon;
02273   return solve (typ, b, info, rcon, 0);
02274 }
02275 
02276 FloatComplexMatrix
02277 FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info,
02278                       float& rcon) const
02279 {
02280   return solve (typ, b, info, rcon, 0);
02281 }
02282 
02283 FloatComplexMatrix
02284 FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info,
02285                       float& rcon, solve_singularity_handler sing_handler,
02286                       bool singular_fallback, blas_trans_type transt) const
02287 {
02288   FloatComplexMatrix tmp (b);
02289   return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt);
02290 }
02291 
02292 FloatComplexMatrix
02293 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b) const
02294 {
02295   octave_idx_type info;
02296   float rcon;
02297   return solve (typ, b, info, rcon, 0);
02298 }
02299 
02300 FloatComplexMatrix
02301 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b,
02302                       octave_idx_type& info) const
02303 {
02304   float rcon;
02305   return solve (typ, b, info, rcon, 0);
02306 }
02307 
02308 FloatComplexMatrix
02309 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b,
02310                       octave_idx_type& info, float& rcon) const
02311 {
02312   return solve (typ, b, info, rcon, 0);
02313 }
02314 
02315 FloatComplexMatrix
02316 FloatComplexMatrix::solve (MatrixType &mattype, const FloatComplexMatrix& b,
02317                       octave_idx_type& info, float& rcon,
02318                       solve_singularity_handler sing_handler,
02319                       bool singular_fallback, blas_trans_type transt) const
02320 {
02321   FloatComplexMatrix retval;
02322   int typ = mattype.type ();
02323 
02324   if (typ == MatrixType::Unknown)
02325     typ = mattype.type (*this);
02326 
02327   // Only calculate the condition number for LU/Cholesky
02328   if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
02329     retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt);
02330   else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
02331     retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt);
02332   else if (transt == blas_trans)
02333     return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback);
02334   else if (transt == blas_conj_trans)
02335     retval = hermitian ().solve (mattype, b, info, rcon, sing_handler, singular_fallback);
02336   else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
02337     retval = fsolve (mattype, b, info, rcon, sing_handler, true);
02338   else if (typ != MatrixType::Rectangular)
02339     {
02340       (*current_liboctave_error_handler) ("unknown matrix type");
02341       return FloatComplexMatrix ();
02342     }
02343 
02344   // Rectangular or one of the above solvers flags a singular matrix
02345   if (singular_fallback && mattype.type () == MatrixType::Rectangular)
02346     {
02347       octave_idx_type rank;
02348       retval = lssolve (b, info, rank, rcon);
02349     }
02350 
02351   return retval;
02352 }
02353 
02354 FloatComplexColumnVector
02355 FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b) const
02356 {
02357   octave_idx_type info;
02358   float rcon;
02359   return solve (typ, FloatComplexColumnVector (b), info, rcon, 0);
02360 }
02361 
02362 FloatComplexColumnVector
02363 FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b,
02364                       octave_idx_type& info) const
02365 {
02366   float rcon;
02367   return solve (typ, FloatComplexColumnVector (b), info, rcon, 0);
02368 }
02369 
02370 FloatComplexColumnVector
02371 FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b,
02372                       octave_idx_type& info, float& rcon) const
02373 {
02374   return solve (typ, FloatComplexColumnVector (b), info, rcon, 0);
02375 }
02376 
02377 FloatComplexColumnVector
02378 FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b,
02379                       octave_idx_type& info, float& rcon,
02380                       solve_singularity_handler sing_handler, blas_trans_type transt) const
02381 {
02382   return solve (typ, FloatComplexColumnVector (b), info, rcon, sing_handler, transt);
02383 }
02384 
02385 FloatComplexColumnVector
02386 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b) const
02387 {
02388   octave_idx_type info;
02389   float rcon;
02390   return solve (typ, b, info, rcon, 0);
02391 }
02392 
02393 FloatComplexColumnVector
02394 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
02395                       octave_idx_type& info) const
02396 {
02397   float rcon;
02398   return solve (typ, b, info, rcon, 0);
02399 }
02400 
02401 FloatComplexColumnVector
02402 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
02403                       octave_idx_type& info, float& rcon) const
02404 {
02405   return solve (typ, b, info, rcon, 0);
02406 }
02407 
02408 FloatComplexColumnVector
02409 FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
02410                       octave_idx_type& info, float& rcon,
02411                       solve_singularity_handler sing_handler, blas_trans_type transt) const
02412 {
02413 
02414   FloatComplexMatrix tmp (b);
02415   tmp = solve (typ, tmp, info, rcon, sing_handler, true, transt);
02416   return tmp.column(static_cast<octave_idx_type> (0));
02417 }
02418 
02419 FloatComplexMatrix
02420 FloatComplexMatrix::solve (const FloatMatrix& b) const
02421 {
02422   octave_idx_type info;
02423   float rcon;
02424   return solve (b, info, rcon, 0);
02425 }
02426 
02427 FloatComplexMatrix
02428 FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info) const
02429 {
02430   float rcon;
02431   return solve (b, info, rcon, 0);
02432 }
02433 
02434 FloatComplexMatrix
02435 FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const
02436 {
02437   return solve (b, info, rcon, 0);
02438 }
02439 
02440 FloatComplexMatrix
02441 FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon,
02442                       solve_singularity_handler sing_handler, blas_trans_type transt) const
02443 {
02444   FloatComplexMatrix tmp (b);
02445   return solve (tmp, info, rcon, sing_handler, transt);
02446 }
02447 
02448 FloatComplexMatrix
02449 FloatComplexMatrix::solve (const FloatComplexMatrix& b) const
02450 {
02451   octave_idx_type info;
02452   float rcon;
02453   return solve (b, info, rcon, 0);
02454 }
02455 
02456 FloatComplexMatrix
02457 FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info) const
02458 {
02459   float rcon;
02460   return solve (b, info, rcon, 0);
02461 }
02462 
02463 FloatComplexMatrix
02464 FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const
02465 {
02466   return solve (b, info, rcon, 0);
02467 }
02468 
02469 FloatComplexMatrix
02470 FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon,
02471                       solve_singularity_handler sing_handler, blas_trans_type transt) const
02472 {
02473   MatrixType mattype (*this);
02474   return solve (mattype, b, info, rcon, sing_handler, true, transt);
02475 }
02476 
02477 FloatComplexColumnVector
02478 FloatComplexMatrix::solve (const FloatColumnVector& b) const
02479 {
02480   octave_idx_type info;
02481   float rcon;
02482   return solve (FloatComplexColumnVector (b), info, rcon, 0);
02483 }
02484 
02485 FloatComplexColumnVector
02486 FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info) const
02487 {
02488   float rcon;
02489   return solve (FloatComplexColumnVector (b), info, rcon, 0);
02490 }
02491 
02492 FloatComplexColumnVector
02493 FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info,
02494                       float& rcon) const
02495 {
02496   return solve (FloatComplexColumnVector (b), info, rcon, 0);
02497 }
02498 
02499 FloatComplexColumnVector
02500 FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info,
02501                       float& rcon,
02502                       solve_singularity_handler sing_handler, blas_trans_type transt) const
02503 {
02504   return solve (FloatComplexColumnVector (b), info, rcon, sing_handler, transt);
02505 }
02506 
02507 FloatComplexColumnVector
02508 FloatComplexMatrix::solve (const FloatComplexColumnVector& b) const
02509 {
02510   octave_idx_type info;
02511   float rcon;
02512   return solve (b, info, rcon, 0);
02513 }
02514 
02515 FloatComplexColumnVector
02516 FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info) const
02517 {
02518   float rcon;
02519   return solve (b, info, rcon, 0);
02520 }
02521 
02522 FloatComplexColumnVector
02523 FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info,
02524                       float& rcon) const
02525 {
02526   return solve (b, info, rcon, 0);
02527 }
02528 
02529 FloatComplexColumnVector
02530 FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info,
02531                       float& rcon,
02532                       solve_singularity_handler sing_handler, blas_trans_type transt) const
02533 {
02534   MatrixType mattype (*this);
02535   return solve (mattype, b, info, rcon, sing_handler, transt);
02536 }
02537 
02538 FloatComplexMatrix
02539 FloatComplexMatrix::lssolve (const FloatMatrix& b) const
02540 {
02541   octave_idx_type info;
02542   octave_idx_type rank;
02543   float rcon;
02544   return lssolve (FloatComplexMatrix (b), info, rank, rcon);
02545 }
02546 
02547 FloatComplexMatrix
02548 FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info) const
02549 {
02550   octave_idx_type rank;
02551   float rcon;
02552   return lssolve (FloatComplexMatrix (b), info, rank, rcon);
02553 }
02554 
02555 FloatComplexMatrix
02556 FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info,
02557                         octave_idx_type& rank) const
02558 {
02559   float rcon;
02560   return lssolve (FloatComplexMatrix (b), info, rank, rcon);
02561 }
02562 
02563 FloatComplexMatrix
02564 FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info,
02565                         octave_idx_type& rank, float& rcon) const
02566 {
02567   return lssolve (FloatComplexMatrix (b), info, rank, rcon);
02568 }
02569 
02570 FloatComplexMatrix
02571 FloatComplexMatrix::lssolve (const FloatComplexMatrix& b) const
02572 {
02573   octave_idx_type info;
02574   octave_idx_type rank;
02575   float rcon;
02576   return lssolve (b, info, rank, rcon);
02577 }
02578 
02579 FloatComplexMatrix
02580 FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info) const
02581 {
02582   octave_idx_type rank;
02583   float rcon;
02584   return lssolve (b, info, rank, rcon);
02585 }
02586 
02587 FloatComplexMatrix
02588 FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info,
02589                         octave_idx_type& rank) const
02590 {
02591   float rcon;
02592   return lssolve (b, info, rank, rcon);
02593 }
02594 
02595 FloatComplexMatrix
02596 FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info,
02597                         octave_idx_type& rank, float& rcon) const
02598 {
02599   FloatComplexMatrix retval;
02600 
02601   octave_idx_type nrhs = b.cols ();
02602 
02603   octave_idx_type m = rows ();
02604   octave_idx_type n = cols ();
02605 
02606   if (m != b.rows ())
02607     (*current_liboctave_error_handler)
02608       ("matrix dimension mismatch solution of linear equations");
02609   else if (m== 0 || n == 0 || b.cols () == 0)
02610     retval = FloatComplexMatrix (n, b.cols (), FloatComplex (0.0, 0.0));
02611   else
02612     {
02613       volatile octave_idx_type minmn = (m < n ? m : n);
02614       octave_idx_type maxmn = m > n ? m : n;
02615       rcon = -1.0;
02616 
02617       if (m != n)
02618         {
02619           retval = FloatComplexMatrix (maxmn, nrhs);
02620 
02621           for (octave_idx_type j = 0; j < nrhs; j++)
02622             for (octave_idx_type i = 0; i < m; i++)
02623               retval.elem (i, j) = b.elem (i, j);
02624         }
02625       else
02626         retval = b;
02627 
02628       FloatComplexMatrix atmp = *this;
02629       FloatComplex *tmp_data = atmp.fortran_vec ();
02630 
02631       FloatComplex *pretval = retval.fortran_vec ();
02632       Array<float> s (dim_vector (minmn, 1));
02633       float *ps = s.fortran_vec ();
02634 
02635       // Ask ZGELSD what the dimension of WORK should be.
02636       octave_idx_type lwork = -1;
02637 
02638       Array<FloatComplex> work (dim_vector (1, 1));
02639 
02640       octave_idx_type smlsiz;
02641       F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("CGELSD", 6),
02642                                    F77_CONST_CHAR_ARG2 (" ", 1),
02643                                    0, 0, 0, 0, smlsiz
02644                                    F77_CHAR_ARG_LEN (6)
02645                                    F77_CHAR_ARG_LEN (1));
02646 
02647       octave_idx_type mnthr;
02648       F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("CGELSD", 6),
02649                                    F77_CONST_CHAR_ARG2 (" ", 1),
02650                                    m, n, nrhs, -1, mnthr
02651                                    F77_CHAR_ARG_LEN (6)
02652                                    F77_CHAR_ARG_LEN (1));
02653 
02654       // We compute the size of rwork and iwork because ZGELSD in
02655       // older versions of LAPACK does not return them on a query
02656       // call.
02657       float dminmn = static_cast<float> (minmn);
02658       float dsmlsizp1 = static_cast<float> (smlsiz+1);
02659 #if defined (HAVE_LOG2)
02660       float tmp = log2 (dminmn / dsmlsizp1);
02661 #else
02662       float tmp = log (dminmn / dsmlsizp1) / log (2.0);
02663 #endif
02664       octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1;
02665       if (nlvl < 0)
02666         nlvl = 0;
02667 
02668       octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
02669         + 3*smlsiz*nrhs + std::max ((smlsiz+1)*(smlsiz+1),
02670                                     n*(1+nrhs) + 2*nrhs);
02671       if (lrwork < 1)
02672         lrwork = 1;
02673       Array<float> rwork (dim_vector (lrwork, 1));
02674       float *prwork = rwork.fortran_vec ();
02675 
02676       octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
02677       if (liwork < 1)
02678         liwork = 1;
02679       Array<octave_idx_type> iwork (dim_vector (liwork, 1));
02680       octave_idx_type* piwork = iwork.fortran_vec ();
02681 
02682       F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
02683                                  ps, rcon, rank, work.fortran_vec (),
02684                                  lwork, prwork, piwork, info));
02685 
02686       // The workspace query is broken in at least LAPACK 3.0.0
02687       // through 3.1.1 when n >= mnthr.  The obtuse formula below
02688       // should provide sufficient workspace for ZGELSD to operate
02689       // efficiently.
02690       if (n > m && n >= mnthr)
02691         {
02692           octave_idx_type addend = m;
02693 
02694           if (2*m-4 > addend)
02695             addend = 2*m-4;
02696 
02697           if (nrhs > addend)
02698             addend = nrhs;
02699 
02700           if (n-3*m > addend)
02701             addend = n-3*m;
02702 
02703           const octave_idx_type lworkaround = 4*m + m*m + addend;
02704 
02705           if (std::real (work(0)) < lworkaround)
02706             work(0) = lworkaround;
02707         }
02708       else if (m >= n)
02709         {
02710           octave_idx_type lworkaround = 2*m + m*nrhs;
02711 
02712           if (std::real (work(0)) < lworkaround)
02713             work(0) = lworkaround;
02714         }
02715 
02716       lwork = static_cast<octave_idx_type> (std::real (work(0)));
02717       work.resize (dim_vector (lwork, 1));
02718 
02719       F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval,
02720                                  maxmn, ps, rcon, rank,
02721                                  work.fortran_vec (), lwork,
02722                                  prwork, piwork, info));
02723 
02724       if (s.elem (0) == 0.0)
02725         rcon = 0.0;
02726       else
02727         rcon = s.elem (minmn - 1) / s.elem (0);
02728 
02729       retval.resize (n, nrhs);
02730     }
02731 
02732   return retval;
02733 }
02734 
02735 FloatComplexColumnVector
02736 FloatComplexMatrix::lssolve (const FloatColumnVector& b) const
02737 {
02738   octave_idx_type info;
02739   octave_idx_type rank;
02740   float rcon;
02741   return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
02742 }
02743 
02744 FloatComplexColumnVector
02745 FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info) const
02746 {
02747   octave_idx_type rank;
02748   float rcon;
02749   return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
02750 }
02751 
02752 FloatComplexColumnVector
02753 FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info,
02754                         octave_idx_type& rank) const
02755 {
02756   float rcon;
02757   return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
02758 }
02759 
02760 FloatComplexColumnVector
02761 FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info,
02762                         octave_idx_type& rank, float& rcon) const
02763 {
02764   return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
02765 }
02766 
02767 FloatComplexColumnVector
02768 FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b) const
02769 {
02770   octave_idx_type info;
02771   octave_idx_type rank;
02772   float rcon;
02773   return lssolve (b, info, rank, rcon);
02774 }
02775 
02776 FloatComplexColumnVector
02777 FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info) const
02778 {
02779   octave_idx_type rank;
02780   float rcon;
02781   return lssolve (b, info, rank, rcon);
02782 }
02783 
02784 FloatComplexColumnVector
02785 FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info,
02786                         octave_idx_type& rank) const
02787 {
02788   float rcon;
02789   return lssolve (b, info, rank, rcon);
02790 
02791 }
02792 
02793 FloatComplexColumnVector
02794 FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info,
02795                         octave_idx_type& rank, float& rcon) const
02796 {
02797   FloatComplexColumnVector retval;
02798 
02799   octave_idx_type nrhs = 1;
02800 
02801   octave_idx_type m = rows ();
02802   octave_idx_type n = cols ();
02803 
02804   if (m != b.length ())
02805     (*current_liboctave_error_handler)
02806       ("matrix dimension mismatch solution of linear equations");
02807   else if (m == 0 || n == 0 || b.cols () == 0)
02808     retval = FloatComplexColumnVector (n, FloatComplex (0.0, 0.0));
02809   else
02810     {
02811       volatile octave_idx_type minmn = (m < n ? m : n);
02812       octave_idx_type maxmn = m > n ? m : n;
02813       rcon = -1.0;
02814 
02815       if (m != n)
02816         {
02817           retval = FloatComplexColumnVector (maxmn);
02818 
02819           for (octave_idx_type i = 0; i < m; i++)
02820             retval.elem (i) = b.elem (i);
02821         }
02822       else
02823         retval = b;
02824 
02825       FloatComplexMatrix atmp = *this;
02826       FloatComplex *tmp_data = atmp.fortran_vec ();
02827 
02828       FloatComplex *pretval = retval.fortran_vec ();
02829       Array<float> s (dim_vector (minmn, 1));
02830       float *ps = s.fortran_vec ();
02831 
02832       // Ask ZGELSD what the dimension of WORK should be.
02833       octave_idx_type lwork = -1;
02834 
02835       Array<FloatComplex> work (dim_vector (1, 1));
02836 
02837       octave_idx_type smlsiz;
02838       F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("CGELSD", 6),
02839                                    F77_CONST_CHAR_ARG2 (" ", 1),
02840                                    0, 0, 0, 0, smlsiz
02841                                    F77_CHAR_ARG_LEN (6)
02842                                    F77_CHAR_ARG_LEN (1));
02843 
02844       // We compute the size of rwork and iwork because ZGELSD in
02845       // older versions of LAPACK does not return them on a query
02846       // call.
02847       float dminmn = static_cast<float> (minmn);
02848       float dsmlsizp1 = static_cast<float> (smlsiz+1);
02849 #if defined (HAVE_LOG2)
02850       float tmp = log2 (dminmn / dsmlsizp1);
02851 #else
02852       float tmp = log (dminmn / dsmlsizp1) / log (2.0);
02853 #endif
02854       octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1;
02855       if (nlvl < 0)
02856         nlvl = 0;
02857 
02858       octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
02859         + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
02860       if (lrwork < 1)
02861         lrwork = 1;
02862       Array<float> rwork (dim_vector (lrwork, 1));
02863       float *prwork = rwork.fortran_vec ();
02864 
02865       octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
02866       if (liwork < 1)
02867         liwork = 1;
02868       Array<octave_idx_type> iwork (dim_vector (liwork, 1));
02869       octave_idx_type* piwork = iwork.fortran_vec ();
02870 
02871       F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
02872                                  ps, rcon, rank, work.fortran_vec (),
02873                                  lwork, prwork, piwork, info));
02874 
02875       lwork = static_cast<octave_idx_type> (std::real (work(0)));
02876       work.resize (dim_vector (lwork, 1));
02877       rwork.resize (dim_vector (static_cast<octave_idx_type> (rwork(0)), 1));
02878       iwork.resize (dim_vector (iwork(0), 1));
02879 
02880       F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval,
02881                                  maxmn, ps, rcon, rank,
02882                                  work.fortran_vec (), lwork,
02883                                  prwork, piwork, info));
02884 
02885       if (rank < minmn)
02886         {
02887           if (s.elem (0) == 0.0)
02888             rcon = 0.0;
02889           else
02890             rcon = s.elem (minmn - 1) / s.elem (0);
02891 
02892           retval.resize (n, nrhs);
02893         }
02894     }
02895 
02896   return retval;
02897 }
02898 
02899 // column vector by row vector -> matrix operations
02900 
02901 FloatComplexMatrix
02902 operator * (const FloatColumnVector& v, const FloatComplexRowVector& a)
02903 {
02904   FloatComplexColumnVector tmp (v);
02905   return tmp * a;
02906 }
02907 
02908 FloatComplexMatrix
02909 operator * (const FloatComplexColumnVector& a, const FloatRowVector& b)
02910 {
02911   FloatComplexRowVector tmp (b);
02912   return a * tmp;
02913 }
02914 
02915 FloatComplexMatrix
02916 operator * (const FloatComplexColumnVector& v, const FloatComplexRowVector& a)
02917 {
02918   FloatComplexMatrix retval;
02919 
02920   octave_idx_type len = v.length ();
02921 
02922   if (len != 0)
02923     {
02924       octave_idx_type a_len = a.length ();
02925 
02926       retval = FloatComplexMatrix (len, a_len);
02927       FloatComplex *c = retval.fortran_vec ();
02928 
02929       F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
02930                                F77_CONST_CHAR_ARG2 ("N", 1),
02931                                len, a_len, 1, 1.0, v.data (), len,
02932                                a.data (), 1, 0.0, c, len
02933                                F77_CHAR_ARG_LEN (1)
02934                                F77_CHAR_ARG_LEN (1)));
02935     }
02936 
02937   return retval;
02938 }
02939 
02940 // matrix by diagonal matrix -> matrix operations
02941 
02942 FloatComplexMatrix&
02943 FloatComplexMatrix::operator += (const FloatDiagMatrix& a)
02944 {
02945   octave_idx_type nr = rows ();
02946   octave_idx_type nc = cols ();
02947 
02948   octave_idx_type a_nr = rows ();
02949   octave_idx_type a_nc = cols ();
02950 
02951   if (nr != a_nr || nc != a_nc)
02952     {
02953       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
02954       return *this;
02955     }
02956 
02957   for (octave_idx_type i = 0; i < a.length (); i++)
02958     elem (i, i) += a.elem (i, i);
02959 
02960   return *this;
02961 }
02962 
02963 FloatComplexMatrix&
02964 FloatComplexMatrix::operator -= (const FloatDiagMatrix& a)
02965 {
02966   octave_idx_type nr = rows ();
02967   octave_idx_type nc = cols ();
02968 
02969   octave_idx_type a_nr = rows ();
02970   octave_idx_type a_nc = cols ();
02971 
02972   if (nr != a_nr || nc != a_nc)
02973     {
02974       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
02975       return *this;
02976     }
02977 
02978   for (octave_idx_type i = 0; i < a.length (); i++)
02979     elem (i, i) -= a.elem (i, i);
02980 
02981   return *this;
02982 }
02983 
02984 FloatComplexMatrix&
02985 FloatComplexMatrix::operator += (const FloatComplexDiagMatrix& a)
02986 {
02987   octave_idx_type nr = rows ();
02988   octave_idx_type nc = cols ();
02989 
02990   octave_idx_type a_nr = rows ();
02991   octave_idx_type a_nc = cols ();
02992 
02993   if (nr != a_nr || nc != a_nc)
02994     {
02995       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
02996       return *this;
02997     }
02998 
02999   for (octave_idx_type i = 0; i < a.length (); i++)
03000     elem (i, i) += a.elem (i, i);
03001 
03002   return *this;
03003 }
03004 
03005 FloatComplexMatrix&
03006 FloatComplexMatrix::operator -= (const FloatComplexDiagMatrix& a)
03007 {
03008   octave_idx_type nr = rows ();
03009   octave_idx_type nc = cols ();
03010 
03011   octave_idx_type a_nr = rows ();
03012   octave_idx_type a_nc = cols ();
03013 
03014   if (nr != a_nr || nc != a_nc)
03015     {
03016       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
03017       return *this;
03018     }
03019 
03020   for (octave_idx_type i = 0; i < a.length (); i++)
03021     elem (i, i) -= a.elem (i, i);
03022 
03023   return *this;
03024 }
03025 
03026 // matrix by matrix -> matrix operations
03027 
03028 FloatComplexMatrix&
03029 FloatComplexMatrix::operator += (const FloatMatrix& a)
03030 {
03031   octave_idx_type nr = rows ();
03032   octave_idx_type nc = cols ();
03033 
03034   octave_idx_type a_nr = a.rows ();
03035   octave_idx_type a_nc = a.cols ();
03036 
03037   if (nr != a_nr || nc != a_nc)
03038     {
03039       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
03040       return *this;
03041     }
03042 
03043   if (nr == 0 || nc == 0)
03044     return *this;
03045 
03046   FloatComplex *d = fortran_vec (); // Ensures only one reference to my privates!
03047 
03048   mx_inline_add2 (length (), d, a.data ());
03049   return *this;
03050 }
03051 
03052 FloatComplexMatrix&
03053 FloatComplexMatrix::operator -= (const FloatMatrix& a)
03054 {
03055   octave_idx_type nr = rows ();
03056   octave_idx_type nc = cols ();
03057 
03058   octave_idx_type a_nr = a.rows ();
03059   octave_idx_type a_nc = a.cols ();
03060 
03061   if (nr != a_nr || nc != a_nc)
03062     {
03063       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
03064       return *this;
03065     }
03066 
03067   if (nr == 0 || nc == 0)
03068     return *this;
03069 
03070   FloatComplex *d = fortran_vec (); // Ensures only one reference to my privates!
03071 
03072   mx_inline_sub2 (length (), d, a.data ());
03073   return *this;
03074 }
03075 
03076 // unary operations
03077 
03078 boolMatrix
03079 FloatComplexMatrix::operator ! (void) const
03080 {
03081   if (any_element_is_nan ())
03082     gripe_nan_to_logical_conversion ();
03083 
03084   return do_mx_unary_op<bool, FloatComplex> (*this, mx_inline_not);
03085 }
03086 
03087 // other operations
03088 
03089 bool
03090 FloatComplexMatrix::any_element_is_nan (void) const
03091 {
03092   return do_mx_check<FloatComplex> (*this, mx_inline_any_nan);
03093 }
03094 
03095 bool
03096 FloatComplexMatrix::any_element_is_inf_or_nan (void) const
03097 {
03098   return ! do_mx_check<FloatComplex> (*this, mx_inline_all_finite);
03099 }
03100 
03101 // Return true if no elements have imaginary components.
03102 
03103 bool
03104 FloatComplexMatrix::all_elements_are_real (void) const
03105 {
03106   return do_mx_check<FloatComplex> (*this, mx_inline_all_real);
03107 }
03108 
03109 // Return nonzero if any element of CM has a non-integer real or
03110 // imaginary part.  Also extract the largest and smallest (real or
03111 // imaginary) values and return them in MAX_VAL and MIN_VAL.
03112 
03113 bool
03114 FloatComplexMatrix::all_integers (float& max_val, float& min_val) const
03115 {
03116   octave_idx_type nr = rows ();
03117   octave_idx_type nc = cols ();
03118 
03119   if (nr > 0 && nc > 0)
03120     {
03121       FloatComplex val = elem (0, 0);
03122 
03123       float r_val = std::real (val);
03124       float i_val = std::imag (val);
03125 
03126       max_val = r_val;
03127       min_val = r_val;
03128 
03129       if (i_val > max_val)
03130         max_val = i_val;
03131 
03132       if (i_val < max_val)
03133         min_val = i_val;
03134     }
03135   else
03136     return false;
03137 
03138   for (octave_idx_type j = 0; j < nc; j++)
03139     for (octave_idx_type i = 0; i < nr; i++)
03140       {
03141         FloatComplex val = elem (i, j);
03142 
03143         float r_val = std::real (val);
03144         float i_val = std::imag (val);
03145 
03146         if (r_val > max_val)
03147           max_val = r_val;
03148 
03149         if (i_val > max_val)
03150           max_val = i_val;
03151 
03152         if (r_val < min_val)
03153           min_val = r_val;
03154 
03155         if (i_val < min_val)
03156           min_val = i_val;
03157 
03158         if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
03159           return false;
03160       }
03161 
03162   return true;
03163 }
03164 
03165 bool
03166 FloatComplexMatrix::too_large_for_float (void) const
03167 {
03168   octave_idx_type nr = rows ();
03169   octave_idx_type nc = cols ();
03170 
03171   for (octave_idx_type j = 0; j < nc; j++)
03172     for (octave_idx_type i = 0; i < nr; i++)
03173       {
03174         FloatComplex val = elem (i, j);
03175 
03176         float r_val = std::real (val);
03177         float i_val = std::imag (val);
03178 
03179         if ((! (xisnan (r_val) || xisinf (r_val))
03180              && fabs (r_val) > FLT_MAX)
03181             || (! (xisnan (i_val) || xisinf (i_val))
03182                 && fabs (i_val) > FLT_MAX))
03183           return true;
03184       }
03185 
03186   return false;
03187 }
03188 
03189 // FIXME Do these really belong here?  Maybe they should be
03190 // in a base class?
03191 
03192 boolMatrix
03193 FloatComplexMatrix::all (int dim) const
03194 {
03195   return do_mx_red_op<bool, FloatComplex> (*this, dim, mx_inline_all);
03196 }
03197 
03198 boolMatrix
03199 FloatComplexMatrix::any (int dim) const
03200 {
03201   return do_mx_red_op<bool, FloatComplex> (*this, dim, mx_inline_any);
03202 }
03203 
03204 FloatComplexMatrix
03205 FloatComplexMatrix::cumprod (int dim) const
03206 {
03207   return do_mx_cum_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_cumprod);
03208 }
03209 
03210 FloatComplexMatrix
03211 FloatComplexMatrix::cumsum (int dim) const
03212 {
03213   return do_mx_cum_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_cumsum);
03214 }
03215 
03216 FloatComplexMatrix
03217 FloatComplexMatrix::prod (int dim) const
03218 {
03219   return do_mx_red_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_prod);
03220 }
03221 
03222 FloatComplexMatrix
03223 FloatComplexMatrix::sum (int dim) const
03224 {
03225   return do_mx_red_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_sum);
03226 }
03227 
03228 FloatComplexMatrix
03229 FloatComplexMatrix::sumsq (int dim) const
03230 {
03231   return do_mx_red_op<float, FloatComplex> (*this, dim, mx_inline_sumsq);
03232 }
03233 
03234 FloatMatrix FloatComplexMatrix::abs (void) const
03235 {
03236   return do_mx_unary_map<float, FloatComplex, std::abs> (*this);
03237 }
03238 
03239 FloatComplexMatrix
03240 FloatComplexMatrix::diag (octave_idx_type k) const
03241 {
03242   return MArray<FloatComplex>::diag (k);
03243 }
03244 
03245 bool
03246 FloatComplexMatrix::row_is_real_only (octave_idx_type i) const
03247 {
03248   bool retval = true;
03249 
03250   octave_idx_type nc = columns ();
03251 
03252   for (octave_idx_type j = 0; j < nc; j++)
03253     {
03254       if (std::imag (elem (i, j)) != 0.0)
03255         {
03256           retval = false;
03257           break;
03258         }
03259     }
03260 
03261   return retval;
03262 }
03263 
03264 bool
03265 FloatComplexMatrix::column_is_real_only (octave_idx_type j) const
03266 {
03267   bool retval = true;
03268 
03269   octave_idx_type nr = rows ();
03270 
03271   for (octave_idx_type i = 0; i < nr; i++)
03272     {
03273       if (std::imag (elem (i, j)) != 0.0)
03274         {
03275           retval = false;
03276           break;
03277         }
03278     }
03279 
03280   return retval;
03281 }
03282 
03283 FloatComplexColumnVector
03284 FloatComplexMatrix::row_min (void) const
03285 {
03286   Array<octave_idx_type> dummy_idx;
03287   return row_min (dummy_idx);
03288 }
03289 
03290 FloatComplexColumnVector
03291 FloatComplexMatrix::row_min (Array<octave_idx_type>& idx_arg) const
03292 {
03293   FloatComplexColumnVector result;
03294 
03295   octave_idx_type nr = rows ();
03296   octave_idx_type nc = cols ();
03297 
03298   if (nr > 0 && nc > 0)
03299     {
03300       result.resize (nr);
03301       idx_arg.resize (dim_vector (nr, 1));
03302 
03303       for (octave_idx_type i = 0; i < nr; i++)
03304         {
03305           bool real_only = row_is_real_only (i);
03306 
03307           octave_idx_type idx_j;
03308 
03309           FloatComplex tmp_min;
03310 
03311           float abs_min = octave_Float_NaN;
03312 
03313           for (idx_j = 0; idx_j < nc; idx_j++)
03314             {
03315               tmp_min = elem (i, idx_j);
03316 
03317               if (! xisnan (tmp_min))
03318                 {
03319                   abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min);
03320                   break;
03321                 }
03322             }
03323 
03324           for (octave_idx_type j = idx_j+1; j < nc; j++)
03325             {
03326               FloatComplex tmp = elem (i, j);
03327 
03328               if (xisnan (tmp))
03329                 continue;
03330 
03331               float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
03332 
03333               if (abs_tmp < abs_min)
03334                 {
03335                   idx_j = j;
03336                   tmp_min = tmp;
03337                   abs_min = abs_tmp;
03338                 }
03339             }
03340 
03341           if (xisnan (tmp_min))
03342             {
03343               result.elem (i) = FloatComplex_NaN_result;
03344               idx_arg.elem (i) = 0;
03345             }
03346           else
03347             {
03348               result.elem (i) = tmp_min;
03349               idx_arg.elem (i) = idx_j;
03350             }
03351         }
03352     }
03353 
03354   return result;
03355 }
03356 
03357 FloatComplexColumnVector
03358 FloatComplexMatrix::row_max (void) const
03359 {
03360   Array<octave_idx_type> dummy_idx;
03361   return row_max (dummy_idx);
03362 }
03363 
03364 FloatComplexColumnVector
03365 FloatComplexMatrix::row_max (Array<octave_idx_type>& idx_arg) const
03366 {
03367   FloatComplexColumnVector result;
03368 
03369   octave_idx_type nr = rows ();
03370   octave_idx_type nc = cols ();
03371 
03372   if (nr > 0 && nc > 0)
03373     {
03374       result.resize (nr);
03375       idx_arg.resize (dim_vector (nr, 1));
03376 
03377       for (octave_idx_type i = 0; i < nr; i++)
03378         {
03379           bool real_only = row_is_real_only (i);
03380 
03381           octave_idx_type idx_j;
03382 
03383           FloatComplex tmp_max;
03384 
03385           float abs_max = octave_Float_NaN;
03386 
03387           for (idx_j = 0; idx_j < nc; idx_j++)
03388             {
03389               tmp_max = elem (i, idx_j);
03390 
03391               if (! xisnan (tmp_max))
03392                 {
03393                   abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max);
03394                   break;
03395                 }
03396             }
03397 
03398           for (octave_idx_type j = idx_j+1; j < nc; j++)
03399             {
03400               FloatComplex tmp = elem (i, j);
03401 
03402               if (xisnan (tmp))
03403                 continue;
03404 
03405               float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
03406 
03407               if (abs_tmp > abs_max)
03408                 {
03409                   idx_j = j;
03410                   tmp_max = tmp;
03411                   abs_max = abs_tmp;
03412                 }
03413             }
03414 
03415           if (xisnan (tmp_max))
03416             {
03417               result.elem (i) = FloatComplex_NaN_result;
03418               idx_arg.elem (i) = 0;
03419             }
03420           else
03421             {
03422               result.elem (i) = tmp_max;
03423               idx_arg.elem (i) = idx_j;
03424             }
03425         }
03426     }
03427 
03428   return result;
03429 }
03430 
03431 FloatComplexRowVector
03432 FloatComplexMatrix::column_min (void) const
03433 {
03434   Array<octave_idx_type> dummy_idx;
03435   return column_min (dummy_idx);
03436 }
03437 
03438 FloatComplexRowVector
03439 FloatComplexMatrix::column_min (Array<octave_idx_type>& idx_arg) const
03440 {
03441   FloatComplexRowVector result;
03442 
03443   octave_idx_type nr = rows ();
03444   octave_idx_type nc = cols ();
03445 
03446   if (nr > 0 && nc > 0)
03447     {
03448       result.resize (nc);
03449       idx_arg.resize (dim_vector (1, nc));
03450 
03451       for (octave_idx_type j = 0; j < nc; j++)
03452         {
03453           bool real_only = column_is_real_only (j);
03454 
03455           octave_idx_type idx_i;
03456 
03457           FloatComplex tmp_min;
03458 
03459           float abs_min = octave_Float_NaN;
03460 
03461           for (idx_i = 0; idx_i < nr; idx_i++)
03462             {
03463               tmp_min = elem (idx_i, j);
03464 
03465               if (! xisnan (tmp_min))
03466                 {
03467                   abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min);
03468                   break;
03469                 }
03470             }
03471 
03472           for (octave_idx_type i = idx_i+1; i < nr; i++)
03473             {
03474               FloatComplex tmp = elem (i, j);
03475 
03476               if (xisnan (tmp))
03477                 continue;
03478 
03479               float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
03480 
03481               if (abs_tmp < abs_min)
03482                 {
03483                   idx_i = i;
03484                   tmp_min = tmp;
03485                   abs_min = abs_tmp;
03486                 }
03487             }
03488 
03489           if (xisnan (tmp_min))
03490             {
03491               result.elem (j) = FloatComplex_NaN_result;
03492               idx_arg.elem (j) = 0;
03493             }
03494           else
03495             {
03496               result.elem (j) = tmp_min;
03497               idx_arg.elem (j) = idx_i;
03498             }
03499         }
03500     }
03501 
03502   return result;
03503 }
03504 
03505 FloatComplexRowVector
03506 FloatComplexMatrix::column_max (void) const
03507 {
03508   Array<octave_idx_type> dummy_idx;
03509   return column_max (dummy_idx);
03510 }
03511 
03512 FloatComplexRowVector
03513 FloatComplexMatrix::column_max (Array<octave_idx_type>& idx_arg) const
03514 {
03515   FloatComplexRowVector result;
03516 
03517   octave_idx_type nr = rows ();
03518   octave_idx_type nc = cols ();
03519 
03520   if (nr > 0 && nc > 0)
03521     {
03522       result.resize (nc);
03523       idx_arg.resize (dim_vector (1, nc));
03524 
03525       for (octave_idx_type j = 0; j < nc; j++)
03526         {
03527           bool real_only = column_is_real_only (j);
03528 
03529           octave_idx_type idx_i;
03530 
03531           FloatComplex tmp_max;
03532 
03533           float abs_max = octave_Float_NaN;
03534 
03535           for (idx_i = 0; idx_i < nr; idx_i++)
03536             {
03537               tmp_max = elem (idx_i, j);
03538 
03539               if (! xisnan (tmp_max))
03540                 {
03541                   abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max);
03542                   break;
03543                 }
03544             }
03545 
03546           for (octave_idx_type i = idx_i+1; i < nr; i++)
03547             {
03548               FloatComplex tmp = elem (i, j);
03549 
03550               if (xisnan (tmp))
03551                 continue;
03552 
03553               float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
03554 
03555               if (abs_tmp > abs_max)
03556                 {
03557                   idx_i = i;
03558                   tmp_max = tmp;
03559                   abs_max = abs_tmp;
03560                 }
03561             }
03562 
03563           if (xisnan (tmp_max))
03564             {
03565               result.elem (j) = FloatComplex_NaN_result;
03566               idx_arg.elem (j) = 0;
03567             }
03568           else
03569             {
03570               result.elem (j) = tmp_max;
03571               idx_arg.elem (j) = idx_i;
03572             }
03573         }
03574     }
03575 
03576   return result;
03577 }
03578 
03579 // i/o
03580 
03581 std::ostream&
03582 operator << (std::ostream& os, const FloatComplexMatrix& a)
03583 {
03584   for (octave_idx_type i = 0; i < a.rows (); i++)
03585     {
03586       for (octave_idx_type j = 0; j < a.cols (); j++)
03587         {
03588           os << " ";
03589           octave_write_complex (os, a.elem (i, j));
03590         }
03591       os << "\n";
03592     }
03593   return os;
03594 }
03595 
03596 std::istream&
03597 operator >> (std::istream& is, FloatComplexMatrix& a)
03598 {
03599   octave_idx_type nr = a.rows ();
03600   octave_idx_type nc = a.cols ();
03601 
03602   if (nr > 0 && nc > 0)
03603     {
03604       FloatComplex tmp;
03605       for (octave_idx_type i = 0; i < nr; i++)
03606         for (octave_idx_type j = 0; j < nc; j++)
03607           {
03608             tmp = octave_read_value<FloatComplex> (is);
03609             if (is)
03610               a.elem (i, j) = tmp;
03611             else
03612               goto done;
03613           }
03614     }
03615 
03616 done:
03617 
03618   return is;
03619 }
03620 
03621 FloatComplexMatrix
03622 Givens (const FloatComplex& x, const FloatComplex& y)
03623 {
03624   float cc;
03625   FloatComplex cs, temp_r;
03626 
03627   F77_FUNC (clartg, CLARTG) (x, y, cc, cs, temp_r);
03628 
03629   FloatComplexMatrix g (2, 2);
03630 
03631   g.elem (0, 0) = cc;
03632   g.elem (1, 1) = cc;
03633   g.elem (0, 1) = cs;
03634   g.elem (1, 0) = -conj (cs);
03635 
03636   return g;
03637 }
03638 
03639 FloatComplexMatrix
03640 Sylvester (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
03641            const FloatComplexMatrix& c)
03642 {
03643   FloatComplexMatrix retval;
03644 
03645   // FIXME -- need to check that a, b, and c are all the same
03646   // size.
03647 
03648   // Compute Schur decompositions
03649 
03650   FloatComplexSCHUR as (a, "U");
03651   FloatComplexSCHUR bs (b, "U");
03652 
03653   // Transform c to new coordinates.
03654 
03655   FloatComplexMatrix ua = as.unitary_matrix ();
03656   FloatComplexMatrix sch_a = as.schur_matrix ();
03657 
03658   FloatComplexMatrix ub = bs.unitary_matrix ();
03659   FloatComplexMatrix sch_b = bs.schur_matrix ();
03660 
03661   FloatComplexMatrix cx = ua.hermitian () * c * ub;
03662 
03663   // Solve the sylvester equation, back-transform, and return the
03664   // solution.
03665 
03666   octave_idx_type a_nr = a.rows ();
03667   octave_idx_type b_nr = b.rows ();
03668 
03669   float scale;
03670   octave_idx_type info;
03671 
03672   FloatComplex *pa = sch_a.fortran_vec ();
03673   FloatComplex *pb = sch_b.fortran_vec ();
03674   FloatComplex *px = cx.fortran_vec ();
03675 
03676   F77_XFCN (ctrsyl, CTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
03677                              F77_CONST_CHAR_ARG2 ("N", 1),
03678                              1, a_nr, b_nr, pa, a_nr, pb,
03679                              b_nr, px, a_nr, scale, info
03680                              F77_CHAR_ARG_LEN (1)
03681                              F77_CHAR_ARG_LEN (1)));
03682 
03683   // FIXME -- check info?
03684 
03685   retval = -ua * cx * ub.hermitian ();
03686 
03687   return retval;
03688 }
03689 
03690 FloatComplexMatrix
03691 operator * (const FloatComplexMatrix& m, const FloatMatrix& a)
03692 {
03693   if (m.columns () > std::min (m.rows (), a.columns ()) / 10)
03694     return FloatComplexMatrix (real (m) * a, imag (m) * a);
03695   else
03696     return m * FloatComplexMatrix (a);
03697 }
03698 
03699 FloatComplexMatrix
03700 operator * (const FloatMatrix& m, const FloatComplexMatrix& a)
03701 {
03702   if (a.rows () > std::min (m.rows (), a.columns ()) / 10)
03703     return FloatComplexMatrix (m * real (a), m * imag (a));
03704   else
03705     return m * FloatComplexMatrix (a);
03706 }
03707 
03708 /* Simple Dot Product, Matrix-Vector and Matrix-Matrix Unit tests
03709 %!assert(single([1+i 2+i 3+i]) * single([ 4+i ; 5+i ; 6+i]), single(29+21i), 5e-7)
03710 %!assert(single([1+i 2+i ; 3+i 4+i ]) * single([5+i ; 6+i]), single([15 + 14i ; 37 + 18i]), 5e-7)
03711 %!assert(single([1+i 2+i ; 3+i 4+i ]) * single([5+i 6+i ; 7+i 8+i]), single([17 + 15i 20 + 17i; 41 + 19i 48 + 21i]), 5e-7)
03712 %!assert(single([1 i])*single([i 0])', single(-i));
03713 */
03714 
03715 /* Test some simple identities
03716 %!shared M, cv, rv
03717 %! M = single(randn(10,10))+i*single(rand(10,10));
03718 %! cv = single(randn(10,1))+i*single(rand(10,1));
03719 %! rv = single(randn(1,10))+i*single(rand(1,10));
03720 %!assert([M*cv,M*cv],M*[cv,cv],5e-6)
03721 %!assert([M.'*cv,M.'*cv],M.'*[cv,cv],5e-6)
03722 %!assert([M'*cv,M'*cv],M'*[cv,cv],5e-6)
03723 %!assert([rv*M;rv*M],[rv;rv]*M,5e-6)
03724 %!assert([rv*M.';rv*M.'],[rv;rv]*M.',5e-6)
03725 %!assert([rv*M';rv*M'],[rv;rv]*M',5e-6)
03726 %!assert(2*rv*cv,[rv,rv]*[cv;cv],5e-6)
03727 */
03728 
03729 static char
03730 get_blas_trans_arg (bool trans, bool conj)
03731 {
03732   return trans ? (conj ? 'C' : 'T') : 'N';
03733 }
03734 
03735 // the general GEMM operation
03736 
03737 FloatComplexMatrix
03738 xgemm (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
03739        blas_trans_type transa, blas_trans_type transb)
03740 {
03741   FloatComplexMatrix retval;
03742 
03743   bool tra = transa != blas_no_trans, trb = transb != blas_no_trans;
03744   bool cja = transa == blas_conj_trans, cjb = transb == blas_conj_trans;
03745 
03746   octave_idx_type a_nr = tra ? a.cols () : a.rows ();
03747   octave_idx_type a_nc = tra ? a.rows () : a.cols ();
03748 
03749   octave_idx_type b_nr = trb ? b.cols () : b.rows ();
03750   octave_idx_type b_nc = trb ? b.rows () : b.cols ();
03751 
03752   if (a_nc != b_nr)
03753     gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
03754   else
03755     {
03756       if (a_nr == 0 || a_nc == 0 || b_nc == 0)
03757         retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
03758       else if (a.data () == b.data () && a_nr == b_nc && tra != trb)
03759         {
03760           octave_idx_type lda = a.rows ();
03761 
03762           // FIXME -- looking at the reference BLAS, it appears that it
03763           // should not be necessary to initialize the output matrix if
03764           // BETA is 0 in the call to CHERK, but ATLAS appears to
03765           // use the result matrix before zeroing the elements.
03766 
03767           retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
03768           FloatComplex *c = retval.fortran_vec ();
03769 
03770           const char ctra = get_blas_trans_arg (tra, cja);
03771           if (cja || cjb)
03772             {
03773               F77_XFCN (cherk, CHERK, (F77_CONST_CHAR_ARG2 ("U", 1),
03774                                        F77_CONST_CHAR_ARG2 (&ctra, 1),
03775                                        a_nr, a_nc, 1.0,
03776                                        a.data (), lda, 0.0, c, a_nr
03777                                        F77_CHAR_ARG_LEN (1)
03778                                        F77_CHAR_ARG_LEN (1)));
03779               for (octave_idx_type j = 0; j < a_nr; j++)
03780                 for (octave_idx_type i = 0; i < j; i++)
03781                   retval.xelem (j,i) = std::conj (retval.xelem (i,j));
03782             }
03783           else
03784             {
03785               F77_XFCN (csyrk, CSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),
03786                                        F77_CONST_CHAR_ARG2 (&ctra, 1),
03787                                        a_nr, a_nc, 1.0,
03788                                        a.data (), lda, 0.0, c, a_nr
03789                                        F77_CHAR_ARG_LEN (1)
03790                                        F77_CHAR_ARG_LEN (1)));
03791               for (octave_idx_type j = 0; j < a_nr; j++)
03792                 for (octave_idx_type i = 0; i < j; i++)
03793                   retval.xelem (j,i) = retval.xelem (i,j);
03794 
03795             }
03796 
03797         }
03798       else
03799         {
03800           octave_idx_type lda = a.rows (), tda = a.cols ();
03801           octave_idx_type ldb = b.rows (), tdb = b.cols ();
03802 
03803           retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
03804           FloatComplex *c = retval.fortran_vec ();
03805 
03806           if (b_nc == 1 && a_nr == 1)
03807             {
03808               if (cja == cjb)
03809                 {
03810                   F77_FUNC (xcdotu, XCDOTU) (a_nc, a.data (), 1, b.data (), 1, *c);
03811                   if (cja) *c = std::conj (*c);
03812                 }
03813               else if (cja)
03814                   F77_FUNC (xcdotc, XCDOTC) (a_nc, a.data (), 1, b.data (), 1, *c);
03815               else
03816                   F77_FUNC (xcdotc, XCDOTC) (a_nc, b.data (), 1, a.data (), 1, *c);
03817             }
03818           else if (b_nc == 1 && ! cjb)
03819             {
03820               const char ctra = get_blas_trans_arg (tra, cja);
03821               F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),
03822                                        lda, tda, 1.0,  a.data (), lda,
03823                                        b.data (), 1, 0.0, c, 1
03824                                        F77_CHAR_ARG_LEN (1)));
03825             }
03826           else if (a_nr == 1 && ! cja && ! cjb)
03827             {
03828               const char crevtrb = get_blas_trans_arg (! trb, cjb);
03829               F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),
03830                                        ldb, tdb, 1.0,  b.data (), ldb,
03831                                        a.data (), 1, 0.0, c, 1
03832                                        F77_CHAR_ARG_LEN (1)));
03833             }
03834           else
03835             {
03836               const char ctra = get_blas_trans_arg (tra, cja);
03837               const char ctrb = get_blas_trans_arg (trb, cjb);
03838               F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),
03839                                        F77_CONST_CHAR_ARG2 (&ctrb, 1),
03840                                        a_nr, b_nc, a_nc, 1.0, a.data (),
03841                                        lda, b.data (), ldb, 0.0, c, a_nr
03842                                        F77_CHAR_ARG_LEN (1)
03843                                        F77_CHAR_ARG_LEN (1)));
03844             }
03845         }
03846     }
03847 
03848   return retval;
03849 }
03850 
03851 FloatComplexMatrix
03852 operator * (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
03853 {
03854   return xgemm (a, b);
03855 }
03856 
03857 // FIXME -- it would be nice to share code among the min/max
03858 // functions below.
03859 
03860 #define EMPTY_RETURN_CHECK(T) \
03861   if (nr == 0 || nc == 0) \
03862     return T (nr, nc);
03863 
03864 FloatComplexMatrix
03865 min (const FloatComplex& c, const FloatComplexMatrix& m)
03866 {
03867   octave_idx_type nr = m.rows ();
03868   octave_idx_type nc = m.columns ();
03869 
03870   EMPTY_RETURN_CHECK (FloatComplexMatrix);
03871 
03872   FloatComplexMatrix result (nr, nc);
03873 
03874   for (octave_idx_type j = 0; j < nc; j++)
03875     for (octave_idx_type i = 0; i < nr; i++)
03876       {
03877         octave_quit ();
03878         result (i, j) = xmin (c, m (i, j));
03879       }
03880 
03881   return result;
03882 }
03883 
03884 FloatComplexMatrix
03885 min (const FloatComplexMatrix& m, const FloatComplex& c)
03886 {
03887   octave_idx_type nr = m.rows ();
03888   octave_idx_type nc = m.columns ();
03889 
03890   EMPTY_RETURN_CHECK (FloatComplexMatrix);
03891 
03892   FloatComplexMatrix result (nr, nc);
03893 
03894   for (octave_idx_type j = 0; j < nc; j++)
03895     for (octave_idx_type i = 0; i < nr; i++)
03896       {
03897         octave_quit ();
03898         result (i, j) = xmin (m (i, j), c);
03899       }
03900 
03901   return result;
03902 }
03903 
03904 FloatComplexMatrix
03905 min (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
03906 {
03907   octave_idx_type nr = a.rows ();
03908   octave_idx_type nc = a.columns ();
03909 
03910   if (nr != b.rows () || nc != b.columns ())
03911     {
03912       (*current_liboctave_error_handler)
03913         ("two-arg min expecting args of same size");
03914       return FloatComplexMatrix ();
03915     }
03916 
03917   EMPTY_RETURN_CHECK (FloatComplexMatrix);
03918 
03919   FloatComplexMatrix result (nr, nc);
03920 
03921   for (octave_idx_type j = 0; j < nc; j++)
03922     {
03923       int columns_are_real_only = 1;
03924       for (octave_idx_type i = 0; i < nr; i++)
03925         {
03926           octave_quit ();
03927           if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0)
03928             {
03929               columns_are_real_only = 0;
03930               break;
03931             }
03932         }
03933 
03934       if (columns_are_real_only)
03935         {
03936           for (octave_idx_type i = 0; i < nr; i++)
03937             result (i, j) = xmin (std::real (a (i, j)), std::real (b (i, j)));
03938         }
03939       else
03940         {
03941           for (octave_idx_type i = 0; i < nr; i++)
03942             {
03943               octave_quit ();
03944               result (i, j) = xmin (a (i, j), b (i, j));
03945             }
03946         }
03947     }
03948 
03949   return result;
03950 }
03951 
03952 FloatComplexMatrix
03953 max (const FloatComplex& c, const FloatComplexMatrix& m)
03954 {
03955   octave_idx_type nr = m.rows ();
03956   octave_idx_type nc = m.columns ();
03957 
03958   EMPTY_RETURN_CHECK (FloatComplexMatrix);
03959 
03960   FloatComplexMatrix result (nr, nc);
03961 
03962   for (octave_idx_type j = 0; j < nc; j++)
03963     for (octave_idx_type i = 0; i < nr; i++)
03964       {
03965         octave_quit ();
03966         result (i, j) = xmax (c, m (i, j));
03967       }
03968 
03969   return result;
03970 }
03971 
03972 FloatComplexMatrix
03973 max (const FloatComplexMatrix& m, const FloatComplex& c)
03974 {
03975   octave_idx_type nr = m.rows ();
03976   octave_idx_type nc = m.columns ();
03977 
03978   EMPTY_RETURN_CHECK (FloatComplexMatrix);
03979 
03980   FloatComplexMatrix result (nr, nc);
03981 
03982   for (octave_idx_type j = 0; j < nc; j++)
03983     for (octave_idx_type i = 0; i < nr; i++)
03984       {
03985         octave_quit ();
03986         result (i, j) = xmax (m (i, j), c);
03987       }
03988 
03989   return result;
03990 }
03991 
03992 FloatComplexMatrix
03993 max (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
03994 {
03995   octave_idx_type nr = a.rows ();
03996   octave_idx_type nc = a.columns ();
03997 
03998   if (nr != b.rows () || nc != b.columns ())
03999     {
04000       (*current_liboctave_error_handler)
04001         ("two-arg max expecting args of same size");
04002       return FloatComplexMatrix ();
04003     }
04004 
04005   EMPTY_RETURN_CHECK (FloatComplexMatrix);
04006 
04007   FloatComplexMatrix result (nr, nc);
04008 
04009   for (octave_idx_type j = 0; j < nc; j++)
04010     {
04011       int columns_are_real_only = 1;
04012       for (octave_idx_type i = 0; i < nr; i++)
04013         {
04014           octave_quit ();
04015           if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0)
04016             {
04017               columns_are_real_only = 0;
04018               break;
04019             }
04020         }
04021 
04022       if (columns_are_real_only)
04023         {
04024           for (octave_idx_type i = 0; i < nr; i++)
04025             {
04026               octave_quit ();
04027               result (i, j) = xmax (std::real (a (i, j)), std::real (b (i, j)));
04028             }
04029         }
04030       else
04031         {
04032           for (octave_idx_type i = 0; i < nr; i++)
04033             {
04034               octave_quit ();
04035               result (i, j) = xmax (a (i, j), b (i, j));
04036             }
04037         }
04038     }
04039 
04040   return result;
04041 }
04042 
04043 FloatComplexMatrix linspace (const FloatComplexColumnVector& x1,
04044                              const FloatComplexColumnVector& x2,
04045                              octave_idx_type n)
04046 
04047 {
04048   if (n < 1) n = 1;
04049 
04050   octave_idx_type m = x1.length ();
04051 
04052   if (x2.length () != m)
04053     (*current_liboctave_error_handler) ("linspace: vectors must be of equal length");
04054 
04055   NoAlias<FloatComplexMatrix> retval;
04056 
04057   retval.clear (m, n);
04058   for (octave_idx_type i = 0; i < m; i++)
04059     retval(i, 0) = x1(i);
04060 
04061   // The last column is not needed while using delta.
04062   FloatComplex *delta = &retval(0, n-1);
04063   for (octave_idx_type i = 0; i < m; i++)
04064     delta[i] = (x2(i) - x1(i)) / (n - 1.0f);
04065 
04066   for (octave_idx_type j = 1; j < n-1; j++)
04067     for (octave_idx_type i = 0; i < m; i++)
04068       retval(i, j) = x1(i) + static_cast<float> (j)*delta[i];
04069 
04070   for (octave_idx_type i = 0; i < m; i++)
04071     retval(i, n-1) = x2(i);
04072 
04073   return retval;
04074 }
04075 
04076 MS_CMP_OPS (FloatComplexMatrix, FloatComplex)
04077 MS_BOOL_OPS (FloatComplexMatrix, FloatComplex)
04078 
04079 SM_CMP_OPS (FloatComplex, FloatComplexMatrix)
04080 SM_BOOL_OPS (FloatComplex, FloatComplexMatrix)
04081 
04082 MM_CMP_OPS (FloatComplexMatrix, FloatComplexMatrix)
04083 MM_BOOL_OPS (FloatComplexMatrix, FloatComplexMatrix)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines