Mercurial > octave-nkf

comparison src/sparse-xpow.cc @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children 23b37da9fd5b
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 /*
2
3 Copyright (C) 2004 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5
6 Octave is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version.
10
11 Octave is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; see the file COPYING. If not, write to the Free
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19
20 */
21
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include <cassert>
27 #include <climits>
28
29 #include "Array-util.h"
30 #include "oct-cmplx.h"
31 #include "quit.h"
32
33 #include "error.h"
34 #include "oct-obj.h"
35 #include "utils.h"
36
37 #include "dSparse.h"
38 #include "CSparse.h"
39 #include "ov-re-sparse.h"
40 #include "ov-cx-sparse.h"
41 #include "sparse-xpow.h"
42
43 static inline int
44 xisint (double x)
45 {
46 return (D_NINT (x) == x
47 && ((x >= 0 && x < INT_MAX)
48 || (x <= 0 && x > INT_MIN)));
49 }
50
51
52 // Safer pow functions. Only two make sense for sparse matrices, the
53 // others should all promote to full matrices.
54
55 octave_value
56 xpow (const SparseMatrix& a, double b)
57 {
58 octave_value retval;
59
60 int nr = a.rows ();
61 int nc = a.cols ();
62
63 if (nr == 0 || nc == 0 || nr != nc)
64 error ("for A^b, A must be square");
65 else
66 {
67 if (static_cast<int> (b) == b)
68 {
69 int btmp = static_cast<int> (b);
70 if (btmp == 0)
71 {
72 SparseMatrix tmp = SparseMatrix (nr, nr, nr);
73 for (int i = 0; i < nr; i++)
74 {
75 tmp.data (i) = 1.0;
76 tmp.ridx (i) = i;
77 }
78 for (int i = 0; i < nr + 1; i++)
79 tmp.cidx (i) = i;
80
81 retval = tmp;
82 }
83 else
84 {
85 SparseMatrix atmp;
86 if (btmp < 0)
87 {
88 btmp = -btmp;
89
90 int info;
91 double rcond = 0.0;
92
93 atmp = a.inverse (info, rcond, 1);
94
95 if (info == -1)
96 warning ("inverse: matrix singular to machine\
97 precision, rcond = %g", rcond);
98 }
99 else
100 atmp = a;
101
102 SparseMatrix result (atmp);
103
104 btmp--;
105
106 while (btmp > 0)
107 {
108 if (btmp & 1)
109 result = result * atmp;
110
111 btmp >>= 1;
112
113 if (btmp > 0)
114 atmp = atmp * atmp;
115 }
116
117 retval = result;
118 }
119 }
120 else
121 error ("use full(a) ^ full(b)");
122 }
123
124 return retval;
125 }
126
127 octave_value
128 xpow (const SparseComplexMatrix& a, double b)
129 {
130 octave_value retval;
131
132 int nr = a.rows ();
133 int nc = a.cols ();
134
135 if (nr == 0 || nc == 0 || nr != nc)
136 error ("for A^b, A must be square");
137 else
138 {
139 if (static_cast<int> (b) == b)
140 {
141 int btmp = static_cast<int> (b);
142 if (btmp == 0)
143 {
144 SparseMatrix tmp = SparseMatrix (nr, nr, nr);
145 for (int i = 0; i < nr; i++)
146 {
147 tmp.data (i) = 1.0;
148 tmp.ridx (i) = i;
149 }
150 for (int i = 0; i < nr + 1; i++)
151 tmp.cidx (i) = i;
152
153 retval = tmp;
154 }
155 else
156 {
157 SparseComplexMatrix atmp;
158 if (btmp < 0)
159 {
160 btmp = -btmp;
161
162 int info;
163 double rcond = 0.0;
164
165 atmp = a.inverse (info, rcond, 1);
166
167 if (info == -1)
168 warning ("inverse: matrix singular to machine\
169 precision, rcond = %g", rcond);
170 }
171 else
172 atmp = a;
173
174 SparseComplexMatrix result (atmp);
175
176 btmp--;
177
178 while (btmp > 0)
179 {
180 if (btmp & 1)
181 result = result * atmp;
182
183 btmp >>= 1;
184
185 if (btmp > 0)
186 atmp = atmp * atmp;
187 }
188
189 retval = result;
190 }
191 }
192 else
193 error ("use full(a) ^ full(b)");
194 }
195
196 return retval;
197 }
198
199 // Safer pow functions that work elementwise for matrices.
200 //
201 // op2 \ op1: s m cs cm
202 // +-- +---+---+----+----+
203 // scalar | | * | 3 | * | 9 |
204 // +---+---+----+----+
205 // matrix | 1 | 4 | 7 | 10 |
206 // +---+---+----+----+
207 // complex_scalar | * | 5 | * | 11 |
208 // +---+---+----+----+
209 // complex_matrix | 2 | 6 | 8 | 12 |
210 // +---+---+----+----+
211 //
212 // * -> not needed.
213
214 // XXX FIXME XXX -- these functions need to be fixed so that things
215 // like
216 //
217 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b
218 //
219 // and
220 //
221 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end
222 //
223 // produce identical results. Also, it would be nice if -1^0.5
224 // produced a pure imaginary result instead of a complex number with a
225 // small real part. But perhaps that's really a problem with the math
226 // library...
227
228 // -*- 1 -*-
229 octave_value
230 elem_xpow (double a, const SparseMatrix& b)
231 {
232 octave_value retval;
233
234 int nr = b.rows ();
235 int nc = b.cols ();
236
237 double d1, d2;
238
239 if (a < 0.0 && ! b.all_integers (d1, d2))
240 {
241 Complex atmp (a);
242 ComplexMatrix result (nr, nc);
243
244 for (int j = 0; j < nc; j++)
245 {
246 for (int i = 0; i < nr; i++)
247 {
248 OCTAVE_QUIT;
249 result (i, j) = pow (atmp, b(i,j));
250 }
251 }
252
253 retval = result;
254 }
255 else
256 {
257 Matrix result (nr, nc);
258
259 for (int j = 0; j < nc; j++)
260 {
261 for (int i = 0; i < nr; i++)
262 {
263 OCTAVE_QUIT;
264 result (i, j) = pow (a, b(i,j));
265 }
266 }
267
268 retval = result;
269 }
270
271 return retval;
272 }
273
274 // -*- 2 -*-
275 octave_value
276 elem_xpow (double a, const SparseComplexMatrix& b)
277 {
278 int nr = b.rows ();
279 int nc = b.cols ();
280
281 Complex atmp (a);
282 ComplexMatrix result (nr, nc);
283
284 for (int j = 0; j < nc; j++)
285 {
286 for (int i = 0; i < nr; i++)
287 {
288 OCTAVE_QUIT;
289 result (i, j) = pow (atmp, b(i,j));
290 }
291 }
292
293 return result;
294 }
295
296 // -*- 3 -*-
297 octave_value
298 elem_xpow (const SparseMatrix& a, double b)
299 {
300 // XXX FIXME XXX What should a .^ 0 give?? Matlab gives a
301 // sparse matrix with same structure as a, which is strictly
302 // incorrect. Keep compatiability.
303
304 octave_value retval;
305
306 int nz = a.nnz ();
307
308 if (b <= 0.0)
309 {
310 int nr = a.rows ();
311 int nc = a.cols ();
312
313 if (static_cast<int> (b) != b && a.any_element_is_negative ())
314 {
315 ComplexMatrix result (nr, nc, Complex (pow (0.0, b)));
316
317 // XXX FIXME XXX -- avoid apparent GNU libm bug by
318 // converting A and B to complex instead of just A.
319 Complex btmp (b);
320
321 for (int j = 0; j < nc; j++)
322 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
323 {
324 OCTAVE_QUIT;
325
326 Complex atmp (a.data (i));
327
328 result (a.ridx(i), j) = pow (atmp, btmp);
329 }
330
331 retval = octave_value (result);
332 }
333 else
334 {
335 Matrix result (nr, nc, (pow (0.0, b)));
336
337 for (int j = 0; j < nc; j++)
338 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
339 {
340 OCTAVE_QUIT;
341 result (a.ridx(i), j) = pow (a.data (i), b);
342 }
343
344 retval = octave_value (result);
345 }
346 }
347 else if (static_cast<int> (b) != b && a.any_element_is_negative ())
348 {
349 SparseComplexMatrix result (a);
350
351 for (int i = 0; i < nz; i++)
352 {
353 OCTAVE_QUIT;
354
355 // XXX FIXME XXX -- avoid apparent GNU libm bug by
356 // converting A and B to complex instead of just A.
357
358 Complex atmp (a.data (i));
359 Complex btmp (b);
360
361 result.data (i) = pow (atmp, btmp);
362 }
363
364 result.maybe_compress (true);
365
366 retval = result;
367 }
368 else
369 {
370 SparseMatrix result (a);
371
372 for (int i = 0; i < nz; i++)
373 {
374 OCTAVE_QUIT;
375 result.data (i) = pow (a.data (i), b);
376 }
377
378 result.maybe_compress (true);
379
380 retval = result;
381 }
382
383 return retval;
384 }
385
386 // -*- 4 -*-
387 octave_value
388 elem_xpow (const SparseMatrix& a, const SparseMatrix& b)
389 {
390 octave_value retval;
391
392 int nr = a.rows ();
393 int nc = a.cols ();
394
395 int b_nr = b.rows ();
396 int b_nc = b.cols ();
397
398 if (nr != b_nr || nc != b_nc)
399 {
400 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
401 return octave_value ();
402 }
403
404 int convert_to_complex = 0;
405 for (int j = 0; j < nc; j++)
406 for (int i = 0; i < nr; i++)
407 {
408 OCTAVE_QUIT;
409 double atmp = a (i, j);
410 double btmp = b (i, j);
411 if (atmp < 0.0 && static_cast<int> (btmp) != btmp)
412 {
413 convert_to_complex = 1;
414 goto done;
415 }
416 }
417
418 done:
419
420 int nel = 0;
421 for (int j = 0; j < nc; j++)
422 for (int i = 0; i < nr; i++)
423 if (!(a.elem (i, j) == 0. && b.elem (i, j) != 0.))
424 nel++;
425
426 if (convert_to_complex)
427 {
428 SparseComplexMatrix complex_result (nr, nc, nel);
429
430 int ii = 0;
431 complex_result.cidx(0) = 0;
432 for (int j = 0; j < nc; j++)
433 {
434 for (int i = 0; i < nr; i++)
435 {
436 OCTAVE_QUIT;
437 Complex atmp (a (i, j));
438 Complex btmp (b (i, j));
439 Complex tmp = pow (atmp, btmp);
440 if (tmp != 0.)
441 {
442 complex_result.data (ii) = tmp;
443 complex_result.ridx (ii++) = i;
444 }
445 }
446 complex_result.cidx (j+1) = ii;
447 }
448 complex_result.maybe_compress ();
449
450 retval = complex_result;
451 }
452 else
453 {
454 SparseMatrix result (nr, nc, nel);
455 int ii = 0;
456
457 result.cidx (0) = 0;
458 for (int j = 0; j < nc; j++)
459 {
460 for (int i = 0; i < nr; i++)
461 {
462 OCTAVE_QUIT;
463 double tmp = pow (a (i, j), b (i, j));
464 if (tmp != 0.)
465 {
466 result.data (ii) = tmp;
467 result.ridx (ii++) = i;
468 }
469 }
470 result.cidx (j+1) = ii;
471 }
472
473 result.maybe_compress ();
474
475 retval = result;
476 }
477
478 return retval;
479 }
480
481 // -*- 5 -*-
482 octave_value
483 elem_xpow (const SparseMatrix& a, const Complex& b)
484 {
485 octave_value retval;
486
487 if (b == 0.0)
488 // Can this case ever happen, due to automatic retyping with maybe_mutate?
489 retval = octave_value (NDArray (a.dims (), 1));
490 else
491 {
492 int nz = a.nnz ();
493 SparseComplexMatrix result (a);
494
495 for (int i = 0; i < nz; i++)
496 {
497 OCTAVE_QUIT;
498 result.data (i) = pow (Complex (a.data (i)), b);
499 }
500
501 result.maybe_compress (true);
502
503 retval = result;
504 }
505
506 return retval;
507 }
508
509 // -*- 6 -*-
510 octave_value
511 elem_xpow (const SparseMatrix& a, const SparseComplexMatrix& b)
512 {
513 int nr = a.rows ();
514 int nc = a.cols ();
515
516 int b_nr = b.rows ();
517 int b_nc = b.cols ();
518
519 if (nr != b_nr || nc != b_nc)
520 {
521 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
522 return octave_value ();
523 }
524
525 int nel = 0;
526 for (int j = 0; j < nc; j++)
527 for (int i = 0; i < nr; i++)
528 if (!(a.elem (i, j) == 0. && b.elem (i, j) != 0.))
529 nel++;
530
531 SparseComplexMatrix result (nr, nc, nel);
532 int ii = 0;
533
534 result.cidx(0) = 0;
535 for (int j = 0; j < nc; j++)
536 {
537 for (int i = 0; i < nr; i++)
538 {
539 OCTAVE_QUIT;
540 Complex tmp = pow (Complex (a (i, j)), b (i, j));
541 if (tmp != 0.)
542 {
543 result.data (ii) = tmp;
544 result.ridx (ii++) = i;
545 }
546 }
547 result.cidx (j+1) = ii;
548 }
549
550 result.maybe_compress ();
551
552 return result;
553 }
554
555 // -*- 7 -*-
556 octave_value
557 elem_xpow (const Complex& a, const SparseMatrix& b)
558 {
559 int nr = b.rows ();
560 int nc = b.cols ();
561
562 ComplexMatrix result (nr, nc);
563
564 for (int j = 0; j < nc; j++)
565 {
566 for (int i = 0; i < nr; i++)
567 {
568 OCTAVE_QUIT;
569 double btmp = b (i, j);
570 if (xisint (btmp))
571 result (i, j) = pow (a, static_cast<int> (btmp));
572 else
573 result (i, j) = pow (a, btmp);
574 }
575 }
576
577 return result;
578 }
579
580 // -*- 8 -*-
581 octave_value
582 elem_xpow (const Complex& a, const SparseComplexMatrix& b)
583 {
584 int nr = b.rows ();
585 int nc = b.cols ();
586
587 ComplexMatrix result (nr, nc);
588 for (int j = 0; j < nc; j++)
589 for (int i = 0; i < nr; i++)
590 {
591 OCTAVE_QUIT;
592 result (i, j) = pow (a, b (i, j));
593 }
594
595 return result;
596 }
597
598 // -*- 9 -*-
599 octave_value
600 elem_xpow (const SparseComplexMatrix& a, double b)
601 {
602 octave_value retval;
603
604 if (b <= 0)
605 {
606 int nr = a.rows ();
607 int nc = a.cols ();
608
609 ComplexMatrix result (nr, nc, Complex (pow (0.0, b)));
610
611 if (xisint (b))
612 {
613 for (int j = 0; j < nc; j++)
614 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
615 {
616 OCTAVE_QUIT;
617 result (a.ridx(i), j) =
618 pow (a.data (i), static_cast<int> (b));
619 }
620 }
621 else
622 {
623 for (int j = 0; j < nc; j++)
624 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
625 {
626 OCTAVE_QUIT;
627 result (a.ridx(i), j) = pow (a.data (i), b);
628 }
629 }
630
631 retval = result;
632 }
633 else
634 {
635 int nz = a.nnz ();
636
637 SparseComplexMatrix result (a);
638
639 if (xisint (b))
640 {
641 for (int i = 0; i < nz; i++)
642 {
643 OCTAVE_QUIT;
644 result.data (i) = pow (a.data (i), static_cast<int> (b));
645 }
646 }
647 else
648 {
649 for (int i = 0; i < nz; i++)
650 {
651 OCTAVE_QUIT;
652 result.data (i) = pow (a.data (i), b);
653 }
654 }
655
656 result.maybe_compress (true);
657
658 retval = result;
659 }
660
661 return retval;
662 }
663
664 // -*- 10 -*-
665 octave_value
666 elem_xpow (const SparseComplexMatrix& a, const SparseMatrix& b)
667 {
668 int nr = a.rows ();
669 int nc = a.cols ();
670
671 int b_nr = b.rows ();
672 int b_nc = b.cols ();
673
674 if (nr != b_nr || nc != b_nc)
675 {
676 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
677 return octave_value ();
678 }
679
680 int nel = 0;
681 for (int j = 0; j < nc; j++)
682 for (int i = 0; i < nr; i++)
683 if (!(a.elem (i, j) == 0. && b.elem (i, j) != 0.))
684 nel++;
685
686 SparseComplexMatrix result (nr, nc, nel);
687 int ii = 0;
688
689 result.cidx (0) = 0;
690 for (int j = 0; j < nc; j++)
691 {
692 for (int i = 0; i < nr; i++)
693 {
694 OCTAVE_QUIT;
695 double btmp = b (i, j);
696 Complex tmp;
697
698 if (xisint (btmp))
699 tmp = pow (a (i, j), static_cast<int> (btmp));
700 else
701 tmp = pow (a (i, j), btmp);
702 if (tmp != 0.)
703 {
704 result.data (ii) = tmp;
705 result.ridx (ii++) = i;
706 }
707 }
708 result.cidx (j+1) = ii;
709 }
710
711 result.maybe_compress ();
712
713 return result;
714 }
715
716 // -*- 11 -*-
717 octave_value
718 elem_xpow (const SparseComplexMatrix& a, const Complex& b)
719 {
720 octave_value retval;
721
722 if (b == 0.0)
723 // Can this case ever happen, due to automatic retyping with maybe_mutate?
724 retval = octave_value (NDArray (a.dims (), 1));
725 else
726 {
727
728 int nz = a.nnz ();
729
730 SparseComplexMatrix result (a);
731
732 for (int i = 0; i < nz; i++)
733 {
734 OCTAVE_QUIT;
735 result.data (i) = pow (a.data (i), b);
736 }
737
738 result.maybe_compress (true);
739
740 retval = result;
741 }
742
743 return retval;
744 }
745
746 // -*- 12 -*-
747 octave_value
748 elem_xpow (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
749 {
750 int nr = a.rows ();
751 int nc = a.cols ();
752
753 int b_nr = b.rows ();
754 int b_nc = b.cols ();
755
756 if (nr != b_nr || nc != b_nc)
757 {
758 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
759 return octave_value ();
760 }
761
762 int nel = 0;
763 for (int j = 0; j < nc; j++)
764 for (int i = 0; i < nr; i++)
765 if (!(a.elem (i, j) == 0. && b.elem (i, j) != 0.))
766 nel++;
767
768 SparseComplexMatrix result (nr, nc, nel);
769 int ii = 0;
770
771 result.cidx (0) = 0;
772 for (int j = 0; j < nc; j++)
773 {
774 for (int i = 0; i < nr; i++)
775 {
776 OCTAVE_QUIT;
777 Complex tmp = pow (a (i, j), b (i, j));
778 if (tmp != 0.)
779 {
780 result.data (ii) = tmp;
781 result.ridx (ii++) = i;
782 }
783 }
784 result.cidx (j+1) = ii;
785 }
786 result.maybe_compress (true);
787
788 return result;
789 }
790
791 /*
792 ;;; Local Variables: ***
793 ;;; mode: C++ ***
794 ;;; End: ***
795 */