← Back to branch summary

~gabriel1984sibiu/octave/octave

« back to all changes in this revision

Viewing changes to liboctave/array/CRowVector.cc

  • Committer: Grevutiu Gabriel
  • Date: 2014-01-02 13:05:54 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20140102130554-3r7ivdjln1ni6kcg
New version (3.8.0) from upstream.

Show diffs side-by-side

added added

removed removed

Lines of Context:
liboctave/array/CRowVector.cc
 
1
// RowVector manipulations.
 
2
/*
 
3
 
 
4
Copyright (C) 1994-2013 John W. Eaton
 
5
 
 
6
This file is part of Octave.
 
7
 
 
8
Octave is free software; you can redistribute it and/or modify it
 
9
under the terms of the GNU General Public License as published by the
 
10
Free Software Foundation; either version 3 of the License, or (at your
 
11
option) any later version.
 
12
 
 
13
Octave is distributed in the hope that it will be useful, but WITHOUT
 
14
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 
15
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 
16
for more details.
 
17
 
 
18
You should have received a copy of the GNU General Public License
 
19
along with Octave; see the file COPYING.  If not, see
 
20
<http://www.gnu.org/licenses/>.
 
21
 
 
22
*/
 
23
 
 
24
#ifdef HAVE_CONFIG_H
 
25
#include <config.h>
 
26
#endif
 
27
 
 
28
#include <iostream>
 
29
 
 
30
#include "Array-util.h"
 
31
#include "f77-fcn.h"
 
32
#include "functor.h"
 
33
#include "lo-error.h"
 
34
#include "mx-base.h"
 
35
#include "mx-inlines.cc"
 
36
#include "oct-cmplx.h"
 
37
 
 
38
// Fortran functions we call.
 
39
 
 
40
extern "C"
 
41
{
 
42
  F77_RET_T
 
43
  F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
 
44
                           const octave_idx_type&, const octave_idx_type&,
 
45
                           const Complex&, const Complex*,
 
46
                           const octave_idx_type&, const Complex*,
 
47
                           const octave_idx_type&, const Complex&, Complex*,
 
48
                           const octave_idx_type&
 
49
                           F77_CHAR_ARG_LEN_DECL);
 
50
 
 
51
  F77_RET_T
 
52
  F77_FUNC (xzdotu, XZDOTU) (const octave_idx_type&, const Complex*,
 
53
                             const octave_idx_type&, const Complex*,
 
54
                             const octave_idx_type&, Complex&);
 
55
}
 
56
 
 
57
// Complex Row Vector class
 
58
 
 
59
bool
 
60
ComplexRowVector::operator == (const ComplexRowVector& a) const
 
61
{
 
62
  octave_idx_type len = length ();
 
63
  if (len != a.length ())
 
64
    return 0;
 
65
  return mx_inline_equal (len, data (), a.data ());
 
66
}
 
67
 
 
68
bool
 
69
ComplexRowVector::operator != (const ComplexRowVector& a) const
 
70
{
 
71
  return !(*this == a);
 
72
}
 
73
 
 
74
// destructive insert/delete/reorder operations
 
75
 
 
76
ComplexRowVector&
 
77
ComplexRowVector::insert (const RowVector& a, octave_idx_type c)
 
78
{
 
79
  octave_idx_type a_len = a.length ();
 
80
 
 
81
  if (c < 0 || c + a_len > length ())
 
82
    {
 
83
      (*current_liboctave_error_handler) ("range error for insert");
 
84
      return *this;
 
85
    }
 
86
 
 
87
  if (a_len > 0)
 
88
    {
 
89
      make_unique ();
 
90
 
 
91
      for (octave_idx_type i = 0; i < a_len; i++)
 
92
        xelem (c+i) = a.elem (i);
 
93
    }
 
94
 
 
95
  return *this;
 
96
}
 
97
 
 
98
ComplexRowVector&
 
99
ComplexRowVector::insert (const ComplexRowVector& a, octave_idx_type c)
 
