Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 161
5df5258332d1
parent 154
03f34ba55685
equal deleted inserted replaced
160:e7920e205785 161:5df5258332d1
54 // { 54 // {
55 // self.make_origin_generator().make_origin() 55 // self.make_origin_generator().make_origin()
56 // } 56 // }
57 57
58 /// Return a similar zero as `x`. 58 /// Return a similar zero as `x`.
59 fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::Owned { 59 fn similar_origin_inst<I: Instance<Self::Owned>>(x: I) -> Self::Owned {
60 x.eval(|xr| xr.similar_origin()) 60 x.eval(|xr| xr.similar_origin())
61 } 61 }
62 } 62 }
63 63
64 /// Efficient in-place summation. 64 /// Efficient in-place summation.
75 + for<'b> SubAssign<&'b Self> 75 + for<'b> SubAssign<&'b Self>
76 where 76 where
77 X: Space, 77 X: Space,
78 { 78 {
79 /// Computes `y = βy + αx`, where `y` is `Self`. 79 /// Computes `y = βy + αx`, where `y` is `Self`.
80 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field); 80 fn axpy<I: Instance<X::OwnedSpace>>(&mut self, α: Self::Field, x: I, β: Self::Field);
81 81
82 /// Copies `x` to `self`. 82 /// Copies `x` to `self`.
83 fn copy_from<I: Instance<X>>(&mut self, x: I) { 83 fn copy_from<I: Instance<X::OwnedSpace>>(&mut self, x: I) {
84 self.axpy(1.0, x, 0.0) 84 self.axpy(1.0, x, 0.0)
85 } 85 }
86 86
87 /// Computes `y = αx`, where `y` is `Self`. 87 /// Computes `y = αx`, where `y` is `Self`.
88 fn scale_from<I: Instance<X>>(&mut self, α: Self::Field, x: I) { 88 fn scale_from<I: Instance<X::OwnedSpace>>(&mut self, α: Self::Field, x: I) {
89 self.axpy(α, x, 0.0) 89 self.axpy(α, x, 0.0)
90 } 90 }
91 91
92 /// Set self to zero. 92 /// Set self to zero.
93 fn set_zero(&mut self); 93 fn set_zero(&mut self);
98 98
99 /// Efficient in-place application for [`Linear`] operators. 99 /// Efficient in-place application for [`Linear`] operators.
100 #[replace_float_literals(F::cast_from(literal))] 100 #[replace_float_literals(F::cast_from(literal))]
101 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> { 101 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> {
102 /// Computes `y = αAx + βy`, where `A` is `Self`. 102 /// Computes `y = αAx + βy`, where `A` is `Self`.
103 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F); 103 fn gemv<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, α: F, x: I, β: F);
104 104
105 #[inline] 105 #[inline]
106 /// Computes `y = Ax`, where `A` is `Self` 106 /// Computes `y = Ax`, where `A` is `Self`
107 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) { 107 fn apply_mut<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, x: I) {
108 self.gemv(y, 1.0, x, 0.0) 108 self.gemv(y, 1.0, x, 0.0)
109 } 109 }
110 110
111 #[inline] 111 #[inline]
112 /// Computes `y += Ax`, where `A` is `Self` 112 /// Computes `y += Ax`, where `A` is `Self`
113 fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) { 113 fn apply_add<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, x: I) {
114 self.gemv(y, 1.0, x, 1.0) 114 self.gemv(y, 1.0, x, 1.0)
115 } 115 }
116 } 116 }
117 117
118 /// Bounded linear operators 118 /// Bounded linear operators
197 } 197 }
198 198
199 impl<X: Space> Mapping<X> for IdOp<X> { 199 impl<X: Space> Mapping<X> for IdOp<X> {
200 type Codomain = X::OwnedSpace; 200 type Codomain = X::OwnedSpace;
201 201
202 fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain { 202 fn apply<I: Instance<X::OwnedSpace>>(&self, x: I) -> Self::Codomain {
203 x.into_owned() 203 x.into_owned()
204 } 204 }
205 } 205 }
206 206
207 impl<X: Space> Linear<X> for IdOp<X> {} 207 impl<X: Space> Linear<X> for IdOp<X> {}
211 where 211 where
212 Y: AXPY<X, Field = F>, 212 Y: AXPY<X, Field = F>,
213 X: Space, 213 X: Space,
214 { 214 {
215 // Computes `y = αAx + βy`, where `A` is `Self`. 215 // Computes `y = αAx + βy`, where `A` is `Self`.
216 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) { 216 fn gemv<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, α: F, x: I, β: F) {
217 y.axpy(α, x, β) 217 y.axpy(α, x, β)
218 } 218 }
219 219
220 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) { 220 fn apply_mut<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, x: I) {
221 y.copy_from(x); 221 y.copy_from(x);
222 } 222 }
223 } 223 }
224 224
225 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> 225 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
262 pub struct SimpleZeroOp; 262 pub struct SimpleZeroOp;
263 263
264 impl<X: VectorSpace> Mapping<X> for SimpleZeroOp { 264 impl<X: VectorSpace> Mapping<X> for SimpleZeroOp {
265 type Codomain = X::Owned; 265 type Codomain = X::Owned;
266 266
267 fn apply<I: Instance<X>>(&self, x: I) -> X::Owned { 267 fn apply<I: Instance<X::OwnedSpace>>(&self, x: I) -> X::Owned {
268 X::similar_origin_inst(x) 268 X::similar_origin_inst(x)
269 } 269 }
270 } 270 }
271 271
272 impl<X: VectorSpace> Linear<X> for SimpleZeroOp {} 272 impl<X: VectorSpace> Linear<X> for SimpleZeroOp {}
277 F: Num, 277 F: Num,
278 Y: AXPY<Field = F>, 278 Y: AXPY<Field = F>,
279 X: VectorSpace<Field = F> + Instance<X>, 279 X: VectorSpace<Field = F> + Instance<X>,
280 { 280 {
281 // Computes `y = αAx + βy`, where `A` is `Self`. 281 // Computes `y = αAx + βy`, where `A` is `Self`.
282 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) { 282 fn gemv<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
283 *y *= β; 283 *y *= β;
284 } 284 }
285 285
286 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) { 286 fn apply_mut<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, _x: I) {
287 y.set_zero(); 287 y.set_zero();
288 } 288 }
289 } 289 }
290 290
291 impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp 291 impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp
412 F: Float, 412 F: Float,
413 OY: OriginGenerator<Y>, 413 OY: OriginGenerator<Y>,
414 { 414 {
415 type Codomain = Y::Owned; 415 type Codomain = Y::Owned;
416 416
417 fn apply<I: Instance<X>>(&self, _x: I) -> Y::Owned { 417 fn apply<I: Instance<X::OwnedSpace>>(&self, _x: I) -> Y::Owned {
418 self.codomain_origin_generator.origin() 418 self.codomain_origin_generator.origin()
419 } 419 }
420 } 420 }
421 421
422 impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F> 422 impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F>
435 Y: AXPY<Field = F, Owned = Y>, 435 Y: AXPY<Field = F, Owned = Y>,
436 F: Float, 436 F: Float,
437 OY: OriginGenerator<Y>, 437 OY: OriginGenerator<Y>,
438 { 438 {
439 // Computes `y = αAx + βy`, where `A` is `Self`. 439 // Computes `y = αAx + βy`, where `A` is `Self`.
440 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) { 440 fn gemv<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
441 *y *= β; 441 *y *= β;
442 } 442 }
443 443
444 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) { 444 fn apply_mut<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, _x: I) {
445 y.set_zero(); 445 y.set_zero();
446 } 446 }
447 } 447 }
448 448
449 impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F> 449 impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F>
502 F: Num, 502 F: Num,
503 X: Space, 503 X: Space,
504 T: Linear<X>, 504 T: Linear<X>,
505 S: GEMV<F, T::Codomain, Y>, 505 S: GEMV<F, T::Codomain, Y>,
506 { 506 {
507 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) { 507 fn gemv<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, α: F, x: I, β: F) {
508 self.outer.gemv(y, α, self.inner.apply(x), β) 508 self.outer.gemv(y, α, self.inner.apply(x), β)
509 } 509 }
510 510
511 /// Computes `y = Ax`, where `A` is `Self` 511 /// Computes `y = Ax`, where `A` is `Self`
512 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) { 512 fn apply_mut<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, x: I) {
513 self.outer.apply_mut(y, self.inner.apply(x)) 513 self.outer.apply_mut(y, self.inner.apply(x))
514 } 514 }
515 515
516 /// Computes `y += Ax`, where `A` is `Self` 516 /// Computes `y += Ax`, where `A` is `Self`
517 fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) { 517 fn apply_add<I: Instance<X::OwnedSpace>>(&self, y: &mut Y, x: I) {
518 self.outer.apply_add(y, self.inner.apply(x)) 518 self.outer.apply_add(y, self.inner.apply(x))
519 } 519 }
520 } 520 }
521 521
522 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp> 522 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
549 S::Codomain: Add<T::Codomain>, 549 S::Codomain: Add<T::Codomain>,
550 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace, 550 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
551 { 551 {
552 type Codomain = <S::Codomain as Add<T::Codomain>>::Output; 552 type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
553 553
554 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain { 554 fn apply<I: Instance<Pair<A::OwnedSpace, B::OwnedSpace>>>(&self, x: I) -> Self::Codomain {
555 x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b)) 555 x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b))
556 } 556 }
557 } 557 }
558 558
559 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T> 559 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
574 S: GEMV<F, U, Y>, 574 S: GEMV<F, U, Y>,
575 T: GEMV<F, V, Y>, 575 T: GEMV<F, V, Y>,
576 F: Num, 576 F: Num,
577 Self: Linear<Pair<U, V>, Codomain = Y>, 577 Self: Linear<Pair<U, V>, Codomain = Y>,
578 { 578 {
579 fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Y, α: F, x: I, β: F) { 579 fn gemv<I: Instance<Pair<U::OwnedSpace, V::OwnedSpace>>>(&self, y: &mut Y, α: F, x: I, β: F) {
580 x.eval_decompose(|Pair(u, v)| { 580 x.eval_decompose(|Pair(u, v)| {
581 self.0.gemv(y, α, u, β); 581 self.0.gemv(y, α, u, β);
582 self.1.gemv(y, α, v, F::ONE); 582 self.1.gemv(y, α, v, F::ONE);
583 }) 583 })
584 } 584 }
585 585
586 fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) { 586 fn apply_mut<I: Instance<Pair<U::OwnedSpace, V::OwnedSpace>>>(&self, y: &mut Y, x: I) {
587 x.eval_decompose(|Pair(u, v)| { 587 x.eval_decompose(|Pair(u, v)| {
588 self.0.apply_mut(y, u); 588 self.0.apply_mut(y, u);
589 self.1.apply_add(y, v); 589 self.1.apply_add(y, v);
590 }) 590 })
591 } 591 }
592 592
593 /// Computes `y += Ax`, where `A` is `Self` 593 /// Computes `y += Ax`, where `A` is `Self`
594 fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) { 594 fn apply_add<I: Instance<Pair<U::OwnedSpace, V::OwnedSpace>>>(&self, y: &mut Y, x: I) {
595 x.eval_decompose(|Pair(u, v)| { 595 x.eval_decompose(|Pair(u, v)| {
596 self.0.apply_add(y, u); 596 self.0.apply_add(y, u);
597 self.1.apply_add(y, v); 597 self.1.apply_add(y, v);
598 }) 598 })
599 } 599 }
609 S: Mapping<A>, 609 S: Mapping<A>,
610 T: Mapping<A>, 610 T: Mapping<A>,
611 { 611 {
612 type Codomain = Pair<S::Codomain, T::Codomain>; 612 type Codomain = Pair<S::Codomain, T::Codomain>;
613 613
614 fn apply<I: Instance<A>>(&self, a: I) -> Self::Codomain { 614 fn apply<I: Instance<A::OwnedSpace>>(&self, a: I) -> Self::Codomain {
615 Pair(a.eval_ref_decompose(|r| self.0.apply(r)), self.1.apply(a)) 615 Pair(a.eval_ref_decompose(|r| self.0.apply(r)), self.1.apply(a))
616 } 616 }
617 } 617 }
618 618
619 impl<A, S, T> Linear<A> for ColOp<S, T> 619 impl<A, S, T> Linear<A> for ColOp<S, T>
630 S: GEMV<F, X, A>, 630 S: GEMV<F, X, A>,
631 T: GEMV<F, X, B>, 631 T: GEMV<F, X, B>,
632 F: Num, 632 F: Num,
633 Self: Linear<X, Codomain = Pair<A, B>>, 633 Self: Linear<X, Codomain = Pair<A, B>>,
634 { 634 {
635 fn gemv<I: Instance<X>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) { 635 fn gemv<I: Instance<X::OwnedSpace>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
636 x.eval_ref_decompose(|r| self.0.gemv(&mut y.0, α, r, β)); 636 x.eval_ref_decompose(|r| self.0.gemv(&mut y.0, α, r, β));
637 self.1.gemv(&mut y.1, α, x, β); 637 self.1.gemv(&mut y.1, α, x, β);
638 } 638 }
639 639
640 fn apply_mut<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) { 640 fn apply_mut<I: Instance<X::OwnedSpace>>(&self, y: &mut Pair<A, B>, x: I) {
641 x.eval_ref_decompose(|r| self.0.apply_mut(&mut y.0, r)); 641 x.eval_ref_decompose(|r| self.0.apply_mut(&mut y.0, r));
642 self.1.apply_mut(&mut y.1, x); 642 self.1.apply_mut(&mut y.1, x);
643 } 643 }
644 644
645 /// Computes `y += Ax`, where `A` is `Self` 645 /// Computes `y += Ax`, where `A` is `Self`
646 fn apply_add<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) { 646 fn apply_add<I: Instance<X::OwnedSpace>>(&self, y: &mut Pair<A, B>, x: I) {
647 x.eval_ref_decompose(|r| self.0.apply_add(&mut y.0, r)); 647 x.eval_ref_decompose(|r| self.0.apply_add(&mut y.0, r));
648 self.1.apply_add(&mut y.1, x); 648 self.1.apply_add(&mut y.1, x);
649 } 649 }
650 } 650 }
651 651
753 S: Mapping<A>, 753 S: Mapping<A>,
754 T: Mapping<B>, 754 T: Mapping<B>,
755 { 755 {
756 type Codomain = Pair<S::Codomain, T::Codomain>; 756 type Codomain = Pair<S::Codomain, T::Codomain>;
757 757
758 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain { 758 fn apply<I: Instance<Pair<A::OwnedSpace, B::OwnedSpace>>>(&self, x: I) -> Self::Codomain {
759 x.eval_decompose(|Pair(a, b)| Pair(self.0.apply(a), self.1.apply(b))) 759 x.eval_decompose(|Pair(a, b)| Pair(self.0.apply(a), self.1.apply(b)))
760 } 760 }
761 } 761 }
762 762
763 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T> 763 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
778 S: GEMV<F, U, A>, 778 S: GEMV<F, U, A>,
779 T: GEMV<F, V, B>, 779 T: GEMV<F, V, B>,
780 F: Num, 780 F: Num,
781 Self: Linear<Pair<U, V>, Codomain = Pair<A, B>>, 781 Self: Linear<Pair<U, V>, Codomain = Pair<A, B>>,
782 { 782 {
783 fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) { 783 fn gemv<I: Instance<Pair<U::OwnedSpace, V::OwnedSpace>>>(
784 &self,
785 y: &mut Pair<A, B>,
786 α: F,
787 x: I,
788 β: F,
789 ) {
784 x.eval_decompose(|Pair(u, v)| { 790 x.eval_decompose(|Pair(u, v)| {
785 self.0.gemv(&mut y.0, α, u, β); 791 self.0.gemv(&mut y.0, α, u, β);
786 self.1.gemv(&mut y.1, α, v, β); 792 self.1.gemv(&mut y.1, α, v, β);
787 }) 793 })
788 } 794 }
789 795
790 fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) { 796 fn apply_mut<I: Instance<Pair<U::OwnedSpace, V::OwnedSpace>>>(&self, y: &mut Pair<A, B>, x: I) {
791 x.eval_decompose(|Pair(u, v)| { 797 x.eval_decompose(|Pair(u, v)| {
792 self.0.apply_mut(&mut y.0, u); 798 self.0.apply_mut(&mut y.0, u);
793 self.1.apply_mut(&mut y.1, v); 799 self.1.apply_mut(&mut y.1, v);
794 }) 800 })
795 } 801 }
796 802
797 /// Computes `y += Ax`, where `A` is `Self` 803 /// Computes `y += Ax`, where `A` is `Self`
798 fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) { 804 fn apply_add<I: Instance<Pair<U::OwnedSpace, V::OwnedSpace>>>(&self, y: &mut Pair<A, B>, x: I) {
799 x.eval_decompose(|Pair(u, v)| { 805 x.eval_decompose(|Pair(u, v)| {
800 self.0.apply_add(&mut y.0, u); 806 self.0.apply_add(&mut y.0, u);
801 self.1.apply_add(&mut y.1, v); 807 self.1.apply_add(&mut y.1, v);
802 }) 808 })
803 } 809 }
923 <Domain as Mul<F>>::Output: ClosedSpace, 929 <Domain as Mul<F>>::Output: ClosedSpace,
924 { 930 {
925 type Codomain = <Domain as Mul<F>>::Output; 931 type Codomain = <Domain as Mul<F>>::Output;
926 932
927 /// Compute the value of `self` at `x`. 933 /// Compute the value of `self` at `x`.
928 fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain { 934 fn apply<I: Instance<Domain::OwnedSpace>>(&self, x: I) -> Self::Codomain {
929 x.own() * self.0 935 x.own() * self.0
930 } 936 }
931 } 937 }
932 938
933 impl<Domain, F> Linear<Domain> for Scaled<F> 939 impl<Domain, F> Linear<Domain> for Scaled<F>

Mercurial Repository: alg_tools

mercurial