Mercurial > hg > octave-lyh

comparison liboctave/CRowVector.cc @ 5275:23b37da9fd5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents deed800e7bef
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
38 38
39 extern "C" 39 extern "C"
40 { 40 {
41 F77_RET_T 41 F77_RET_T
42 F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL, 42 F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
43 const int&, const int&, const Complex&, 43 const octave_idx_type&, const octave_idx_type&, const Complex&,
44 const Complex*, const int&, const Complex*, 44 const Complex*, const octave_idx_type&, const Complex*,
45 const int&, const Complex&, Complex*, const int& 45 const octave_idx_type&, const Complex&, Complex*, const octave_idx_type&
46 F77_CHAR_ARG_LEN_DECL); 46 F77_CHAR_ARG_LEN_DECL);
47 } 47 }
48 48
49 // Complex Row Vector class 49 // Complex Row Vector class
50 50
51 ComplexRowVector::ComplexRowVector (const RowVector& a) 51 ComplexRowVector::ComplexRowVector (const RowVector& a)
52 : MArray<Complex> (a.length ()) 52 : MArray<Complex> (a.length ())
53 { 53 {
54 for (int i = 0; i < length (); i++) 54 for (octave_idx_type i = 0; i < length (); i++)
55 elem (i) = a.elem (i); 55 elem (i) = a.elem (i);
56 } 56 }
57 57
58 bool 58 bool
59 ComplexRowVector::operator == (const ComplexRowVector& a) const 59 ComplexRowVector::operator == (const ComplexRowVector& a) const
60 { 60 {
61 int len = length (); 61 octave_idx_type len = length ();
62 if (len != a.length ()) 62 if (len != a.length ())
63 return 0; 63 return 0;
64 return mx_inline_equal (data (), a.data (), len); 64 return mx_inline_equal (data (), a.data (), len);
65 } 65 }
66 66
71 } 71 }
72 72
73 // destructive insert/delete/reorder operations 73 // destructive insert/delete/reorder operations
74 74
75 ComplexRowVector& 75 ComplexRowVector&
76 ComplexRowVector::insert (const RowVector& a, int c) 76 ComplexRowVector::insert (const RowVector& a, octave_idx_type c)
77 { 77 {
78 int a_len = a.length (); 78 octave_idx_type a_len = a.length ();
79 79
80 if (c < 0 || c + a_len > length ()) 80 if (c < 0 || c + a_len > length ())
81 { 81 {
82 (*current_liboctave_error_handler) ("range error for insert"); 82 (*current_liboctave_error_handler) ("range error for insert");
83 return *this; 83 return *this;
85 85
86 if (a_len > 0) 86 if (a_len > 0)
87 { 87 {
88 make_unique (); 88 make_unique ();
89 89
90 for (int i = 0; i < a_len; i++) 90 for (octave_idx_type i = 0; i < a_len; i++)
91 xelem (c+i) = a.elem (i); 91 xelem (c+i) = a.elem (i);
92 } 92 }
93 93
94 return *this; 94 return *this;
95 } 95 }
96 96
97 ComplexRowVector& 97 ComplexRowVector&
98 ComplexRowVector::insert (const ComplexRowVector& a, int c) 98 ComplexRowVector::insert (const ComplexRowVector& a, octave_idx_type c)
99 { 99 {
100 int a_len = a.length (); 100 octave_idx_type a_len = a.length ();
101 101
102 if (c < 0 || c + a_len > length ()) 102 if (c < 0 || c + a_len > length ())
103 { 103 {
104 (*current_liboctave_error_handler) ("range error for insert"); 104 (*current_liboctave_error_handler) ("range error for insert");
105 return *this; 105 return *this;
107 107
108 if (a_len > 0) 108 if (a_len > 0)
109 { 109 {
110 make_unique (); 110 make_unique ();
111 111
112 for (int i = 0; i < a_len; i++) 112 for (octave_idx_type i = 0; i < a_len; i++)
113 xelem (c+i) = a.elem (i); 113 xelem (c+i) = a.elem (i);
114 } 114 }
115 115
116 return *this; 116 return *this;
117 } 117 }
118 118
119 ComplexRowVector& 119 ComplexRowVector&
120 ComplexRowVector::fill (double val) 120 ComplexRowVector::fill (double val)
121 { 121 {
122 int len = length (); 122 octave_idx_type len = length ();
123 123
124 if (len > 0) 124 if (len > 0)
125 { 125 {
126 make_unique (); 126 make_unique ();
127 127
128 for (int i = 0; i < len; i++) 128 for (octave_idx_type i = 0; i < len; i++)
129 xelem (i) = val; 129 xelem (i) = val;
130 } 130 }
131 131
132 return *this; 132 return *this;
133 } 133 }
134 134
135 ComplexRowVector& 135 ComplexRowVector&
136 ComplexRowVector::fill (const Complex& val) 136 ComplexRowVector::fill (const Complex& val)
137 { 137 {
138 int len = length (); 138 octave_idx_type len = length ();
139 139
140 if (len > 0) 140 if (len > 0)
141 { 141 {
142 make_unique (); 142 make_unique ();
143 143
144 for (int i = 0; i < len; i++) 144 for (octave_idx_type i = 0; i < len; i++)
145 xelem (i) = val; 145 xelem (i) = val;
146 } 146 }
147 147
148 return *this; 148 return *this;
149 } 149 }
150 150
151 ComplexRowVector& 151 ComplexRowVector&
152 ComplexRowVector::fill (double val, int c1, int c2) 152 ComplexRowVector::fill (double val, octave_idx_type c1, octave_idx_type c2)
153 { 153 {
154 int len = length (); 154 octave_idx_type len = length ();
155 155
156 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len) 156 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
157 { 157 {
158 (*current_liboctave_error_handler) ("range error for fill"); 158 (*current_liboctave_error_handler) ("range error for fill");
159 return *this; 159 return *this;
160 } 160 }
161 161
162 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } 162 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
163 163
164 if (c2 >= c1) 164 if (c2 >= c1)
165 { 165 {
166 make_unique (); 166 make_unique ();
167 167
168 for (int i = c1; i <= c2; i++) 168 for (octave_idx_type i = c1; i <= c2; i++)
169 xelem (i) = val; 169 xelem (i) = val;
170 } 170 }
171 171
172 return *this; 172 return *this;
173 } 173 }
174 174
175 ComplexRowVector& 175 ComplexRowVector&
176 ComplexRowVector::fill (const Complex& val, int c1, int c2) 176 ComplexRowVector::fill (const Complex& val, octave_idx_type c1, octave_idx_type c2)
177 { 177 {
178 int len = length (); 178 octave_idx_type len = length ();
179 179
180 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len) 180 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
181 { 181 {
182 (*current_liboctave_error_handler) ("range error for fill"); 182 (*current_liboctave_error_handler) ("range error for fill");
183 return *this; 183 return *this;
184 } 184 }
185 185
186 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } 186 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
187 187
188 if (c2 >= c1) 188 if (c2 >= c1)
189 { 189 {
190 make_unique (); 190 make_unique ();
191 191
192 for (int i = c1; i <= c2; i++) 192 for (octave_idx_type i = c1; i <= c2; i++)
193 xelem (i) = val; 193 xelem (i) = val;
194 } 194 }
195 195
196 return *this; 196 return *this;
197 } 197 }
198 198
199 ComplexRowVector 199 ComplexRowVector
200 ComplexRowVector::append (const RowVector& a) const 200 ComplexRowVector::append (const RowVector& a) const
201 { 201 {
202 int len = length (); 202 octave_idx_type len = length ();
203 int nc_insert = len; 203 octave_idx_type nc_insert = len;
204 ComplexRowVector retval (len + a.length ()); 204 ComplexRowVector retval (len + a.length ());
205 retval.insert (*this, 0); 205 retval.insert (*this, 0);
206 retval.insert (a, nc_insert); 206 retval.insert (a, nc_insert);
207 return retval; 207 return retval;
208 } 208 }
209 209
210 ComplexRowVector 210 ComplexRowVector
211 ComplexRowVector::append (const ComplexRowVector& a) const 211 ComplexRowVector::append (const ComplexRowVector& a) const
212 { 212 {
213 int len = length (); 213 octave_idx_type len = length ();
214 int nc_insert = len; 214 octave_idx_type nc_insert = len;
215 ComplexRowVector retval (len + a.length ()); 215 ComplexRowVector retval (len + a.length ());
216 retval.insert (*this, 0); 216 retval.insert (*this, 0);
217 retval.insert (a, nc_insert); 217 retval.insert (a, nc_insert);
218 return retval; 218 return retval;
219 } 219 }
220 220
221 ComplexColumnVector 221 ComplexColumnVector
222 ComplexRowVector::hermitian (void) const 222 ComplexRowVector::hermitian (void) const
223 { 223 {
224 int len = length (); 224 octave_idx_type len = length ();
225 return ComplexColumnVector (mx_inline_conj_dup (data (), len), len); 225 return ComplexColumnVector (mx_inline_conj_dup (data (), len), len);
226 } 226 }
227 227
228 ComplexColumnVector 228 ComplexColumnVector
229 ComplexRowVector::transpose (void) const 229 ComplexRowVector::transpose (void) const
232 } 232 }
233 233
234 ComplexRowVector 234 ComplexRowVector
235 conj (const ComplexRowVector& a) 235 conj (const ComplexRowVector& a)
236 { 236 {
237 int a_len = a.length (); 237 octave_idx_type a_len = a.length ();
238 ComplexRowVector retval; 238 ComplexRowVector retval;
239 if (a_len > 0) 239 if (a_len > 0)
240 retval = ComplexRowVector (mx_inline_conj_dup (a.data (), a_len), a_len); 240 retval = ComplexRowVector (mx_inline_conj_dup (a.data (), a_len), a_len);
241 return retval; 241 return retval;
242 } 242 }
243 243
244 // resize is the destructive equivalent for this one 244 // resize is the destructive equivalent for this one
245 245
246 ComplexRowVector 246 ComplexRowVector
247 ComplexRowVector::extract (int c1, int c2) const 247 ComplexRowVector::extract (octave_idx_type c1, octave_idx_type c2) const
248 { 248 {
249 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } 249 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
250 250
251 int new_c = c2 - c1 + 1; 251 octave_idx_type new_c = c2 - c1 + 1;
252 252
253 ComplexRowVector result (new_c); 253 ComplexRowVector result (new_c);
254 254
255 for (int i = 0; i < new_c; i++) 255 for (octave_idx_type i = 0; i < new_c; i++)
256 result.elem (i) = elem (c1+i); 256 result.elem (i) = elem (c1+i);
257 257
258 return result; 258 return result;
259 } 259 }
260 260
261 ComplexRowVector 261 ComplexRowVector
262 ComplexRowVector::extract_n (int r1, int n) const 262 ComplexRowVector::extract_n (octave_idx_type r1, octave_idx_type n) const
263 { 263 {
264 ComplexRowVector result (n); 264 ComplexRowVector result (n);
265 265
266 for (int i = 0; i < n; i++) 266 for (octave_idx_type i = 0; i < n; i++)
267 result.elem (i) = elem (r1+i); 267 result.elem (i) = elem (r1+i);
268 268
269 return result; 269 return result;
270 } 270 }
271 271
272 // row vector by row vector -> row vector operations 272 // row vector by row vector -> row vector operations
273 273
274 ComplexRowVector& 274 ComplexRowVector&
275 ComplexRowVector::operator += (const RowVector& a) 275 ComplexRowVector::operator += (const RowVector& a)
276 { 276 {
277 int len = length (); 277 octave_idx_type len = length ();
278 278
279 int a_len = a.length (); 279 octave_idx_type a_len = a.length ();
280 280
281 if (len != a_len) 281 if (len != a_len)
282 { 282 {
283 gripe_nonconformant ("operator +=", len, a_len); 283 gripe_nonconformant ("operator +=", len, a_len);
284 return *this; 284 return *this;
294 } 294 }
295 295
296 ComplexRowVector& 296 ComplexRowVector&
297 ComplexRowVector::operator -= (const RowVector& a) 297 ComplexRowVector::operator -= (const RowVector& a)
298 { 298 {
299 int len = length (); 299 octave_idx_type len = length ();
300 300
301 int a_len = a.length (); 301 octave_idx_type a_len = a.length ();
302 302
303 if (len != a_len) 303 if (len != a_len)
304 { 304 {
305 gripe_nonconformant ("operator -=", len, a_len); 305 gripe_nonconformant ("operator -=", len, a_len);
306 return *this; 306 return *this;
320 ComplexRowVector 320 ComplexRowVector
321 operator * (const ComplexRowVector& v, const ComplexMatrix& a) 321 operator * (const ComplexRowVector& v, const ComplexMatrix& a)
322 { 322 {
323 ComplexRowVector retval; 323 ComplexRowVector retval;
324 324
325 int len = v.length (); 325 octave_idx_type len = v.length ();
326 326
327 int a_nr = a.rows (); 327 octave_idx_type a_nr = a.rows ();
328 int a_nc = a.cols (); 328 octave_idx_type a_nc = a.cols ();
329 329
330 if (a_nr != len) 330 if (a_nr != len)
331 gripe_nonconformant ("operator *", 1, len, a_nr, a_nc); 331 gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
332 else 332 else
333 { 333 {
335 retval.resize (a_nc, 0.0); 335 retval.resize (a_nc, 0.0);
336 else 336 else
337 { 337 {
338 // Transpose A to form A'*x == (x'*A)' 338 // Transpose A to form A'*x == (x'*A)'
339 339
340 int ld = a_nr; 340 octave_idx_type ld = a_nr;
341 341
342 retval.resize (a_nc); 342 retval.resize (a_nc);
343 Complex *y = retval.fortran_vec (); 343 Complex *y = retval.fortran_vec ();
344 344
345 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1), 345 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
375 RowVector 375 RowVector
376 ComplexRowVector::map (d_c_Mapper f) const 376 ComplexRowVector::map (d_c_Mapper f) const
377 { 377 {
378 const Complex *d = data (); 378 const Complex *d = data ();
379 379
380 int len = length (); 380 octave_idx_type len = length ();
381 381
382 RowVector retval (len); 382 RowVector retval (len);
383 383
384 double *r = retval.fortran_vec (); 384 double *r = retval.fortran_vec ();
385 385
386 for (int i = 0; i < len; i++) 386 for (octave_idx_type i = 0; i < len; i++)
387 r[i] = f (d[i]); 387 r[i] = f (d[i]);
388 388
389 return retval; 389 return retval;
390 } 390 }
391 391
392 ComplexRowVector& 392 ComplexRowVector&
393 ComplexRowVector::apply (c_c_Mapper f) 393 ComplexRowVector::apply (c_c_Mapper f)
394 { 394 {
395 Complex *d = fortran_vec (); // Ensures only one reference to my privates! 395 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
396 396
397 for (int i = 0; i < length (); i++) 397 for (octave_idx_type i = 0; i < length (); i++)
398 d[i] = f (d[i]); 398 d[i] = f (d[i]);
399 399
400 return *this; 400 return *this;
401 } 401 }
402 402
403 Complex 403 Complex
404 ComplexRowVector::min (void) const 404 ComplexRowVector::min (void) const
405 { 405 {
406 int len = length (); 406 octave_idx_type len = length ();
407 if (len == 0) 407 if (len == 0)
408 return Complex (0.0); 408 return Complex (0.0);
409 409
410 Complex res = elem (0); 410 Complex res = elem (0);
411 double absres = std::abs (res); 411 double absres = std::abs (res);
412 412
413 for (int i = 1; i < len; i++) 413 for (octave_idx_type i = 1; i < len; i++)
414 if (std::abs (elem (i)) < absres) 414 if (std::abs (elem (i)) < absres)
415 { 415 {
416 res = elem (i); 416 res = elem (i);
417 absres = std::abs (res); 417 absres = std::abs (res);
418 } 418 }
421 } 421 }
422 422
423 Complex 423 Complex
424 ComplexRowVector::max (void) const 424 ComplexRowVector::max (void) const
425 { 425 {
426 int len = length (); 426 octave_idx_type len = length ();
427 if (len == 0) 427 if (len == 0)
428 return Complex (0.0); 428 return Complex (0.0);
429 429
430 Complex res = elem (0); 430 Complex res = elem (0);
431 double absres = std::abs (res); 431 double absres = std::abs (res);
432 432
433 for (int i = 1; i < len; i++) 433 for (octave_idx_type i = 1; i < len; i++)
434 if (std::abs (elem (i)) > absres) 434 if (std::abs (elem (i)) > absres)
435 { 435 {
436 res = elem (i); 436 res = elem (i);
437 absres = std::abs (res); 437 absres = std::abs (res);
438 } 438 }
444 444
445 std::ostream& 445 std::ostream&
446 operator << (std::ostream& os, const ComplexRowVector& a) 446 operator << (std::ostream& os, const ComplexRowVector& a)
447 { 447 {
448 // int field_width = os.precision () + 7; 448 // int field_width = os.precision () + 7;
449 for (int i = 0; i < a.length (); i++) 449 for (octave_idx_type i = 0; i < a.length (); i++)
450 os << " " /* setw (field_width) */ << a.elem (i); 450 os << " " /* setw (field_width) */ << a.elem (i);
451 return os; 451 return os;
452 } 452 }
453 453
454 std::istream& 454 std::istream&
455 operator >> (std::istream& is, ComplexRowVector& a) 455 operator >> (std::istream& is, ComplexRowVector& a)
456 { 456 {
457 int len = a.length(); 457 octave_idx_type len = a.length();
458 458
459 if (len < 1) 459 if (len < 1)
460 is.clear (std::ios::badbit); 460 is.clear (std::ios::badbit);
461 else 461 else
462 { 462 {
463 Complex tmp; 463 Complex tmp;
464 for (int i = 0; i < len; i++) 464 for (octave_idx_type i = 0; i < len; i++)
465 { 465 {
466 is >> tmp; 466 is >> tmp;
467 if (is) 467 if (is)
468 a.elem (i) = tmp; 468 a.elem (i) = tmp;
469 else 469 else
485 } 485 }
486 486
487 Complex 487 Complex
488 operator * (const ComplexRowVector& v, const ComplexColumnVector& a) 488 operator * (const ComplexRowVector& v, const ComplexColumnVector& a)
489 { 489 {
490 int len = v.length (); 490 octave_idx_type len = v.length ();
491 491
492 int a_len = a.length (); 492 octave_idx_type a_len = a.length ();
493 493
494 if (len != a_len) 494 if (len != a_len)
495 { 495 {
496 gripe_nonconformant ("operator *", len, a_len); 496 gripe_nonconformant ("operator *", len, a_len);
497 return 0.0; 497 return 0.0;
498 } 498 }
499 499
500 Complex retval (0.0, 0.0); 500 Complex retval (0.0, 0.0);
501 501
502 for (int i = 0; i < len; i++) 502 for (octave_idx_type i = 0; i < len; i++)
503 retval += v.elem (i) * a.elem (i); 503 retval += v.elem (i) * a.elem (i);
504 504
505 return retval; 505 return retval;
506 } 506 }
507 507
508 // other operations 508 // other operations
509 509
510 ComplexRowVector 510 ComplexRowVector
511 linspace (const Complex& x1, const Complex& x2, int n) 511 linspace (const Complex& x1, const Complex& x2, octave_idx_type n)
512 { 512 {
513 ComplexRowVector retval; 513 ComplexRowVector retval;
514 514
515 if (n > 0) 515 if (n > 0)
516 { 516 {
517 retval.resize (n); 517 retval.resize (n);
518 Complex delta = (x2 - x1) / (n - 1.0); 518 Complex delta = (x2 - x1) / (n - 1.0);
519 retval.elem (0) = x1; 519 retval.elem (0) = x1;
520 for (int i = 1; i < n-1; i++) 520 for (octave_idx_type i = 1; i < n-1; i++)
521 retval.elem (i) = x1 + 1.0 * i * delta; 521 retval.elem (i) = x1 + 1.0 * i * delta;
522 retval.elem (n-1) = x2; 522 retval.elem (n-1) = x2;
523 } 523 }
524 else if (n == 1) 524 else if (n == 1)
525 { 525 {