100
{
 
101
  octave_idx_type a_len = a.length ();
 
102
 
 
103
  if (c < 0 || c + a_len > length ())
 
104
    {
 
105
      (*current_liboctave_error_handler) ("range error for insert");
 
106
      return *this;
 
107
    }
 
108
 
 
109
  if (a_len > 0)
 
110
    {
 
111
      make_unique ();
 
112
 
 
113
      for (octave_idx_type i = 0; i < a_len; i++)
 
114
        xelem (c+i) = a.elem (i);
 
115
    }
 
116
 
 
117
  return *this;
 
118
}
 
119
 
 
120
ComplexRowVector&
 
121
ComplexRowVector::fill (double val)
 
122
{
 
123
  octave_idx_type len = length ();
 
124
 
 
125
  if (len > 0)
 
126
    {
 
127
      make_unique ();
 
128
 
 
129
      for (octave_idx_type i = 0; i < len; i++)
 
130
        xelem (i) = val;
 
131
    }
 
132
 
 
133
  return *this;
 
134
}
 
135
 
 
136
ComplexRowVector&
 
137
ComplexRowVector::fill (const Complex& val)
 
138
{
 
139
  octave_idx_type len = length ();
 
140
 
 
141
  if (len > 0)
 
142
    {
 
143
      make_unique ();
 
144
 
 
145
      for (octave_idx_type i = 0; i < len; i++)
 
146
        xelem (i) = val;
 
147
    }
 
148
 
 
149
  return *this;
 
150
}
 
151
 
 
152
ComplexRowVector&
 
153
ComplexRowVector::fill (double val, octave_idx_type c1, octave_idx_type c2)
 
154
{
 
155
  octave_idx_type len = length ();
 
156
 
 
157
  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
 
158
    {
 
159
      (*current_liboctave_error_handler) ("range error for fill");
 
160
      return *this;
 
161
    }
 
162
 
 
163
  if (c1 > c2) { std::swap (c1, c2); }
 
164
 
 
165
  if (c2 >= c1)
 
166
    {
 
167
      make_unique ();
 
168
 
 
169
      for (octave_idx_type i = c1; i <= c2; i++)
 
170
        xelem (i) = val;
 
171
    }
 
172
 
 
173
  return *this;
 
174
}
 
175
 
 
176
ComplexRowVector&
 
177
ComplexRowVector::fill (const Complex& val,
 
178
                        octave_idx_type c1, octave_idx_type c2)
 
179
{
 
180
  octave_idx_type len = length ();
 
181
 
 
182
  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
 
183
    {
 
184
      (*current_liboctave_error_handler) ("range error for fill");
 
185
      return *this;
 
186
    }
 
187
 
 
188
  if (c1 > c2) { std::swap (c1, c2); }
 
189
 
 
190
  if (c2 >= c1)
 
191
    {
 
192
      make_unique ();
 
193
 
 
194
      for (octave_idx_type i = c1; i <= c2; i++)
 
195
        xelem (i) = val;
 
196
    }
 
197
 
 
198
  return *this;
 
199
}
 
200
 
 
201
ComplexRowVector
 
202
ComplexRowVector::append (const RowVector& a) const
 
203
{
 
204
  octave_idx_type len = length ();
 
205
  octave_idx_type nc_insert = len;
 
206
  ComplexRowVector retval (len + a.length ());
 
207
  retval.insert (*this, 0);
 
208
  retval.insert (a, nc_insert);
 
209
  return retval;
 
210
}
 
211
 
 
212
ComplexRowVector
 
213
ComplexRowVector::append (const ComplexRowVector& a) const
 
214
{
 
215
  octave_idx_type len = length ();
 
216
  octave_idx_type nc_insert = len;
 
217
  ComplexRowVector retval (len + a.length ());
 
218
  retval.insert (*this, 0);
 
219
  retval.insert (a, nc_insert);
 
220
  return retval;
 
221
}
 
222
 
 
223
ComplexColumnVector
 
224
ComplexRowVector::hermitian (void) const
 
225
{
 
226
  return MArray<Complex>::hermitian (std::conj);
 
227
}
 
228
 
 
229
ComplexColumnVector
 
230
ComplexRowVector::transpose (void) const
 
