メインページ   クラス階層   構成   ファイル一覧   構成メンバ   ファイルメンバ  

Array.cc

解説を見る。
00001 // Template array classes
00002 /*
00003 
00004 Copyright (C) 1996, 1997 John W. Eaton
00005 
00006 This file is part of Octave.
00007 
00008 Octave is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 2, or (at your option) any
00011 later version.
00012 
00013 Octave is distributed in the hope that it will be useful, but WITHOUT
00014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00016 for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with Octave; see the file COPYING.  If not, write to the Free
00020 Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00021 
00022 */
00023 
00024 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00025 #pragma implementation
00026 #endif
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <cassert>
00033 #include <climits>
00034 
00035 #include <iostream>
00036 
00037 #include "Array.h"
00038 #include "Array-flags.h"
00039 #include "Array-util.h"
00040 #include "Range.h"
00041 #include "idx-vector.h"
00042 #include "lo-error.h"
00043 #include "lo-sstream.h"
00044 
00045 // One dimensional array class.  Handles the reference counting for
00046 // all the derived classes.
00047 
00048 template <class T>
00049 Array<T>::Array (const Array<T>& a, const dim_vector& dv)
00050   : rep (a.rep), dimensions (dv), idx (0), idx_count (0)
00051 {
00052   rep->count++;
00053 
00054   if (a.numel () < dv.numel ())
00055     (*current_liboctave_error_handler)
00056       ("Array::Array (const Array&, const dim_vector&): dimension mismatch");
00057 }
00058 
00059 template <class T>
00060 Array<T>::~Array (void)
00061 {
00062   if (--rep->count <= 0)
00063     delete rep;
00064 
00065   delete [] idx;
00066 }
00067 
00068 template <class T>
00069 Array<T>
00070 Array<T>::squeeze (void) const
00071 {
00072   Array<T> retval = *this;
00073 
00074   if (ndims () > 2)
00075     {
00076       bool dims_changed = false;
00077 
00078       dim_vector new_dimensions = dimensions;
00079 
00080       int k = 0;
00081 
00082       for (int i = 0; i < ndims (); i++)
00083         {
00084           if (dimensions(i) == 1)
00085             dims_changed = true;
00086           else
00087             new_dimensions(k++) = dimensions(i);
00088         }
00089 
00090       if (dims_changed)
00091         {
00092           switch (k)
00093             {
00094             case 0:
00095               new_dimensions = dim_vector (1, 1);
00096               break;
00097 
00098             case 1:
00099               {
00100                 int tmp = new_dimensions(0);
00101 
00102                 new_dimensions.resize (2);
00103 
00104                 new_dimensions(0) = tmp;
00105                 new_dimensions(1) = 1;
00106               }
00107               break;
00108 
00109             default:
00110               new_dimensions.resize (k);
00111               break;
00112             }
00113         }
00114 
00115       // XXX FIXME XXX -- it would be better if we did not have to do
00116       // this, so we could share the data while still having different
00117       // dimension vectors.
00118 
00119       retval.make_unique ();
00120 
00121       retval.dimensions = new_dimensions;
00122     }
00123 
00124   return retval;
00125 }
00126 
00127 // A guess (should be quite conservative).
00128 #define MALLOC_OVERHEAD 1024
00129 
00130 template <class T>
00131 int
00132 Array<T>::get_size (int r, int c)
00133 {
00134   // XXX KLUGE XXX
00135 
00136   // If an allocation of an array with r * c elements of type T
00137   // would cause an overflow in the allocator when computing the
00138   // size of the allocation, then return a value which, although
00139   // not equivalent to the actual request, should be too large for
00140   // most current hardware, but not so large to cause the
00141   // allocator to barf on computing retval * sizeof (T).
00142 
00143   static int nl;
00144   static double dl
00145     = frexp (static_cast<double>
00146              (INT_MAX - MALLOC_OVERHEAD) / sizeof (T), &nl);
00147 
00148   // This value should be an integer.  If we return this value and
00149   // things work the way we expect, we should be paying a visit to
00150   // new_handler in no time flat.
00151   static int max_items = static_cast<int> (ldexp (dl, nl));
00152 
00153   int nr, nc;
00154   double dr = frexp (static_cast<double> (r), &nr);
00155   double dc = frexp (static_cast<double> (c), &nc);
00156 
00157   int nt = nr + nc;
00158   double dt = dr * dc;
00159 
00160   if (dt < 0.5)
00161     {
00162       nt--;
00163       dt *= 2;
00164     }
00165 
00166   return (nt < nl || (nt == nl && dt < dl)) ? r * c : max_items;
00167 }
00168 
00169 template <class T>
00170 int
00171 Array<T>::get_size (int r, int c, int p)
00172 {
00173   // XXX KLUGE XXX
00174 
00175   // If an allocation of an array with r * c * p elements of type T
00176   // would cause an overflow in the allocator when computing the
00177   // size of the allocation, then return a value which, although
00178   // not equivalent to the actual request, should be too large for
00179   // most current hardware, but not so large to cause the
00180   // allocator to barf on computing retval * sizeof (T).
00181 
00182   static int nl;
00183   static double dl
00184     = frexp (static_cast<double>
00185              (INT_MAX - MALLOC_OVERHEAD) / sizeof (T), &nl);
00186 
00187   // This value should be an integer.  If we return this value and
00188   // things work the way we expect, we should be paying a visit to
00189   // new_handler in no time flat.
00190   static int max_items = static_cast<int> (ldexp (dl, nl));
00191 
00192   int nr, nc, np;
00193   double dr = frexp (static_cast<double> (r), &nr);
00194   double dc = frexp (static_cast<double> (c), &nc);
00195   double dp = frexp (static_cast<double> (p), &np);
00196 
00197   int nt = nr + nc + np;
00198   double dt = dr * dc * dp;
00199 
00200   if (dt < 0.5)
00201     {
00202       nt--;
00203       dt *= 2;
00204 
00205       if (dt < 0.5)
00206         {
00207           nt--;
00208           dt *= 2;
00209         }
00210     }
00211 
00212   return (nt < nl || (nt == nl && dt < dl)) ? r * c * p : max_items;
00213 }
00214 
00215 template <class T>
00216 int
00217 Array<T>::get_size (const dim_vector& ra_idx)
00218 {
00219   // XXX KLUGE XXX
00220 
00221   // If an allocation of an array with r * c elements of type T
00222   // would cause an overflow in the allocator when computing the
00223   // size of the allocation, then return a value which, although
00224   // not equivalent to the actual request, should be too large for
00225   // most current hardware, but not so large to cause the
00226   // allocator to barf on computing retval * sizeof (T).
00227 
00228   static int nl;
00229   static double dl
00230     = frexp (static_cast<double>
00231              (INT_MAX - MALLOC_OVERHEAD) / sizeof (T), &nl);
00232 
00233   // This value should be an integer.  If we return this value and
00234   // things work the way we expect, we should be paying a visit to
00235   // new_handler in no time flat.
00236 
00237   static int max_items = static_cast<int> (ldexp (dl, nl));
00238 
00239   int retval = max_items;
00240 
00241   int n = ra_idx.length ();
00242 
00243   int nt = 0;
00244   double dt = 1;
00245 
00246   for (int i = 0; i < n; i++)
00247     {
00248       int nra_idx;
00249       double dra_idx = frexp (static_cast<double> (ra_idx(i)), &nra_idx);
00250 
00251       nt += nra_idx;
00252       dt *= dra_idx;
00253 
00254       if (dt < 0.5)
00255         {
00256           nt--;
00257           dt *= 2;
00258         }
00259     }
00260 
00261   if (nt < nl || (nt == nl && dt < dl))
00262     {
00263       retval = 1;
00264 
00265       for (int i = 0; i < n; i++)
00266         retval *= ra_idx(i);
00267     }
00268 
00269   return retval;
00270 }
00271 
00272 #undef MALLOC_OVERHEAD
00273 
00274 template <class T>
00275 int
00276 Array<T>::compute_index (const Array<int>& ra_idx) const
00277 {
00278   int retval = -1;
00279 
00280   int n = dimensions.length ();
00281 
00282   if (n > 0 && n == ra_idx.length ())
00283     {
00284       retval = ra_idx(--n);
00285 
00286       while (--n >= 0)
00287         {
00288           retval *= dimensions(n);
00289           retval += ra_idx(n);
00290         }
00291     }
00292   else
00293     (*current_liboctave_error_handler)
00294       ("Array<T>::compute_index: invalid ra_idxing operation");
00295 
00296   return retval;
00297 }
00298 
00299 template <class T>
00300 T
00301 Array<T>::range_error (const char *fcn, int n) const
00302 {
00303   (*current_liboctave_error_handler) ("%s (%d): range error", fcn, n);
00304   return T ();
00305 }
00306 
00307 template <class T>
00308 T&
00309 Array<T>::range_error (const char *fcn, int n)
00310 {
00311   (*current_liboctave_error_handler) ("%s (%d): range error", fcn, n);
00312   static T foo;
00313   return foo;
00314 }
00315 
00316 template <class T>
00317 T
00318 Array<T>::range_error (const char *fcn, int i, int j) const
00319 {
00320   (*current_liboctave_error_handler)
00321     ("%s (%d, %d): range error", fcn, i, j);
00322   return T ();
00323 }
00324 
00325 template <class T>
00326 T&
00327 Array<T>::range_error (const char *fcn, int i, int j)
00328 {
00329   (*current_liboctave_error_handler)
00330     ("%s (%d, %d): range error", fcn, i, j);
00331   static T foo;
00332   return foo;
00333 }
00334 
00335 template <class T>
00336 T
00337 Array<T>::range_error (const char *fcn, int i, int j, int k) const
00338 {
00339   (*current_liboctave_error_handler)
00340     ("%s (%d, %d, %d): range error", fcn, i, j, k);
00341   return T ();
00342 }
00343 
00344 template <class T>
00345 T&
00346 Array<T>::range_error (const char *fcn, int i, int j, int k)
00347 {
00348   (*current_liboctave_error_handler)
00349     ("%s (%d, %d, %d): range error", fcn, i, j, k);
00350   static T foo;
00351   return foo;
00352 }
00353 
00354 template <class T>
00355 T
00356 Array<T>::range_error (const char *fcn, const Array<int>& ra_idx) const
00357 {
00358   OSSTREAM buf;
00359 
00360   buf << fcn << " (";
00361 
00362   int n = ra_idx.length ();
00363 
00364   if (n > 0)
00365     buf << ra_idx(0);
00366 
00367   for (int i = 1; i < n; i++)
00368     buf << ", " << ra_idx(i);
00369 
00370   buf << "): range error";
00371 
00372   buf << OSSTREAM_ENDS;
00373 
00374   (*current_liboctave_error_handler) (OSSTREAM_C_STR (buf));
00375 
00376   OSSTREAM_FREEZE (buf);
00377 
00378   return T ();
00379 }
00380 
00381 template <class T>
00382 T&
00383 Array<T>::range_error (const char *fcn, const Array<int>& ra_idx)
00384 {
00385   OSSTREAM buf;
00386 
00387   buf << fcn << " (";
00388 
00389   int n = ra_idx.length ();
00390 
00391   if (n > 0)
00392     buf << ra_idx(0);
00393 
00394   for (int i = 1; i < n; i++)
00395     buf << ", " << ra_idx(i);
00396 
00397   buf << "): range error";
00398 
00399   buf << OSSTREAM_ENDS;
00400 
00401   (*current_liboctave_error_handler) (OSSTREAM_C_STR (buf));
00402 
00403   OSSTREAM_FREEZE (buf);
00404 
00405   static T foo;
00406   return foo;
00407 }
00408 
00409 template <class T>
00410 Array<T>
00411 Array<T>::reshape (const dim_vector& new_dims) const
00412 {
00413   Array<T> retval;
00414 
00415   if (dimensions != new_dims)
00416     {
00417       if (dimensions.numel () == new_dims.numel ())
00418         retval = Array<T> (*this, new_dims);
00419       else
00420         (*current_liboctave_error_handler) ("reshape: size mismatch");
00421     }
00422   else
00423     retval = *this;
00424 
00425   return retval;
00426 }
00427 
00428 template <class T>
00429 Array<T>
00430 Array<T>::permute (const Array<int>& perm_vec, bool inv) const
00431 {
00432   Array<T> retval;
00433 
00434   dim_vector dv = dims ();
00435   dim_vector dv_new;
00436 
00437   int nd = dv.length ();
00438 
00439   dv_new.resize (nd);
00440 
00441   // Need this array to check for identical elements in permutation array.
00442   Array<bool> checked (nd, false);
00443 
00444   // Find dimension vector of permuted array.
00445   for (int i = 0; i < nd; i++)
00446     {
00447       int perm_el = perm_vec.elem (i);
00448 
00449       if (perm_el > dv.length () || perm_el < 1)
00450         {
00451           (*current_liboctave_error_handler)
00452             ("permutation vector contains an invalid element");
00453 
00454           return retval;
00455         }
00456 
00457       if (checked.elem(perm_el - 1))
00458         {
00459           (*current_liboctave_error_handler)
00460             ("PERM cannot contain identical elements");
00461 
00462           return retval;
00463         }
00464       else
00465         checked.elem(perm_el - 1) = true;
00466 
00467       dv_new (i) = dv (perm_el - 1);
00468     }
00469 
00470   retval.resize (dv_new);
00471 
00472   // Index array to the original array.
00473   Array<int> old_idx (nd, 0);
00474 
00475   // Number of elements in Array (should be the same for
00476   // both the permuted array and original array).
00477   int n = retval.length ();
00478 
00479   // Permute array.
00480   for (int i = 0; i < n; i++)
00481     {
00482       // Get the idx of permuted array.
00483       Array<int> new_idx = calc_permutated_idx (old_idx, perm_vec, inv);
00484 
00485       retval.elem (new_idx) = elem (old_idx);
00486 
00487       increment_index (old_idx, dv);
00488     }
00489 
00490   return retval;
00491 }
00492 
00493 template <class T>
00494 void
00495 Array<T>::resize_no_fill (int n)
00496 {
00497   if (n < 0)
00498     {
00499       (*current_liboctave_error_handler)
00500         ("can't resize to negative dimension");
00501       return;
00502     }
00503 
00504   if (n == length ())
00505     return;
00506 
00507   typename Array<T>::ArrayRep *old_rep = rep;
00508   const T *old_data = data ();
00509   int old_len = length ();
00510 
00511   rep = new typename Array<T>::ArrayRep (n);
00512 
00513   dimensions = dim_vector (n);
00514 
00515   if (n > 0 && old_data && old_len > 0)
00516     {
00517       int min_len = old_len < n ? old_len : n;
00518 
00519       for (int i = 0; i < min_len; i++)
00520         xelem (i) = old_data[i];
00521     }
00522 
00523   if (--old_rep->count <= 0)
00524     delete old_rep;
00525 }
00526 
00527 template <class T>
00528 void
00529 Array<T>::resize_no_fill (const dim_vector& dv)
00530 {
00531   int n = dv.length ();
00532 
00533   for (int i = 0; i < n; i++)
00534     {
00535       if (dv(i) < 0)
00536         {
00537           (*current_liboctave_error_handler)
00538             ("can't resize to negative dimension");
00539           return;
00540         }
00541     }
00542 
00543   bool same_size = true;
00544 
00545   if (dimensions.length () != n)
00546     {
00547       same_size = false;
00548     }
00549   else
00550     {
00551       for (int i = 0; i < n; i++)
00552         {
00553           if (dv(i) != dimensions(i))
00554             {
00555               same_size = false;
00556               break;
00557             }
00558         }
00559     }
00560 
00561   if (same_size)
00562     return;
00563 
00564   typename Array<T>::ArrayRep *old_rep = rep;
00565   const T *old_data = data ();
00566 
00567   int ts = get_size (dv);
00568 
00569   rep = new typename Array<T>::ArrayRep (ts);
00570 
00571   dim_vector dv_old = dimensions;
00572   int dv_old_orig_len = dv_old.length ();
00573   dimensions = dv;
00574   int ts_old = get_size (dv_old);
00575 
00576   if (ts > 0 && ts_old > 0 && dv_old_orig_len > 0)
00577     {
00578       Array<int> ra_idx (dimensions.length (), 0);
00579 
00580       if (n > dv_old_orig_len)
00581         {
00582           dv_old.resize (n);
00583 
00584           for (int i = dv_old_orig_len; i < n; i++)
00585             dv_old.elem (i) = 1;
00586         }
00587 
00588       for (int i = 0; i < ts; i++)
00589         {
00590           if (index_in_bounds (ra_idx, dv_old))
00591             rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)];
00592 
00593           increment_index (ra_idx, dimensions);
00594         }
00595     }
00596 
00597   if (--old_rep->count <= 0)
00598     delete old_rep;
00599 }
00600 
00601 template <class T>
00602 void
00603 Array<T>::resize_no_fill (int r, int c)
00604 {
00605   if (r < 0 || c < 0)
00606     {
00607       (*current_liboctave_error_handler)
00608         ("can't resize to negative dimension");
00609       return;
00610     }
00611 
00612   int n = ndims ();
00613 
00614   if (n == 0)
00615     dimensions = dim_vector (0, 0);
00616 
00617   assert (ndims () == 2);
00618 
00619   if (r == dim1 () && c == dim2 ())
00620     return;
00621 
00622   typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
00623   const T *old_data = data ();
00624 
00625   int old_d1 = dim1 ();
00626   int old_d2 = dim2 ();
00627   int old_len = length ();
00628 
00629   int ts = get_size (r, c);
00630 
00631   rep = new typename Array<T>::ArrayRep (ts);
00632 
00633   dimensions = dim_vector (r, c);
00634 
00635   if (ts > 0 && old_data && old_len > 0)
00636     {
00637       int min_r = old_d1 < r ? old_d1 : r;
00638       int min_c = old_d2 < c ? old_d2 : c;
00639 
00640       for (int j = 0; j < min_c; j++)
00641         for (int i = 0; i < min_r; i++)
00642           xelem (i, j) = old_data[old_d1*j+i];
00643     }
00644 
00645   if (--old_rep->count <= 0)
00646     delete old_rep;
00647 }
00648 
00649 template <class T>
00650 void
00651 Array<T>::resize_no_fill (int r, int c, int p)
00652 {
00653   if (r < 0 || c < 0 || p < 0)
00654     {
00655       (*current_liboctave_error_handler)
00656         ("can't resize to negative dimension");
00657       return;
00658     }
00659 
00660   int n = ndims ();
00661 
00662   if (n == 0)
00663     dimensions = dim_vector (0, 0, 0);
00664 
00665   assert (ndims () == 3);
00666 
00667   if (r == dim1 () && c == dim2 () && p == dim3 ())
00668     return;
00669 
00670   typename Array<T>::ArrayRep *old_rep = rep;
00671   const T *old_data = data ();
00672 
00673   int old_d1 = dim1 ();
00674   int old_d2 = dim2 ();
00675   int old_d3 = dim3 ();
00676   int old_len = length ();
00677 
00678   int ts = get_size (get_size (r, c), p);
00679 
00680   rep = new typename Array<T>::ArrayRep (ts);
00681 
00682   dimensions = dim_vector (r, c, p);
00683 
00684   if (ts > 0 && old_data && old_len > 0)
00685     {
00686       int min_r = old_d1 < r ? old_d1 : r;
00687       int min_c = old_d2 < c ? old_d2 : c;
00688       int min_p = old_d3 < p ? old_d3 : p;
00689 
00690       for (int k = 0; k < min_p; k++)
00691         for (int j = 0; j < min_c; j++)
00692           for (int i = 0; i < min_r; i++)
00693             xelem (i, j, k) = old_data[old_d1*(old_d2*k+j)+i];
00694     }
00695 
00696   if (--old_rep->count <= 0)
00697     delete old_rep;
00698 }
00699 
00700 template <class T>
00701 void
00702 Array<T>::resize_and_fill (int n, const T& val)
00703 {
00704   if (n < 0)
00705     {
00706       (*current_liboctave_error_handler)
00707         ("can't resize to negative dimension");
00708       return;
00709     }
00710 
00711   if (n == length ())
00712     return;
00713 
00714   typename Array<T>::ArrayRep *old_rep = rep;
00715   const T *old_data = data ();
00716   int old_len = length ();
00717 
00718   rep = new typename Array<T>::ArrayRep (n);
00719 
00720   dimensions = dim_vector (n);
00721 
00722   if (n > 0)
00723     {
00724       int min_len = old_len < n ? old_len : n;
00725 
00726       if (old_data && old_len > 0)
00727         {
00728           for (int i = 0; i < min_len; i++)
00729             xelem (i) = old_data[i];
00730         }
00731 
00732       for (int i = old_len; i < n; i++)
00733         xelem (i) = val;
00734     }
00735 
00736   if (--old_rep->count <= 0)
00737     delete old_rep;
00738 }
00739 
00740 template <class T>
00741 void
00742 Array<T>::resize_and_fill (int r, int c, const T& val)
00743 {
00744   if (r < 0 || c < 0)
00745     {
00746       (*current_liboctave_error_handler)
00747         ("can't resize to negative dimension");
00748       return;
00749     }
00750 
00751   if (ndims () == 0)
00752     dimensions = dim_vector (0, 0);
00753 
00754   assert (ndims () == 2);
00755 
00756   if (r == dim1 () && c == dim2 ())
00757     return;
00758 
00759   typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
00760   const T *old_data = data ();
00761 
00762   int old_d1 = dim1 ();
00763   int old_d2 = dim2 ();
00764   int old_len = length ();
00765 
00766   int ts = get_size (r, c);
00767 
00768   rep = new typename Array<T>::ArrayRep (ts);
00769 
00770   dimensions = dim_vector (r, c);
00771 
00772   if (ts > 0)
00773     {
00774       int min_r = old_d1 < r ? old_d1 : r;
00775       int min_c = old_d2 < c ? old_d2 : c;
00776 
00777       if (old_data && old_len > 0)
00778         {
00779           for (int j = 0; j < min_c; j++)
00780             for (int i = 0; i < min_r; i++)
00781               xelem (i, j) = old_data[old_d1*j+i];
00782         }
00783 
00784       for (int j = 0; j < min_c; j++)
00785         for (int i = min_r; i < r; i++)
00786           xelem (i, j) = val;
00787 
00788       for (int j = min_c; j < c; j++)
00789         for (int i = 0; i < r; i++)
00790           xelem (i, j) = val;
00791     }
00792 
00793   if (--old_rep->count <= 0)
00794     delete old_rep;
00795 }
00796 
00797 template <class T>
00798 void
00799 Array<T>::resize_and_fill (int r, int c, int p, const T& val)
00800 {
00801   if (r < 0 || c < 0 || p < 0)
00802     {
00803       (*current_liboctave_error_handler)
00804         ("can't resize to negative dimension");
00805       return;
00806     }
00807 
00808   if (ndims () == 0)
00809     dimensions = dim_vector (0, 0, 0);
00810 
00811   assert (ndims () == 3);
00812 
00813   if (r == dim1 () && c == dim2 () && p == dim3 ())
00814     return;
00815 
00816   typename Array<T>::ArrayRep *old_rep = rep;
00817   const T *old_data = data ();
00818 
00819   int old_d1 = dim1 ();
00820   int old_d2 = dim2 ();
00821   int old_d3 = dim3 ();
00822 
00823   int old_len = length ();
00824 
00825   int ts = get_size (get_size (r, c), p);
00826 
00827   rep = new typename Array<T>::ArrayRep (ts);
00828 
00829   dimensions = dim_vector (r, c, p);
00830 
00831   if (ts > 0)
00832     {
00833       int min_r = old_d1 < r ? old_d1 : r;
00834       int min_c = old_d2 < c ? old_d2 : c;
00835       int min_p = old_d3 < p ? old_d3 : p;
00836 
00837       if (old_data && old_len > 0)
00838         for (int k = 0; k < min_p; k++)
00839           for (int j = 0; j < min_c; j++)
00840             for (int i = 0; i < min_r; i++)
00841               xelem (i, j, k) = old_data[old_d1*(old_d2*k+j)+i];
00842 
00843       // XXX FIXME XXX -- if the copy constructor is expensive, this
00844       // may win.  Otherwise, it may make more sense to just copy the
00845       // value everywhere when making the new ArrayRep.
00846 
00847       for (int k = 0; k < min_p; k++)
00848         for (int j = min_c; j < c; j++)
00849           for (int i = 0; i < min_r; i++)
00850             xelem (i, j, k) = val;
00851 
00852       for (int k = 0; k < min_p; k++)
00853         for (int j = 0; j < c; j++)
00854           for (int i = min_r; i < r; i++)
00855             xelem (i, j, k) = val;
00856 
00857       for (int k = min_p; k < p; k++)
00858         for (int j = 0; j < c; j++)
00859           for (int i = 0; i < r; i++)
00860             xelem (i, j, k) = val;
00861     }
00862 
00863   if (--old_rep->count <= 0)
00864     delete old_rep;
00865 }
00866 
00867 template <class T>
00868 void
00869 Array<T>::resize_and_fill (const dim_vector& dv, const T& val)
00870 {
00871   int n = dv.length ();
00872 
00873   for (int i = 0; i < n; i++)
00874     {
00875       if (dv(i) < 0)
00876         {
00877           (*current_liboctave_error_handler)
00878             ("can't resize to negative dimension");
00879           return;
00880         }
00881     }
00882 
00883   bool same_size = true;
00884 
00885   if (dimensions.length () != n)
00886     {
00887       same_size = false;
00888     }
00889   else
00890     {
00891       for (int i = 0; i < n; i++)
00892         {
00893           if (dv(i) != dimensions(i))
00894             {
00895               same_size = false;
00896               break;
00897             }
00898         }
00899     }
00900 
00901   if (same_size)
00902     return;
00903 
00904   typename Array<T>::ArrayRep *old_rep = rep;
00905   const T *old_data = data ();
00906 
00907   int len = get_size (dv);
00908 
00909   rep = new typename Array<T>::ArrayRep (len);
00910 
00911   dim_vector dv_old = dimensions;
00912   int dv_old_orig_len = dv_old.length ();
00913   dimensions = dv;
00914 
00915   if (len > 0 && dv_old_orig_len > 0)
00916     {
00917       Array<int> ra_idx (dimensions.length (), 0);
00918       
00919       if (n > dv_old_orig_len)
00920         {
00921           dv_old.resize (n);
00922 
00923           for (int i = dv_old_orig_len; i < n; i++)
00924             dv_old.elem (i) = 1;
00925         }
00926 
00927       for (int i = 0; i < len; i++)
00928         {
00929           if (index_in_bounds (ra_idx, dv_old))
00930             rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)];
00931           else
00932             rep->elem (i) = val;
00933           
00934           increment_index (ra_idx, dimensions);
00935         }
00936     }
00937   else
00938     for (int i = 0; i < len; i++)
00939       rep->elem (i) = val;
00940 
00941   if (--old_rep->count <= 0)
00942     delete old_rep;
00943 }
00944 
00945 template <class T>
00946 Array<T>&
00947 Array<T>::insert (const Array<T>& a, int r, int c)
00948 {
00949   if (ndims () == 2 && a.ndims () == 2)
00950     insert2 (a, r, c);
00951   else
00952     insertN (a, r, c);
00953 
00954   return *this;
00955 }
00956 
00957 
00958 template <class T>
00959 Array<T>&
00960 Array<T>::insert2 (const Array<T>& a, int r, int c)
00961 {
00962   int a_rows = a.rows ();
00963   int a_cols = a.cols ();
00964 
00965   if (r < 0 || r + a_rows > rows () || c < 0 || c + a_cols > cols ())
00966     {
00967       (*current_liboctave_error_handler) ("range error for insert");
00968       return *this;
00969     }
00970 
00971   for (int j = 0; j < a_cols; j++)
00972     for (int i = 0; i < a_rows; i++)
00973       elem (r+i, c+j) = a.elem (i, j);
00974 
00975   return *this;
00976 }
00977 
00978 template <class T>
00979 Array<T>&
00980 Array<T>::insertN (const Array<T>& a, int r, int c)
00981 {
00982   dim_vector dv = dims ();
00983 
00984   dim_vector a_dv = a.dims ();
00985 
00986   int n = a_dv.length ();
00987 
00988   if (n == dimensions.length ())
00989     {
00990       Array<int> a_ra_idx (a_dv.length (), 0);
00991 
00992       a_ra_idx.elem (0) = r;
00993       a_ra_idx.elem (1) = c;
00994 
00995       for (int i = 0; i < n; i++)
00996         {
00997           if (a_ra_idx(i) < 0 || (a_ra_idx(i) + a_dv(i)) > dv(i))
00998             {
00999               (*current_liboctave_error_handler)
01000                 ("Array<T>::insert: range error for insert");
01001               return *this;
01002             }
01003         }
01004 
01005       int n_elt = a.numel ();
01006       
01007       const T *a_data = a.data ();   
01008    
01009       int iidx = 0;
01010           
01011       int a_rows = a_dv(0);
01012 
01013       int this_rows = dv(0);
01014           
01015       int numel_page = a_dv(0) * a_dv(1);         
01016 
01017       int count_pages = 0;
01018           
01019       for (int i = 0; i < n_elt; i++)
01020         {
01021           if (i != 0 && i % a_rows == 0)
01022             iidx += (this_rows - a_rows);             
01023           
01024           if (i % numel_page == 0)
01025             iidx = c * dv(0) + r + dv(0) * dv(1) * count_pages++;
01026 
01027           elem (iidx++) = a_data[i];
01028         }
01029     }
01030   else
01031     (*current_liboctave_error_handler)
01032       ("Array<T>::insert: invalid indexing operation");
01033 
01034   return *this;
01035 }
01036 
01037 template <class T>
01038 Array<T>&
01039 Array<T>::insert (const Array<T>& a, const Array<int>& ra_idx)
01040 {
01041   int n = ra_idx.length ();
01042 
01043   if (n == dimensions.length ())
01044     {
01045       dim_vector dva = a.dims ();
01046       dim_vector dv = dims ();
01047       int len_a = dva.length ();
01048 
01049       for (int i = 0; i < n; i++)
01050         {
01051           if (ra_idx(i) < 0 || (ra_idx(i) + 
01052                                 (i < len_a ? dva(i) : 1)) > dimensions(i))
01053             {
01054               (*current_liboctave_error_handler)
01055                 ("Array<T>::insert: range error for insert");
01056               return *this;
01057             }
01058         }
01059 
01060       if (dva.numel ())
01061         {
01062           const T *a_data = a.data ();   
01063           int numel_to_move = dva (0);
01064           int skip = dv (0);
01065           for (int i = 0; i < len_a - 1; i++)
01066             if (ra_idx(i) == 0 && dva(i) == dv(i))
01067               {
01068                 numel_to_move *= dva(i+1);
01069                 skip *= dv(i+1);
01070               }
01071             else
01072               {
01073                 skip -= dva(i);
01074                 break;
01075               }
01076 
01077           int jidx = ra_idx (n - 1);
01078           for (int i = n-2; i >= 0; i--)
01079             {
01080               jidx *= dv (i);
01081               jidx += ra_idx (i);
01082             }
01083 
01084           int iidx = 0;
01085           int moves = dva.numel () / numel_to_move;
01086           for (int i = 0; i < moves; i++)
01087             {
01088               for (int j = 0; j < numel_to_move; j++)
01089                 elem (jidx++) = a_data[iidx++];
01090               jidx += skip;
01091             }
01092         }
01093     }
01094   else
01095     (*current_liboctave_error_handler)
01096       ("Array<T>::insert: invalid indexing operation");
01097 
01098   return *this;
01099 }
01100 
01101 template <class T>
01102 Array<T>
01103 Array<T>::transpose (void) const
01104 {
01105   assert (ndims () == 2);
01106 
01107   int nr = dim1 ();
01108   int nc = dim2 ();
01109 
01110   if (nr > 1 && nc > 1)
01111     {
01112       Array<T> result (dim_vector (nc, nr));
01113 
01114       for (int j = 0; j < nc; j++)
01115         for (int i = 0; i < nr; i++)
01116           result.xelem (j, i) = xelem (i, j);
01117 
01118       return result;
01119     }
01120   else
01121     {
01122       // Fast transpose for vectors and empty matrices
01123       return Array<T> (*this, dim_vector (nc, nr));
01124     }
01125 }
01126 
01127 template <class T>
01128 T *
01129 Array<T>::fortran_vec (void)
01130 {
01131   if (rep->count > 1)
01132     {
01133       --rep->count;
01134       rep = new typename Array<T>::ArrayRep (*rep);
01135     }
01136   return rep->data;
01137 }
01138 
01139 template <class T>
01140 void
01141 Array<T>::maybe_delete_dims (void)
01142 {
01143   int nd = dimensions.length ();
01144 
01145   dim_vector new_dims (1, 1);
01146 
01147   bool delete_dims = true;
01148 
01149   for (int i = nd - 1; i >= 0; i--)
01150     {
01151       if (delete_dims)
01152         {
01153           if (dimensions(i) != 1)
01154             {
01155               delete_dims = false;
01156 
01157               new_dims = dim_vector (i + 1, dimensions(i));
01158             }
01159         }
01160       else
01161         new_dims(i) = dimensions(i);
01162     }
01163 
01164   if (nd != new_dims.length ())
01165     dimensions = new_dims;
01166 }
01167 
01168 template <class T>
01169 void
01170 Array<T>::clear_index (void)
01171 {
01172   delete [] idx;
01173   idx = 0;
01174   idx_count = 0;
01175 }
01176 
01177 template <class T>
01178 void
01179 Array<T>::set_index (const idx_vector& idx_arg)
01180 {
01181   int nd = ndims ();
01182 
01183   if (! idx && nd > 0)
01184     idx = new idx_vector [nd];
01185 
01186   if (idx_count < nd)
01187     {
01188       idx[idx_count++] = idx_arg;
01189     }
01190   else
01191     {
01192       idx_vector *new_idx = new idx_vector [idx_count+1];
01193 
01194       for (int i = 0; i < idx_count; i++)
01195         new_idx[i] = idx[i];
01196 
01197       new_idx[idx_count++] = idx_arg;
01198 
01199       delete [] idx;
01200 
01201       idx = new_idx;
01202     }
01203 }
01204 
01205 template <class T>
01206 void
01207 Array<T>::maybe_delete_elements (idx_vector& idx_arg)
01208 {
01209   switch (ndims ())
01210     {
01211     case 1:
01212       maybe_delete_elements_1 (idx_arg);
01213       break;
01214 
01215     case 2:
01216       maybe_delete_elements_2 (idx_arg);
01217       break;
01218 
01219     default:
01220       (*current_liboctave_error_handler)
01221         ("Array<T>::maybe_delete_elements: invalid operation");
01222       break;
01223     }
01224 }
01225 
01226 template <class T>
01227 void
01228 Array<T>::maybe_delete_elements_1 (idx_vector& idx_arg)
01229 {
01230   int len = length ();
01231 
01232   if (len == 0)
01233     return;
01234 
01235   if (idx_arg.is_colon_equiv (len, 1))
01236     resize_no_fill (0);
01237   else
01238     {
01239       int num_to_delete = idx_arg.length (len);
01240 
01241       if (num_to_delete != 0)
01242         {
01243           int new_len = len;
01244 
01245           int iidx = 0;
01246 
01247           for (int i = 0; i < len; i++)
01248             if (i == idx_arg.elem (iidx))
01249               {
01250                 iidx++;
01251                 new_len--;
01252 
01253                 if (iidx == num_to_delete)
01254                   break;
01255               }
01256 
01257           if (new_len > 0)
01258             {
01259               T *new_data = new T [new_len];
01260 
01261               int ii = 0;
01262               iidx = 0;
01263               for (int i = 0; i < len; i++)
01264                 {
01265                   if (iidx < num_to_delete && i == idx_arg.elem (iidx))
01266                     iidx++;
01267                   else
01268                     {
01269                       new_data[ii] = elem (i);
01270                       ii++;
01271                     }
01272                 }
01273 
01274               if (--rep->count <= 0)
01275                 delete rep;
01276 
01277               rep = new typename Array<T>::ArrayRep (new_data, new_len);
01278 
01279               dimensions.resize (1);
01280               dimensions(0) = new_len;
01281             }
01282           else
01283             (*current_liboctave_error_handler)
01284               ("A(idx) = []: index out of range");
01285         }
01286     }
01287 }
01288 
01289 template <class T>
01290 void
01291 Array<T>::maybe_delete_elements_2 (idx_vector& idx_arg)
01292 {
01293   assert (ndims () == 2);
01294 
01295   int nr = dim1 ();
01296   int nc = dim2 ();
01297 
01298   if (nr == 0 && nc == 0)
01299     return;
01300 
01301   int n;
01302   if (nr == 1)
01303     n = nc;
01304   else if (nc == 1)
01305     n = nr;
01306   else
01307     {
01308       // Reshape to row vector for Matlab compatibility.
01309 
01310       n = nr * nc;
01311       nr = 1;
01312       nc = n;
01313     }
01314 
01315   if (idx_arg.is_colon_equiv (n, 1))
01316     {
01317       // Either A(:) = [] or A(idx) = [] with idx enumerating all
01318       // elements, so we delete all elements and return [](0x0).  To
01319       // preserve the orientation of the vector, you have to use
01320       // A(idx,:) = [] (delete rows) or A(:,idx) (delete columns).
01321 
01322       resize_no_fill (0, 0);
01323       return;
01324     }
01325 
01326   idx_arg.sort (true);
01327 
01328   int num_to_delete = idx_arg.length (n);
01329 
01330   if (num_to_delete != 0)
01331     {
01332       int new_n = n;
01333 
01334       int iidx = 0;
01335 
01336       for (int i = 0; i < n; i++)
01337         if (i == idx_arg.elem (iidx))
01338           {
01339             iidx++;
01340             new_n--;
01341 
01342             if (iidx == num_to_delete)
01343               break;
01344           }
01345 
01346       if (new_n > 0)
01347         {
01348           T *new_data = new T [new_n];
01349 
01350           int ii = 0;
01351           iidx = 0;
01352           for (int i = 0; i < n; i++)
01353             {
01354               if (iidx < num_to_delete && i == idx_arg.elem (iidx))
01355                 iidx++;
01356               else
01357                 {
01358                   new_data[ii] = elem (i);
01359 
01360                   ii++;
01361                 }
01362             }
01363 
01364           if (--(Array<T>::rep)->count <= 0)
01365             delete Array<T>::rep;
01366 
01367           Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_n);
01368 
01369           dimensions.resize (2);
01370 
01371           if (nr == 1)
01372             {
01373               dimensions(0) = 1;
01374               dimensions(1) = new_n;
01375             }
01376           else
01377             {
01378               dimensions(0) = new_n;
01379               dimensions(1) = 1;
01380             }
01381         }
01382       else
01383         (*current_liboctave_error_handler)
01384           ("A(idx) = []: index out of range");
01385     }
01386 }
01387 
01388 template <class T>
01389 void
01390 Array<T>::maybe_delete_elements (idx_vector& idx_i, idx_vector& idx_j)
01391 {
01392   assert (ndims () == 2);
01393 
01394   int nr = dim1 ();
01395   int nc = dim2 ();
01396 
01397   if (nr == 0 && nc == 0)
01398     return;
01399 
01400   if (idx_i.is_colon ())
01401     {
01402       if (idx_j.is_colon ())
01403         {
01404           // A(:,:) -- We are deleting columns and rows, so the result
01405           // is [](0x0).
01406 
01407           resize_no_fill (0, 0);
01408           return;
01409         }
01410 
01411       if (idx_j.is_colon_equiv (nc, 1))
01412         {
01413           // A(:,j) -- We are deleting columns by enumerating them,
01414           // If we enumerate all of them, we should have zero columns
01415           // with the same number of rows that we started with.
01416 
01417           resize_no_fill (nr, 0);
01418           return;
01419         }
01420     }
01421 
01422   if (idx_j.is_colon () && idx_i.is_colon_equiv (nr, 1))
01423     {
01424       // A(i,:) -- We are deleting rows by enumerating them.  If we
01425       // enumerate all of them, we should have zero rows with the
01426       // same number of columns that we started with.
01427 
01428       resize_no_fill (0, nc);
01429       return;
01430     }
01431 
01432   if (idx_i.is_colon_equiv (nr, 1))
01433     {
01434       if (idx_j.is_colon_equiv (nc, 1))
01435         resize_no_fill (0, 0);
01436       else
01437         {
01438           idx_j.sort (true);
01439 
01440           int num_to_delete = idx_j.length (nc);
01441 
01442           if (num_to_delete != 0)
01443             {
01444               if (nr == 1 && num_to_delete == nc)
01445                 resize_no_fill (0, 0);
01446               else
01447                 {
01448                   int new_nc = nc;
01449 
01450                   int iidx = 0;
01451 
01452                   for (int j = 0; j < nc; j++)
01453                     if (j == idx_j.elem (iidx))
01454                       {
01455                         iidx++;
01456                         new_nc--;
01457 
01458                         if (iidx == num_to_delete)
01459                           break;
01460                       }
01461 
01462                   if (new_nc > 0)
01463                     {
01464                       T *new_data = new T [nr * new_nc];
01465 
01466                       int jj = 0;
01467                       iidx = 0;
01468                       for (int j = 0; j < nc; j++)
01469                         {
01470                           if (iidx < num_to_delete && j == idx_j.elem (iidx))
01471                             iidx++;
01472                           else
01473                             {
01474                               for (int i = 0; i < nr; i++)
01475                                 new_data[nr*jj+i] = elem (i, j);
01476                               jj++;
01477                             }
01478                         }
01479 
01480                       if (--(Array<T>::rep)->count <= 0)
01481                         delete Array<T>::rep;
01482 
01483                       Array<T>::rep = new typename Array<T>::ArrayRep (new_data, nr * new_nc);
01484 
01485                       dimensions.resize (2);
01486                       dimensions(1) = new_nc;
01487                     }
01488                   else
01489                     (*current_liboctave_error_handler)
01490                       ("A(idx) = []: index out of range");
01491                 }
01492             }
01493         }
01494     }
01495   else if (idx_j.is_colon_equiv (nc, 1))
01496     {
01497       if (idx_i.is_colon_equiv (nr, 1))
01498         resize_no_fill (0, 0);
01499       else
01500         {
01501           idx_i.sort (true);
01502 
01503           int num_to_delete = idx_i.length (nr);
01504 
01505           if (num_to_delete != 0)
01506             {
01507               if (nc == 1 && num_to_delete == nr)
01508                 resize_no_fill (0, 0);
01509               else
01510                 {
01511                   int new_nr = nr;
01512 
01513                   int iidx = 0;
01514 
01515                   for (int i = 0; i < nr; i++)
01516                     if (i == idx_i.elem (iidx))
01517                       {
01518                         iidx++;
01519                         new_nr--;
01520 
01521                         if (iidx == num_to_delete)
01522                           break;
01523                       }
01524 
01525                   if (new_nr > 0)
01526                     {
01527                       T *new_data = new T [new_nr * nc];
01528 
01529                       int ii = 0;
01530                       iidx = 0;
01531                       for (int i = 0; i < nr; i++)
01532                         {
01533                           if (iidx < num_to_delete && i == idx_i.elem (iidx))
01534                             iidx++;
01535                           else
01536                             {
01537                               for (int j = 0; j < nc; j++)
01538                                 new_data[new_nr*j+ii] = elem (i, j);
01539                               ii++;
01540                             }
01541                         }
01542 
01543                       if (--(Array<T>::rep)->count <= 0)
01544                         delete Array<T>::rep;
01545 
01546                       Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_nr * nc);
01547 
01548                       dimensions.resize (2);
01549                       dimensions(0) = new_nr;
01550                     }
01551                   else
01552                     (*current_liboctave_error_handler)
01553                       ("A(idx) = []: index out of range");
01554                 }
01555             }
01556         }
01557     }
01558 }
01559 
01560 template <class T>
01561 void
01562 Array<T>::maybe_delete_elements (idx_vector&, idx_vector&, idx_vector&)
01563 {
01564   assert (0);
01565 }
01566 
01567 template <class T>
01568 void
01569 Array<T>::maybe_delete_elements (Array<idx_vector>& ra_idx, const T& rfv)
01570 {
01571   int n_idx = ra_idx.length ();
01572 
01573   dim_vector lhs_dims = dims ();
01574 
01575   if (lhs_dims.all_zero ())
01576     return;
01577 
01578   int n_lhs_dims = lhs_dims.length ();
01579 
01580   Array<int> idx_is_colon (n_idx, 0);
01581 
01582   Array<int> idx_is_colon_equiv (n_idx, 0);
01583 
01584   // Initialization of colon arrays.
01585 
01586   for (int i = 0; i < n_idx; i++)
01587     {
01588       idx_is_colon_equiv(i) = ra_idx(i).is_colon_equiv (lhs_dims(i), 1);
01589 
01590       idx_is_colon(i) = ra_idx(i).is_colon ();
01591     }
01592 
01593   bool idx_ok = true;
01594 
01595   // Check for index out of bounds.
01596 
01597   for (int i = 0 ; i < n_idx - 1; i++)
01598     {
01599       if (! (idx_is_colon(i) || idx_is_colon_equiv(i)))
01600         {
01601           ra_idx(i).sort (true);
01602 
01603           if (ra_idx(i).max () > lhs_dims(i))
01604             {
01605               (*current_liboctave_error_handler)
01606                 ("index exceeds array dimensions");
01607 
01608               idx_ok = false;
01609               break;
01610             }
01611           else if (ra_idx(i).min () < 0) // I believe this is checked elsewhere
01612             {
01613               (*current_liboctave_error_handler)
01614                 ("index must be one or larger");
01615 
01616               idx_ok = false;
01617               break;
01618             }
01619         }
01620     }
01621 
01622   if (n_idx <= n_lhs_dims)
01623     {
01624       int last_idx = ra_idx(n_idx-1).max ();
01625 
01626       int sum_el = lhs_dims(n_idx-1);
01627 
01628       for (int i = n_idx; i < n_lhs_dims; i++)
01629           sum_el *= lhs_dims(i);
01630 
01631       if (last_idx > sum_el - 1)
01632         {
01633           (*current_liboctave_error_handler)
01634             ("index exceeds array dimensions");
01635 
01636           idx_ok = false;
01637         }
01638     }
01639 
01640   if (idx_ok)
01641     {
01642       if (n_idx > 1
01643           && (all_ones (idx_is_colon) || all_ones (idx_is_colon_equiv)))
01644         {
01645           // A(:,:,:) -- we are deleting elements in all dimensions, so
01646           // the result is [](0x0x0).
01647 
01648           dim_vector zeros;
01649           zeros.resize (n_idx);
01650 
01651           for (int i = 0; i < n_idx; i++)
01652             zeros(i) = 0;
01653 
01654           resize (zeros, rfv);
01655         }
01656 
01657       else if (n_idx > 1
01658                && num_ones (idx_is_colon) == n_idx - 1
01659                && num_ones (idx_is_colon_equiv) == n_idx)
01660         {
01661           // A(:,:,j) -- we are deleting elements in one dimension by
01662           // enumerating them.
01663           //
01664           // If we enumerate all of the elements, we should have zero
01665           // elements in that dimension with the same number of elements
01666           // in the other dimensions that we started with.
01667 
01668           dim_vector temp_dims;
01669           temp_dims.resize (n_idx);
01670 
01671           for (int i = 0; i < n_idx; i++)
01672             {
01673               if (idx_is_colon (i))
01674                 temp_dims(i) =  lhs_dims(i);
01675               else
01676                 temp_dims(i) = 0;
01677             }
01678 
01679           resize (temp_dims);
01680         }
01681       else if (n_idx > 1 && num_ones (idx_is_colon) == n_idx - 1)
01682         {
01683           // We have colons in all indices except for one.
01684           // This index tells us which slice to delete
01685 
01686           if (n_idx < n_lhs_dims)
01687             {
01688               // Collapse dimensions beyond last index.
01689 
01690               if (liboctave_wfi_flag && ! (ra_idx(n_idx-1).is_colon ()))
01691                 (*current_liboctave_warning_handler)
01692                   ("fewer indices than dimensions for N-d array");
01693 
01694               for (int i = n_idx; i < n_lhs_dims; i++)
01695                 lhs_dims(n_idx-1) *= lhs_dims(i);
01696 
01697               lhs_dims.resize (n_idx);
01698 
01699               // Reshape *this.
01700               dimensions = lhs_dims;
01701             }
01702 
01703           int non_col = 0;
01704 
01705           // Find the non-colon column.
01706 
01707           for (int i = 0; i < n_idx; i++)
01708             {
01709               if (! idx_is_colon(i))
01710                 non_col = i;
01711             }
01712 
01713           // The length of the non-colon dimension.
01714 
01715           int non_col_dim = lhs_dims (non_col);
01716 
01717           int num_to_delete = ra_idx(non_col).length (lhs_dims (non_col));
01718 
01719           if (num_to_delete > 0)
01720             {
01721               int temp = lhs_dims.num_ones ();
01722 
01723               if (non_col_dim == 1)
01724                 temp--;
01725 
01726               if (temp == n_idx - 1 && num_to_delete == non_col_dim)
01727                 {
01728                   // We have A with (1x1x4), where A(1,:,1:4)
01729                   // Delete all (0x0x0)
01730 
01731                   dim_vector zero_dims (n_idx, 0);
01732 
01733                   resize (zero_dims, rfv);
01734                 }
01735               else
01736                 {
01737                   // New length of non-colon dimension
01738                   // (calculated in the next for loop)
01739 
01740                   int new_dim = non_col_dim;
01741 
01742                   int iidx = 0;
01743 
01744                   for (int j = 0; j < non_col_dim; j++)
01745                     if (j == ra_idx(non_col).elem (iidx))
01746                       {
01747                         iidx++;
01748 
01749                         new_dim--;
01750 
01751                         if (iidx == num_to_delete)
01752                           break;
01753                       }
01754 
01755                   // Creating the new nd array after deletions.
01756 
01757                   if (new_dim > 0)
01758                     {
01759                       // Calculate number of elements in new array.
01760 
01761                       int num_new_elem=1;
01762 
01763                       for (int i = 0; i < n_idx; i++)
01764                         {
01765                           if (i == non_col)
01766                             num_new_elem *= new_dim;
01767 
01768                           else
01769                             num_new_elem *= lhs_dims(i);
01770                         }
01771 
01772                       T *new_data = new T [num_new_elem];
01773 
01774                       Array<int> result_idx (n_lhs_dims, 0);
01775 
01776                       dim_vector new_lhs_dim = lhs_dims;
01777 
01778                       new_lhs_dim(non_col) = new_dim;
01779 
01780                       int num_elem = 1;
01781 
01782                       int numidx = 0;
01783 
01784                       int n = length ();
01785 
01786                       for (int i = 0; i < n_lhs_dims; i++)
01787                         if (i != non_col)
01788                           num_elem *= lhs_dims(i);
01789 
01790                       num_elem *= ra_idx(non_col).capacity ();
01791 
01792                       for (int i = 0; i < n; i++)
01793                         {
01794                           if (numidx < num_elem
01795                               && is_in (result_idx(non_col), ra_idx(non_col)))
01796                             numidx++;
01797 
01798                           else
01799                             {
01800                               Array<int> temp_result_idx = result_idx;
01801 
01802                               int num_lgt = how_many_lgt (result_idx(non_col),
01803                                                           ra_idx(non_col));
01804 
01805                               temp_result_idx(non_col) -= num_lgt;
01806 
01807                               int kidx
01808                                 = ::compute_index (temp_result_idx, new_lhs_dim);
01809 
01810                               new_data[kidx] = elem (result_idx);
01811                             }
01812 
01813                           increment_index (result_idx, lhs_dims);
01814                         }
01815 
01816                       if (--rep->count <= 0)
01817                         delete rep;
01818 
01819                       rep = new typename Array<T>::ArrayRep (new_data,
01820                                                              num_new_elem);
01821 
01822                       dimensions = new_lhs_dim;
01823                     }
01824                 }
01825             }
01826         }
01827       else if (n_idx == 1)
01828         {
01829           // This handle cases where we only have one index (not
01830           // colon).  The index denotes which elements we should
01831           // delete in the array which can be of any dimension. We
01832           // return a column vector, except for the case where we are
01833           // operating on a row vector. The elements are numerated
01834           // column by column.
01835           //
01836           // A(3,3,3)=2;
01837           // A(3:5) = []; A(6)=[]
01838 
01839           int lhs_numel = numel ();
01840 
01841           idx_vector idx_vec = ra_idx(0);
01842 
01843           idx_vec.freeze (lhs_numel, 0, true, liboctave_wrore_flag);
01844       
01845           idx_vec.sort (true);
01846 
01847           int num_to_delete = idx_vec.length (lhs_numel);
01848 
01849           if (num_to_delete > 0)
01850             {
01851               int new_numel = lhs_numel - num_to_delete;
01852 
01853               T *new_data = new T[new_numel];
01854 
01855               Array<int> lhs_ra_idx (ndims (), 0);
01856 
01857               int ii = 0;
01858               int iidx = 0;
01859 
01860               for (int i = 0; i < lhs_numel; i++)
01861                 {
01862                   if (iidx < num_to_delete && i == idx_vec.elem (iidx))
01863                     {
01864                       iidx++;
01865                     }
01866                   else
01867                     {
01868                       new_data[ii++] = elem (lhs_ra_idx);
01869                     }
01870 
01871                   increment_index (lhs_ra_idx, lhs_dims);
01872                 }
01873 
01874               if (--(Array<T>::rep)->count <= 0)
01875                 delete Array<T>::rep;
01876 
01877               Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_numel);
01878 
01879               dimensions.resize (2);
01880 
01881               if (lhs_dims.length () == 2 && lhs_dims(1) == 1)
01882                 {
01883                   dimensions(0) = new_numel;
01884                   dimensions(1) = 1;
01885                 }
01886               else
01887                 {
01888                   dimensions(0) = 1;
01889                   dimensions(1) = new_numel;
01890                 }
01891             }
01892         }
01893       else if (num_ones (idx_is_colon) < n_idx)
01894         {
01895           (*current_liboctave_error_handler)
01896             ("a null assignment can have only one non-colon index");
01897         }
01898     }
01899 }
01900 
01901 template <class T>
01902 Array<T>
01903 Array<T>::value (void)
01904 {
01905   Array<T> retval;
01906 
01907   int n_idx = index_count ();
01908 
01909   if (n_idx == 2)
01910     {
01911       idx_vector *tmp = get_idx ();
01912 
01913       idx_vector idx_i = tmp[0];
01914       idx_vector idx_j = tmp[1];
01915 
01916       retval = index (idx_i, idx_j);
01917     }
01918   else if (n_idx == 1)
01919     {
01920       retval = index (idx[0]);
01921     }
01922   else
01923     (*current_liboctave_error_handler)
01924       ("Array<T>::value: invalid number of indices specified");
01925 
01926   clear_index ();
01927 
01928   return retval;
01929 }
01930 
01931 template <class T>
01932 Array<T>
01933 Array<T>::index (idx_vector& idx_arg, int resize_ok, const T& rfv) const
01934 {
01935   Array<T> retval;
01936 
01937   dim_vector dv = idx_arg.orig_dimensions ();
01938 
01939   if (dv.length () > 2 || ndims () > 2)
01940     retval = indexN (idx_arg, resize_ok, rfv);
01941   else
01942     {
01943       switch (ndims ())
01944         {
01945         case 1:
01946           retval = index1 (idx_arg, resize_ok, rfv);
01947           break;
01948 
01949         case 2:
01950           retval = index2 (idx_arg, resize_ok, rfv);
01951           break;
01952 
01953         default:
01954           (*current_liboctave_error_handler)
01955             ("invalid array (internal error)");
01956           break;
01957         }
01958     }
01959 
01960   return retval;
01961 }
01962 
01963 template <class T>
01964 Array<T>
01965 Array<T>::index1 (idx_vector& idx_arg, int resize_ok, const T& rfv) const
01966 {
01967   Array<T> retval;
01968 
01969   int len = length ();
01970 
01971   int n = idx_arg.freeze (len, "vector", resize_ok);
01972 
01973   if (idx_arg)
01974     {
01975       if (idx_arg.is_colon_equiv (len))
01976         {
01977           retval = *this;
01978         }
01979       else if (n == 0)
01980         {
01981           retval.resize_no_fill (0);
01982         }
01983       else if (len == 1 && n > 1
01984                && idx_arg.one_zero_only ()
01985                && idx_arg.ones_count () == n)
01986         {
01987           retval.resize_and_fill (n, elem (0));
01988         }
01989       else
01990         {
01991           retval.resize_no_fill (n);
01992 
01993           for (int i = 0; i < n; i++)
01994             {
01995               int ii = idx_arg.elem (i);
01996               if (ii >= len)
01997                 retval.elem (i) = rfv;
01998               else
01999                 retval.elem (i) = elem (ii);
02000             }
02001         }
02002     }
02003 
02004   // idx_vector::freeze() printed an error message for us.
02005 
02006   return retval;
02007 }
02008 
02009 template <class T>
02010 Array<T>
02011 Array<T>::index2 (idx_vector& idx_arg, int resize_ok, const T& rfv) const
02012 {
02013   Array<T> retval;
02014 
02015   assert (ndims () == 2);
02016 
02017   int nr = dim1 ();
02018   int nc = dim2 ();
02019 
02020   int orig_len = nr * nc;
02021 
02022   dim_vector idx_orig_dims = idx_arg.orig_dimensions ();
02023 
02024   int idx_orig_rows = idx_arg.orig_rows ();
02025   int idx_orig_columns = idx_arg.orig_columns ();
02026 
02027   if (idx_arg.is_colon ())
02028     {
02029       // Fast magic colon processing.
02030 
02031       int result_nr = nr * nc;
02032       int result_nc = 1;
02033 
02034       retval = Array<T> (*this, dim_vector (result_nr, result_nc));
02035     }
02036   else if (nr == 1 && nc == 1)
02037     {
02038       Array<T> tmp = Array<T>::index1 (idx_arg, resize_ok);
02039 
02040       int len = tmp.length ();
02041 
02042       if (len == 0 && idx_arg.one_zero_only ())
02043         retval = Array<T> (tmp, dim_vector (0, 0));
02044       else if (len >= idx_orig_dims.numel ())
02045         retval = Array<T> (tmp, idx_orig_dims);
02046     }
02047   else if (nr == 1 || nc == 1)
02048     {
02049       // If indexing a vector with a matrix, return value has same
02050       // shape as the index.  Otherwise, it has same orientation as
02051       // indexed object.
02052 
02053       Array<T> tmp = Array<T>::index1 (idx_arg, resize_ok);
02054 
02055       int len = tmp.length ();
02056 
02057       if ((len != 0 && idx_arg.one_zero_only ())
02058           || idx_orig_rows == 1 || idx_orig_columns == 1)
02059         {
02060           if (nr == 1)
02061             retval = Array<T> (tmp, dim_vector (1, len));
02062           else
02063             retval = Array<T> (tmp, dim_vector (len, 1));
02064         }
02065       else if (len >= idx_orig_dims.numel ())
02066         retval = Array<T> (tmp, idx_orig_dims);
02067     }
02068   else
02069     {
02070       if (liboctave_wfi_flag
02071           && ! (idx_arg.one_zero_only ()
02072                 && idx_orig_rows == nr
02073                 && idx_orig_columns == nc))
02074         (*current_liboctave_warning_handler) ("single index used for matrix");
02075 
02076       // This code is only for indexing matrices.  The vector
02077       // cases are handled above.
02078 
02079       idx_arg.freeze (nr * nc, "matrix", resize_ok);
02080 
02081       if (idx_arg)
02082         {
02083           int result_nr = idx_orig_rows;
02084           int result_nc = idx_orig_columns;
02085 
02086           if (idx_arg.one_zero_only ())
02087             {
02088               result_nr = idx_arg.ones_count ();
02089               result_nc = (result_nr > 0 ? 1 : 0);
02090             }
02091 
02092           retval.resize_no_fill (result_nr, result_nc);
02093 
02094           int k = 0;
02095           for (int j = 0; j < result_nc; j++)
02096             {
02097               for (int i = 0; i < result_nr; i++)
02098                 {
02099                   int ii = idx_arg.elem (k++);
02100                   if (ii >= orig_len)
02101                     retval.elem (i, j) = rfv;
02102                   else
02103                     {
02104                       int fr = ii % nr;
02105                       int fc = (ii - fr) / nr;
02106                       retval.elem (i, j) = elem (fr, fc);
02107                     }
02108                 }
02109             }
02110         }
02111       // idx_vector::freeze() printed an error message for us.
02112     }
02113 
02114   return retval;
02115 }
02116 
02117 template <class T>
02118 Array<T>
02119 Array<T>::indexN (idx_vector& ra_idx, int resize_ok, const T& rfv) const
02120 {
02121   Array<T> retval;
02122 
02123   int n_dims = dims().length ();
02124 
02125   int orig_len = dims().numel ();
02126 
02127   dim_vector idx_orig_dims = ra_idx.orig_dimensions ();
02128 
02129   if (ra_idx.is_colon ())
02130     {
02131       // Fast magic colon processing.
02132 
02133       retval = Array<T> (*this, dim_vector (orig_len, 1));
02134     }
02135   else if (length () == 1)
02136     {
02137       // Only one element in array.
02138 
02139       Array<T> tmp = Array<T>::index (ra_idx, resize_ok);
02140 
02141       int len = tmp.length ();
02142 
02143       if (len != 0)
02144         {
02145           if (len >= idx_orig_dims.numel ())
02146             retval = Array<T> (tmp, idx_orig_dims);
02147         }
02148       else
02149         retval = Array<T> (tmp, dim_vector (0, 0));
02150     }
02151   else if (vector_equivalent (dims ()))
02152     {
02153       // We're getting elements from a vector equivalent i.e. (1x4x1).
02154 
02155       Array<T> tmp = Array<T>::index (ra_idx, resize_ok);
02156 
02157       int len = tmp.length ();
02158 
02159       if (len == 0)
02160         {
02161           if (idx_orig_dims.any_zero ())
02162             retval = Array<T> (idx_orig_dims);
02163           else
02164             {
02165               dim_vector new_dims;
02166 
02167               new_dims.resize (n_dims);
02168 
02169               for (int i = 0; i < n_dims; i++)
02170                 {
02171                   if ((dims ())(i) == 1)
02172                     new_dims(i) = 1;
02173                 }
02174 
02175               new_dims.chop_trailing_singletons ();
02176 
02177               retval = Array<T> (new_dims);
02178             }
02179         }
02180       else
02181         {
02182           if (vector_equivalent (idx_orig_dims))
02183             {
02184               // Array<int> index (n_dims, len);
02185               dim_vector new_dims;
02186 
02187               new_dims.resize (n_dims);
02188 
02189               for (int i = 0; i < n_dims; i++)
02190                 {
02191                   if ((dims ())(i) == 1)
02192                     new_dims(i) = 1;
02193                 }
02194 
02195               new_dims.chop_trailing_singletons ();
02196 
02197               retval = Array<T> (tmp, new_dims);
02198             }
02199           else if (tmp.length () >= idx_orig_dims.numel ())
02200             retval = Array<T> (tmp, idx_orig_dims);
02201 
02202           (*current_liboctave_error_handler)
02203             ("I do not know what to do here yet!");
02204         }
02205     }
02206   else
02207     {
02208       if (liboctave_wfi_flag
02209           && ! (ra_idx.is_colon ()
02210                 || (ra_idx.one_zero_only () && idx_orig_dims == dims ())))
02211         (*current_liboctave_warning_handler)
02212           ("single index used for N-d array");
02213 
02214       ra_idx.freeze (orig_len, "nd-array", resize_ok);
02215 
02216       if (ra_idx)
02217         {
02218           dim_vector result_dims (idx_orig_dims);
02219 
02220           if (ra_idx.one_zero_only ())
02221             {
02222               result_dims.resize (2);
02223               int ntot = ra_idx.ones_count ();
02224               result_dims(0) = ntot;
02225               result_dims(1) = (ntot > 0 ? 1 : 0);
02226             }
02227 
02228           result_dims.chop_trailing_singletons ();
02229 
02230           retval.resize (result_dims);
02231 
02232           int n = result_dims.numel ();
02233 
02234           int r_dims = result_dims.length ();
02235 
02236           Array<int> iidx (r_dims, 0);
02237 
02238           int k = 0;
02239 
02240           for (int i = 0; i < n; i++)
02241             {
02242               int ii = ra_idx.elem (k++);
02243 
02244               if (ii >= orig_len)
02245                 retval.elem (iidx) = rfv;
02246               else
02247                 {
02248                   Array<int> temp = get_ra_idx (ii, dims ());
02249 
02250                   retval.elem (iidx) = elem (temp);
02251                 }
02252               if (i != n - 1)
02253                 increment_index (iidx, result_dims);
02254             }
02255         }
02256     }
02257 
02258   return retval;
02259 }
02260 
02261 template <class T>
02262 Array<T>
02263 Array<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok,
02264                  const T& rfv) const
02265 {
02266   Array<T> retval;
02267 
02268   assert (ndims () == 2);
02269 
02270   int nr = dim1 ();
02271   int nc = dim2 ();
02272 
02273   int n = idx_i.freeze (nr, "row", resize_ok);
02274   int m = idx_j.freeze (nc, "column", resize_ok);
02275 
02276   if (idx_i && idx_j)
02277     {
02278       if (idx_i.orig_empty () || idx_j.orig_empty () || n == 0 || m == 0)
02279         {
02280           retval.resize_no_fill (n, m);
02281         }
02282       else if (idx_i.is_colon_equiv (nr) && idx_j.is_colon_equiv (nc))
02283         {
02284           retval = *this;
02285         }
02286       else
02287         {
02288           retval.resize_no_fill (n, m);
02289 
02290           for (int j = 0; j < m; j++)
02291             {
02292               int jj = idx_j.elem (j);
02293               for (int i = 0; i < n; i++)
02294                 {
02295                   int ii = idx_i.elem (i);
02296                   if (ii >= nr || jj >= nc)
02297                     retval.elem (i, j) = rfv;
02298                   else
02299                     retval.elem (i, j) = elem (ii, jj);
02300                 }
02301             }
02302         }
02303     }
02304 
02305   // idx_vector::freeze() printed an error message for us.
02306 
02307   return retval;
02308 }
02309 
02310 template <class T>
02311 Array<T>
02312 Array<T>::index (Array<idx_vector>& ra_idx, int resize_ok, const T&) const
02313 {
02314   // This function handles all calls with more than one idx.
02315   // For (3x3x3), the call can be A(2,5), A(2,:,:), A(3,2,3) etc.
02316 
02317   Array<T> retval;
02318 
02319   int n_dims = dimensions.length ();
02320 
02321   // Remove trailing singletons in ra_idx, but leave at least ndims
02322   // elements.
02323 
02324   int ra_idx_len = ra_idx.length ();
02325 
02326   bool trim_trailing_singletons = true;
02327   for (int j = ra_idx_len; j > n_dims; j--)
02328     {
02329       idx_vector iidx = ra_idx (ra_idx_len-1);
02330       if (iidx.capacity () == 1 && trim_trailing_singletons)
02331         ra_idx_len--;
02332       else
02333         trim_trailing_singletons = false;
02334 
02335       for (int i = 0; i < iidx.capacity (); i++)
02336         if (iidx (i) != 0)
02337           {
02338             (*current_liboctave_error_handler)
02339               ("index exceeds N-d array dimensions");
02340             
02341             return retval;
02342           }
02343     }
02344 
02345   ra_idx.resize (ra_idx_len);
02346 
02347   dim_vector new_dims = dims ();
02348   dim_vector frozen_lengths;
02349 
02350   if (! any_orig_empty (ra_idx) && ra_idx_len < n_dims)
02351     frozen_lengths = short_freeze (ra_idx, dimensions, resize_ok);
02352   else
02353     {
02354       new_dims.resize (ra_idx_len, 1);
02355       frozen_lengths = freeze (ra_idx, new_dims, resize_ok);
02356     }
02357 
02358   if (all_ok (ra_idx))
02359     {
02360       if (any_orig_empty (ra_idx) || frozen_lengths.any_zero ())
02361         {
02362           frozen_lengths.chop_trailing_singletons ();
02363 
02364           retval.resize (frozen_lengths);
02365         }
02366       else if (frozen_lengths.length () == n_dims
02367                && all_colon_equiv (ra_idx, dimensions))
02368         {
02369           retval = *this;
02370         }
02371       else
02372         {
02373           dim_vector frozen_lengths_for_resize = frozen_lengths;
02374 
02375           frozen_lengths_for_resize.chop_trailing_singletons ();
02376 
02377           retval.resize (frozen_lengths_for_resize);
02378 
02379           int n = retval.length ();
02380 
02381           Array<int> result_idx (ra_idx.length (), 0);
02382 
02383           Array<int> elt_idx;
02384 
02385           for (int i = 0; i < n; i++)
02386             {
02387               elt_idx = get_elt_idx (ra_idx, result_idx);
02388 
02389               int numelem_elt = get_scalar_idx (elt_idx, new_dims);
02390 
02391               if (numelem_elt > length () || numelem_elt < 0)
02392                 (*current_liboctave_error_handler)
02393                   ("invalid N-d array index");
02394               else
02395                 retval.elem (i) = elem (numelem_elt);
02396 
02397               increment_index (result_idx, frozen_lengths);
02398 
02399             }
02400         }
02401     }
02402 
02403   return retval;
02404 }
02405 
02406 // XXX FIXME XXX -- this is a mess.
02407 
02408 template <class LT, class RT>
02409 int
02410 assign (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
02411 {
02412   int retval = 0;
02413 
02414   switch (lhs.ndims ())
02415     {
02416     case 0:
02417       {
02418         if (lhs.index_count () < 3)
02419           {
02420             // kluge...
02421             lhs.resize_no_fill (0, 0);
02422             retval = assign2 (lhs, rhs, rfv);
02423           }
02424         else
02425           retval = assignN (lhs, rhs, rfv);
02426       }
02427       break;
02428 
02429     case 1:
02430       {
02431         if (lhs.index_count () > 1)
02432           retval = assignN (lhs, rhs, rfv);
02433         else
02434           retval = assign1 (lhs, rhs, rfv);
02435       }
02436       break;
02437 
02438     case 2:
02439       {
02440         if (lhs.index_count () > 2)
02441           retval = assignN (lhs, rhs, rfv);
02442         else
02443           retval = assign2 (lhs, rhs, rfv);
02444       }
02445       break;
02446 
02447     default:
02448       retval = assignN (lhs, rhs, rfv);
02449       break;
02450     }
02451 
02452   return retval;
02453 }
02454 
02455 template <class LT, class RT>
02456 int
02457 assign1 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
02458 {
02459   int retval = 1;
02460 
02461   idx_vector *tmp = lhs.get_idx ();
02462 
02463   idx_vector lhs_idx = tmp[0];
02464 
02465   int lhs_len = lhs.length ();
02466   int rhs_len = rhs.length ();
02467 
02468   int n = lhs_idx.freeze (lhs_len, "vector", true, liboctave_wrore_flag);
02469 
02470   if (n != 0)
02471     {
02472       if (rhs_len == n || rhs_len == 1)
02473         {
02474           int max_idx = lhs_idx.max () + 1;
02475           if (max_idx > lhs_len)
02476             lhs.resize_and_fill (max_idx, rfv);
02477         }
02478 
02479       if (rhs_len == n)
02480         {
02481           for (int i = 0; i < n; i++)
02482             {
02483               int ii = lhs_idx.elem (i);
02484               lhs.elem (ii) = rhs.elem (i);
02485             }
02486         }
02487       else if (rhs_len == 1)
02488         {
02489           RT scalar = rhs.elem (0);
02490 
02491           for (int i = 0; i < n; i++)
02492             {
02493               int ii = lhs_idx.elem (i);
02494               lhs.elem (ii) = scalar;
02495             }
02496         }
02497       else
02498         {
02499           (*current_liboctave_error_handler)
02500             ("A(I) = X: X must be a scalar or a vector with same length as I");
02501 
02502           retval = 0;
02503         }
02504     }
02505   else if (lhs_idx.is_colon ())
02506     {
02507       if (lhs_len == 0)
02508         {
02509           lhs.resize_no_fill (rhs_len);
02510 
02511           for (int i = 0; i < rhs_len; i++)
02512             lhs.elem (i) = rhs.elem (i);
02513         }
02514       else
02515         (*current_liboctave_error_handler)
02516           ("A(:) = X: A must be the same size as X");
02517     }
02518   else if (! (rhs_len == 1 || rhs_len == 0))
02519     {
02520       (*current_liboctave_error_handler)
02521         ("A([]) = X: X must also be an empty matrix or a scalar");
02522 
02523       retval = 0;
02524     }
02525 
02526   lhs.clear_index ();
02527 
02528   return retval;
02529 }
02530 
02531 #define MAYBE_RESIZE_LHS \
02532   do \
02533     { \
02534       int max_row_idx = idx_i_is_colon ? rhs_nr : idx_i.max () + 1; \
02535       int max_col_idx = idx_j_is_colon ? rhs_nc : idx_j.max () + 1; \
02536  \
02537       int new_nr = max_row_idx > lhs_nr ? max_row_idx : lhs_nr; \
02538       int new_nc = max_col_idx > lhs_nc ? max_col_idx : lhs_nc; \
02539  \
02540       lhs.resize_and_fill (new_nr, new_nc, rfv); \
02541     } \
02542   while (0)
02543 
02544 template <class LT, class RT>
02545 int
02546 assign2 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
02547 {
02548   int retval = 1;
02549 
02550   int n_idx = lhs.index_count ();
02551 
02552   int lhs_nr = lhs.rows ();
02553   int lhs_nc = lhs.cols ();
02554 
02555   Array<RT> xrhs = rhs;
02556 
02557   int rhs_nr = xrhs.rows ();
02558   int rhs_nc = xrhs.cols ();
02559 
02560   if (xrhs.ndims () > 2)
02561     {
02562       xrhs = xrhs.squeeze ();
02563 
02564       dim_vector dv_tmp = xrhs.dims ();
02565 
02566       switch (dv_tmp.length ())
02567         {
02568         case 1:
02569           // XXX FIXME XXX -- this case should be unnecessary, because
02570           // squeeze should always return an object with 2 dimensions.
02571           if (rhs_nr == 1)
02572             rhs_nc = dv_tmp.elem (0);
02573           break;
02574 
02575         case 2:
02576           rhs_nr = dv_tmp.elem (0);
02577           rhs_nc = dv_tmp.elem (1);
02578           break;
02579 
02580         default:
02581           (*current_liboctave_error_handler)
02582             ("Array<T>::assign2: Dimension mismatch");
02583           return 0;
02584         }
02585     }
02586 
02587   idx_vector *tmp = lhs.get_idx ();
02588 
02589   idx_vector idx_i;
02590   idx_vector idx_j;
02591 
02592   if (n_idx > 1)
02593     idx_j = tmp[1];
02594 
02595   if (n_idx > 0)
02596     idx_i = tmp[0];
02597 
02598   if (n_idx == 2)
02599     {
02600       int n = idx_i.freeze (lhs_nr, "row", true, liboctave_wrore_flag);
02601 
02602       int m = idx_j.freeze (lhs_nc, "column", true, liboctave_wrore_flag);
02603 
02604       int idx_i_is_colon = idx_i.is_colon ();
02605       int idx_j_is_colon = idx_j.is_colon ();
02606 
02607       if (idx_i_is_colon)
02608         n = lhs_nr > 0 ? lhs_nr : rhs_nr;
02609 
02610       if (idx_j_is_colon)
02611         m = lhs_nc > 0 ? lhs_nc : rhs_nc;
02612 
02613       if (idx_i && idx_j)
02614         {
02615           if (rhs_nr == 0 && rhs_nc == 0)
02616             {
02617               lhs.maybe_delete_elements (idx_i, idx_j);
02618             }
02619           else
02620             {
02621               if (rhs_nr == 1 && rhs_nc == 1 && n >= 0 && m >= 0)
02622                 {
02623                   // No need to do anything if either of the indices
02624                   // are empty.
02625 
02626                   if (n > 0 && m > 0)
02627                     {
02628                       MAYBE_RESIZE_LHS;
02629 
02630                       RT scalar = xrhs.elem (0, 0);
02631 
02632                       for (int j = 0; j < m; j++)
02633                         {
02634                           int jj = idx_j.elem (j);
02635                           for (int i = 0; i < n; i++)
02636                             {
02637                               int ii = idx_i.elem (i);
02638                               lhs.elem (ii, jj) = scalar;
02639                             }
02640                         }
02641                     }
02642                 }
02643               else if (n == rhs_nr && m == rhs_nc)
02644                 {
02645                   if (n > 0 && m > 0)
02646                     {
02647                       MAYBE_RESIZE_LHS;
02648 
02649                       for (int j = 0; j < m; j++)
02650                         {
02651                           int jj = idx_j.elem (j);
02652                           for (int i = 0; i < n; i++)
02653                             {
02654                               int ii = idx_i.elem (i);
02655                               lhs.elem (ii, jj) = xrhs.elem (i, j);
02656                             }
02657                         }
02658                     }
02659                 }
02660               else if (n == 0 && m == 0)
02661                 {
02662                   if (! ((rhs_nr == 1 && rhs_nc == 1)
02663                          || (rhs_nr == 0 || rhs_nc == 0)))
02664                     {
02665                       (*current_liboctave_error_handler)
02666                 ("A([], []) = X: X must be an empty matrix or a scalar");
02667 
02668                       retval = 0;
02669                     }
02670                 }
02671               else
02672                 {
02673                   (*current_liboctave_error_handler)
02674     ("A(I, J) = X: X must be a scalar or the number of elements in I must");
02675                   (*current_liboctave_error_handler)
02676     ("match the number of rows in X and the number of elements in J must");
02677                   (*current_liboctave_error_handler)
02678     ("match the number of columns in X");
02679 
02680                   retval = 0;
02681                 }
02682             }
02683         }
02684       // idx_vector::freeze() printed an error message for us.
02685     }
02686   else if (n_idx == 1)
02687     {
02688       int lhs_is_empty = lhs_nr == 0 || lhs_nc == 0;
02689 
02690       if (lhs_is_empty || (lhs_nr == 1 && lhs_nc == 1))
02691         {
02692           int lhs_len = lhs.length ();
02693 
02694           int n = idx_i.freeze (lhs_len, 0, true, liboctave_wrore_flag);
02695 
02696           if (idx_i)
02697             {
02698               if (rhs_nr == 0 && rhs_nc == 0)
02699                 {
02700                   if (n != 0 && (lhs_nr != 0 || lhs_nc != 0))
02701                     lhs.maybe_delete_elements (idx_i);
02702                 }
02703               else
02704                 {
02705                   if (liboctave_wfi_flag)
02706                     {
02707                       if (lhs_is_empty
02708                           && idx_i.is_colon ()
02709                           && ! (rhs_nr == 1 || rhs_nc == 1))
02710                         {
02711                           (*current_liboctave_warning_handler)
02712                             ("A(:) = X: X is not a vector or scalar");
02713                         }
02714                       else
02715                         {
02716                           int idx_nr = idx_i.orig_rows ();
02717                           int idx_nc = idx_i.orig_columns ();
02718 
02719                           if (! (rhs_nr == idx_nr && rhs_nc == idx_nc))
02720                             (*current_liboctave_warning_handler)
02721                               ("A(I) = X: X does not have same shape as I");
02722                         }
02723                     }
02724 
02725                   if (assign1 (lhs, xrhs, rfv))
02726                     {
02727                       int len = lhs.length ();
02728 
02729                       if (len > 0)
02730                         {
02731                           // The following behavior is much simplified
02732                           // over previous versions of Octave.  It
02733                           // seems to be compatible with Matlab.
02734 
02735                           lhs.dimensions = dim_vector (1, lhs.length ());
02736                         }
02737                       else
02738                         lhs.dimensions = dim_vector (0, 0);
02739                     }
02740                   else
02741                     retval = 0;
02742                 }
02743             }
02744           // idx_vector::freeze() printed an error message for us.
02745         }
02746       else if (lhs_nr == 1)
02747         {
02748           idx_i.freeze (lhs_nc, "vector", true, liboctave_wrore_flag);
02749 
02750           if (idx_i)
02751             {
02752               if (rhs_nr == 0 && rhs_nc == 0)
02753                 lhs.maybe_delete_elements (idx_i);
02754               else
02755                 {
02756                   if (assign1 (lhs, xrhs, rfv))
02757                     lhs.dimensions = dim_vector (1, lhs.length ());
02758                   else
02759                     retval = 0;
02760                 }
02761             }
02762           // idx_vector::freeze() printed an error message for us.
02763         }
02764       else if (lhs_nc == 1)
02765         {
02766           idx_i.freeze (lhs_nr, "vector", true, liboctave_wrore_flag);
02767 
02768           if (idx_i)
02769             {
02770               if (rhs_nr == 0 && rhs_nc == 0)
02771                 lhs.maybe_delete_elements (idx_i);
02772               else
02773                 {
02774                   if (assign1 (lhs, xrhs, rfv))
02775                     lhs.dimensions = dim_vector (lhs.length (), 1);
02776                   else
02777                     retval = 0;
02778                 }
02779             }
02780           // idx_vector::freeze() printed an error message for us.
02781         }
02782       else
02783         {
02784           if (liboctave_wfi_flag
02785               && ! (idx_i.is_colon ()
02786                     || (idx_i.one_zero_only ()
02787                         && idx_i.orig_rows () == lhs_nr
02788                         && idx_i.orig_columns () == lhs_nc)))
02789             (*current_liboctave_warning_handler)
02790               ("single index used for matrix");
02791 
02792           int len = idx_i.freeze (lhs_nr * lhs_nc, "matrix");
02793 
02794           if (idx_i)
02795             {
02796               if (rhs_nr == 0 && rhs_nc == 0)
02797                 lhs.maybe_delete_elements (idx_i);
02798               else if (len == 0)
02799                 {
02800                   if (! ((rhs_nr == 1 && rhs_nc == 1)
02801                          || (rhs_nr == 0 || rhs_nc == 0)))
02802                     (*current_liboctave_error_handler)
02803                       ("A([]) = X: X must be an empty matrix or scalar");
02804                 }
02805               else if (len == rhs_nr * rhs_nc)
02806                 {
02807                   int k = 0;
02808                   for (int j = 0; j < rhs_nc; j++)
02809                     {
02810                       for (int i = 0; i < rhs_nr; i++)
02811                         {
02812                           int ii = idx_i.elem (k++);
02813                           int fr = ii % lhs_nr;
02814                           int fc = (ii - fr) / lhs_nr;
02815                           lhs.elem (fr, fc) = xrhs.elem (i, j);
02816                         }
02817                     }
02818                 }
02819               else if (rhs_nr == 1 && rhs_nc == 1)
02820                 {
02821                   RT scalar = rhs.elem (0, 0);
02822 
02823                   for (int i = 0; i < len; i++)
02824                     {
02825                       int ii = idx_i.elem (i);
02826                       lhs.elem (ii) = scalar;
02827                     }
02828                 }
02829               else
02830                 {
02831                   (*current_liboctave_error_handler)
02832       ("A(I) = X: X must be a scalar or a matrix with the same size as I");
02833 
02834                   retval = 0;
02835                 }
02836             }
02837           // idx_vector::freeze() printed an error message for us.
02838         }
02839     }
02840   else
02841     {
02842       (*current_liboctave_error_handler)
02843         ("invalid number of indices for matrix expression");
02844 
02845       retval = 0;
02846     }
02847 
02848   lhs.clear_index ();
02849 
02850   return retval;
02851 }
02852 
02853 template <class LT, class RT>
02854 int
02855 assignN (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
02856 {
02857   int retval = 1;
02858 
02859   dim_vector rhs_dims = rhs.dims ();
02860 
02861   int rhs_dims_len = rhs_dims.length ();
02862 
02863   bool rhs_is_scalar = is_scalar (rhs_dims);
02864 
02865   int n_idx = lhs.index_count ();
02866 
02867   idx_vector *idx_vex = lhs.get_idx ();
02868 
02869   Array<idx_vector> idx = conv_to_array (idx_vex, n_idx);
02870 
02871   if (rhs_dims_len == 2 && rhs_dims(0) == 0 && rhs_dims(1) == 0)
02872     {
02873       lhs.maybe_delete_elements (idx, rfv);
02874     }
02875   else if (n_idx == 1)
02876     {
02877       idx_vector iidx = idx(0);
02878 
02879       if (liboctave_wfi_flag
02880           && ! (iidx.is_colon ()
02881                 || (iidx.one_zero_only ()
02882                     && iidx.orig_dimensions () == lhs.dims ())))
02883         (*current_liboctave_warning_handler)
02884           ("single index used for N-d array");
02885 
02886       int lhs_len = lhs.length ();
02887 
02888       int len = iidx.freeze (lhs_len, "N-d arrray");
02889 
02890       if (iidx)
02891         {
02892           if (len == 0)
02893             {
02894               if (! (rhs_dims.all_ones () || rhs_dims.any_zero ()))
02895                 {
02896                   (*current_liboctave_error_handler)
02897                     ("A([]) = X: X must be an empty matrix or scalar");
02898 
02899                   retval = 0;
02900                 }
02901             }
02902           else if (len == rhs.length ())
02903             {
02904               for (int i = 0; i < len; i++)
02905                 {
02906                   int ii = iidx.elem (i);
02907 
02908                   lhs.elem (ii) = rhs.elem (i);
02909                 }
02910             }
02911           else if (rhs_is_scalar)
02912             {
02913               RT scalar = rhs.elem (0);
02914 
02915               for (int i = 0; i < len; i++)
02916                 {
02917                   int ii = iidx.elem (i);
02918 
02919                   lhs.elem (ii) = scalar;
02920                 }
02921             }
02922           else
02923             {
02924               (*current_liboctave_error_handler)
02925                 ("A(I) = X: X must be a scalar or a matrix with the same size as I");
02926 
02927               retval = 0;
02928             }
02929 
02930           // idx_vector::freeze() printed an error message for us.
02931         }
02932     }
02933   else
02934     {
02935       // Maybe expand to more dimensions.
02936 
02937       dim_vector lhs_dims = lhs.dims ();
02938 
02939       int lhs_dims_len = lhs_dims.length ();
02940 
02941       dim_vector final_lhs_dims = lhs_dims;
02942 
02943       dim_vector frozen_len;
02944 
02945       int orig_lhs_dims_len = lhs_dims_len;
02946 
02947       bool orig_empty = lhs_dims.all_zero ();
02948 
02949       if (n_idx < lhs_dims_len)
02950         {
02951           // Collapse dimensions beyond last index.  Note that we
02952           // delay resizing LHS until we know that the assignment will
02953           // succeed.
02954 
02955           if (liboctave_wfi_flag && ! (idx(n_idx-1).is_colon ()))
02956             (*current_liboctave_warning_handler)
02957               ("fewer indices than dimensions for N-d array");
02958 
02959           for (int i = n_idx; i < lhs_dims_len; i++)
02960             lhs_dims(n_idx-1) *= lhs_dims(i);
02961 
02962           lhs_dims.resize (n_idx);
02963 
02964           lhs_dims_len = lhs_dims.length ();
02965         }
02966 
02967       // Resize.
02968 
02969       dim_vector new_dims;
02970       new_dims.resize (n_idx);
02971 
02972       for (int i = 0; i < n_idx; i++)
02973         {
02974           if (orig_empty)
02975             {
02976               // If index is a colon, resizing to RHS dimensions is
02977               // allowed because we started out empty.
02978 
02979               new_dims(i)
02980                 = (i < rhs_dims.length () && idx(i).is_colon ())
02981                 ? rhs_dims(i) : idx(i).max () + 1;
02982             }
02983           else
02984             {
02985               // We didn't start out with all zero dimensions, so if
02986               // index is a colon, it refers to the current LHS
02987               // dimension.  Otherwise, it is OK to enlarge to a
02988               // dimension given by the largest index, but if that 
02989               // index is a colon the new dimension is singleton.
02990 
02991               if (i < lhs_dims_len
02992                   && (idx(i).is_colon () || idx(i).max () < lhs_dims(i)))
02993                 new_dims(i) = lhs_dims(i);
02994               else if (! idx(i).is_colon ())
02995                 new_dims(i) = idx(i).max () + 1;
02996               else
02997                 new_dims(i) = 1;
02998             }
02999         }
03000 
03001       if (retval != 0)
03002         {
03003           if (! orig_empty
03004               && n_idx < orig_lhs_dims_len
03005               && new_dims(n_idx-1) != lhs_dims(n_idx-1))
03006             {
03007               // We reshaped and the last dimension changed.  This has to
03008               // be an error, because we don't know how to undo that
03009               // later...
03010 
03011               (*current_liboctave_error_handler)
03012                 ("array index %d (= %d) for assignment requires invalid resizing operation",
03013                  n_idx, new_dims(n_idx-1));
03014 
03015               retval = 0;
03016             }
03017           else
03018             {
03019               // Determine final dimensions for LHS and reset the
03020               // current size of the LHS.  Note that we delay actually
03021               // resizing LHS until we know that the assignment will
03022               // succeed.
03023 
03024               if (n_idx < orig_lhs_dims_len)
03025                 {
03026                   for (int i = 0; i < n_idx-1; i++)
03027                     final_lhs_dims(i) = new_dims(i);
03028                 }
03029               else
03030                 final_lhs_dims = new_dims;
03031 
03032               lhs_dims = new_dims;
03033 
03034               lhs_dims_len = lhs_dims.length ();
03035 
03036               frozen_len = freeze (idx, lhs_dims, true);
03037 
03038               if (rhs_is_scalar)
03039                 {
03040                   lhs.resize_and_fill (new_dims, rfv);
03041 
03042                   if  (! final_lhs_dims.any_zero ())
03043                     {
03044                       int n = Array<LT>::get_size (frozen_len);
03045 
03046                       Array<int> result_idx (lhs_dims_len, 0);
03047 
03048                       RT scalar = rhs.elem (0);
03049 
03050                       for (int i = 0; i < n; i++)
03051                         {
03052                           Array<int> elt_idx = get_elt_idx (idx, result_idx);
03053 
03054                           lhs.elem (elt_idx) = scalar;
03055 
03056                           increment_index (result_idx, frozen_len);
03057                         }
03058                     }
03059                 }
03060               else
03061                 {
03062                   // RHS is matrix or higher dimension.
03063 
03064                   // Check that non-singleton RHS dimensions conform to
03065                   // non-singleton LHS index dimensions.
03066 
03067                   dim_vector t_rhs_dims = rhs_dims.squeeze ();
03068                   dim_vector t_frozen_len = frozen_len.squeeze ();
03069 
03070                   // If after sqeezing out singleton dimensions, RHS is
03071                   // vector and LHS is vector, force them to have the same
03072                   // orientation so that operations like
03073                   //
03074                   //   a = zeros (3, 3, 3);
03075                   //   a(1:3,1,1) = [1,2,3];
03076                   //
03077                   // will work.
03078 
03079                   if (t_rhs_dims.length () == 2 && t_frozen_len.length () == 2
03080                       && ((t_rhs_dims.elem(1) == 1
03081                            && t_frozen_len.elem(0) == 1)
03082                           || (t_rhs_dims.elem(0) == 1
03083                               && t_frozen_len.elem(1) == 1)))
03084                     {
03085                       int t0 = t_rhs_dims.elem(0);
03086                       t_rhs_dims.elem(0) = t_rhs_dims.elem(1);
03087                       t_rhs_dims.elem(1) = t0;
03088                     }
03089 
03090                   if (t_rhs_dims != t_frozen_len)
03091                     {
03092                       (*current_liboctave_error_handler)
03093                         ("A(IDX-LIST) = X: X must be a scalar or size of X must equal number of elements indexed by IDX-LIST");
03094 
03095                           retval = 0;
03096                     }
03097                   else
03098                     {
03099                       lhs.resize_and_fill (new_dims, rfv);
03100 
03101                       if  (! final_lhs_dims.any_zero ())
03102                         {
03103                           int n = Array<LT>::get_size (frozen_len);
03104 
03105                           Array<int> result_idx (lhs_dims_len, 0);
03106 
03107                           for (int i = 0; i < n; i++)
03108                             {
03109                               Array<int> elt_idx = get_elt_idx (idx, result_idx);
03110 
03111                               lhs.elem (elt_idx) = rhs.elem (i);
03112 
03113                               increment_index (result_idx, frozen_len);
03114                             }
03115                         }
03116                     }
03117                 }
03118             }
03119         }
03120 
03121       if (retval != 0)
03122         lhs.resize (final_lhs_dims);
03123     }
03124 
03125   if (retval != 0)
03126     lhs.chop_trailing_singletons ();
03127 
03128   lhs.clear_index ();
03129 
03130   return retval;
03131 }
03132 
03133 template <class T>
03134 void
03135 Array<T>::print_info (std::ostream& os, const std::string& prefix) const
03136 {
03137   os << prefix << "rep address: " << rep << "\n"
03138      << prefix << "rep->len:    " << rep->len << "\n"
03139      << prefix << "rep->data:   " << static_cast<void *> (rep->data) << "\n"
03140      << prefix << "rep->count:  " << rep->count << "\n";
03141 
03142   // 2D info:
03143   //
03144   //     << pefix << "rows: " << rows () << "\n"
03145   //     << prefix << "cols: " << cols () << "\n";
03146 }
03147 
03148 /*
03149 ;;; Local Variables: ***
03150 ;;; mode: C++ ***
03151 ;;; End: ***
03152 */

Wed Dec 29 11:50:53 2004に生成されました。 doxygen1.2.18
SEO [PR] 爆速!無料ブログ 無料ホームページ開設 無料ライブ放送