231
{
 
232
  return MArray<Complex>::transpose ();
 
233
}
 
234
 
 
235
ComplexRowVector
 
236
conj (const ComplexRowVector& a)
 
237
{
 
238
  return do_mx_unary_map<Complex, Complex, std::conj<double> > (a);
 
239
}
 
240
 
 
241
// resize is the destructive equivalent for this one
 
242
 
 
243
ComplexRowVector
 
244
ComplexRowVector::extract (octave_idx_type c1, octave_idx_type c2) const
 
245
{
 
246
  if (c1 > c2) { std::swap (c1, c2); }
 
247
 
 
248
  octave_idx_type new_c = c2 - c1 + 1;
 
249
 
 
250
  ComplexRowVector result (new_c);
 
251
 
 
252
  for (octave_idx_type i = 0; i < new_c; i++)
 
253
    result.elem (i) = elem (c1+i);
 
254
 
 
255
  return result;
 
256
}
 
257
 
 
258
ComplexRowVector
 
259
ComplexRowVector::extract_n (octave_idx_type r1, octave_idx_type n) const
 
260
{
 
261
  ComplexRowVector result (n);
 
262
 
 
263
  for (octave_idx_type i = 0; i < n; i++)
 
264
    result.elem (i) = elem (r1+i);
 
265
 
 
266
  return result;
 
267
}
 
268
 
 
269
// row vector by row vector -> row vector operations
 
270
 
 
271
ComplexRowVector&
 
272
ComplexRowVector::operator += (const RowVector& a)
 
273
{
 
274
  octave_idx_type len = length ();
 
275
 
 
276
  octave_idx_type a_len = a.length ();
 
277
 
 
278
  if (len != a_len)
 
279
    {
 
280
      gripe_nonconformant ("operator +=", len, a_len);
 
281
      return *this;
 
282
    }
 
283
 
 
284
  if (len == 0)
 
285
    return *this;
 
286
 
 
287
  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
 
288
 
 
289
  mx_inline_add2 (len, d, a.data ());
 
290
  return *this;
 
291
}
 
292
 
 
293
ComplexRowVector&
 
294
ComplexRowVector::operator -= (const RowVector& a)
 
295
{
 
296
  octave_idx_type len = length ();
 
297
 
 
298
  octave_idx_type a_len = a.length ();
 
299
 
 
300
  if (len != a_len)
 
301
    {
 
302
      gripe_nonconformant ("operator -=", len, a_len);
 
303
      return *this;
 
304
    }
 
305
 
 
306
  if (len == 0)
 
307
    return *this;
 
308
 
 
309
  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
 
310
 
 
311
  mx_inline_sub2 (len, d, a.data ());
 
312
  return *this;
 
313
}
 
314
 
 
315
// row vector by matrix -> row vector
 
316
 
 
317
ComplexRowVector
 
318
operator * (const ComplexRowVector& v, const ComplexMatrix& a)
 
319
{
 
320
  ComplexRowVector retval;
 
321
 
 
322
  octave_idx_type len = v.length ();
 
323
 
 
324
  octave_idx_type a_nr = a.rows ();
 
325
  octave_idx_type a_nc = a.cols ();
 
326
 
 
327
  if (a_nr != len)
 
328
    gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
 
329
  else
 
330
    {
 
331
      if (len == 0)
 
332
        retval.resize (a_nc, 0.0);
 
333
      else
 
334
        {
 
335
          // Transpose A to form A'*x == (x'*A)'
 
336
 
 
337
          octave_idx_type ld = a_nr;
 
338
 
 
339
          retval.resize (a_nc);
 
340
          Complex *y = retval.fortran_vec ();
 
341
 
 
342
          F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
 
343
                                   a_nr, a_nc, 1.0, a.data (),
 
344
                                   ld, v.data (), 1, 0.0, y, 1
 
345
                                   F77_CHAR_ARG_LEN (1)));
 
346
        }
 
347
    }
 
348
 
 
349
  return retval;
 
350
}
 
351
 
 
352
ComplexRowVector
 
353
operator * (const RowVector& v, const ComplexMatrix& a)
 
354
{
 
355
  ComplexRowVector tmp (v);
 
356
  return tmp * a;
 
357
}
 
358
 
 
359
// other operations
 
360
 
 
361
Complex
 
362
ComplexRowVector::min (void) const
 
363
{
 
364
  octave_idx_type len = length ();
 
365
  if (len == 0)
 
366
    return Complex (0.0);
 
367
 
 
368
  Complex res = elem (0);
 
369
  double absres = std::abs (res);
 
370
 
 
371
  for (octave_idx_type i = 1; i < len; i++)
 
372
    if (std::abs (elem (i)) < absres)
 
373
      {
 
374
        res = elem (i);
 
375
        absres = std::abs (res);
 
376
      }
 
377
 
 
378
  return res;
 
379
}
 
380
 
 
381
Complex
 
382
ComplexRowVector::max (void) const
 
383
{
 
384
  octave_idx_type len = length ();
 
385
  if (len == 0)
 
386
    return Complex (0.0);
 
387
 
 
388
  Complex res = elem (0);
 
389
  double absres = std::abs (res);
 
390
 
 
391
  for (octave_idx_type i = 1; i < len; i++)
 
392
    if (std::abs (elem (i)) > absres)
 
393
      {
 
394
        res = elem (i);
 
395
        absres = std::abs (res);
 
396
      }
 
397
 
 
398
  return res;
 
399
}
 
400
 
 
401
// i/o
 
402
 
 
403
std::ostream&
 
404
operator << (std::ostream& os, const ComplexRowVector& a)
 
405
{
 
406
//  int field_width = os.precision () + 7;
 
407
  for (octave_idx_type i = 0; i < a.length (); i++)
 
408
    os << " " /* setw (field_width) */ << a.elem (i);
 
409
  return os;
 
410
}
 
411
 
 
412
std::istream&
 
413
operator >> (std::istream& is, ComplexRowVector& a)
 
414
{
 
415
  octave_idx_type len = a.length ();
 
416
 
 
417
  if (len > 0)
 
418
    {
 
419
      Complex tmp;
 
420
      for (octave_idx_type i = 0; i < len; i++)
 
421
        {
 
422
          is >> tmp;
 
423
          if (is)
 
424
            a.elem (i) = tmp;
 
425
          else
 
426
            break;
 
427
        }
 
428
    }
 
429
  return is;
 
430
}
 
431
 
 
432
// row vector by column vector -> scalar
 
433
 
 
434
// row vector by column vector -> scalar
 
435
 
 
436
Complex
 
437
operator * (const ComplexRowVector& v, const ColumnVector& a)
 
438
{
 
439
  ComplexColumnVector tmp (a);
 
440
  return v * tmp;
 
441
}
 
442
 
 
443
Complex
 
444
operator * (const ComplexRowVector& v, const ComplexColumnVector& a)
 
445
{
 
446
  Complex retval (0.0, 0.0);
 
447
 
 
448
  octave_idx_type len = v.length ();
 
449
 
 
450
  octave_idx_type a_len = a.length ();
 
451
 
 
452
  if (len != a_len)
 
453
    gripe_nonconformant ("operator *", len, a_len);
 
454
  else if (len != 0)
 
455
    F77_FUNC (xzdotu, XZDOTU) (len, v.data (), 1, a.data (), 1, retval);
 
456
 
 
457
  return retval;
 
458
}
 
459
 
 
460
// other operations
 
461
 
 
462
ComplexRowVector
 
463
linspace (const Complex& x1, const Complex& x2, octave_idx_type n)
 
464
{
 
465
  if (n < 1) n = 1;
 
466
 
 
467
  NoAlias<ComplexRowVector> retval (n);
 
468
 
 
469
  Complex delta = (x2 - x1) / (n - 1.0);
 
470
  retval(0) = x1;
 
471
  for (octave_idx_type i = 1; i < n-1; i++)
 
472
    retval(i) = x1 + static_cast<double> (i)*delta;
 
473
  retval(n-1) = x2;
 
474
 
 
475
  return retval;
 
476
}