Mercurial > repos > alg_tools / file diff

--- a/src/linops.rs	Sun Nov 10 09:02:57 2024 -0500
+++ b/src/linops.rs	Tue Dec 31 09:12:43 2024 -0500
@@ -10,7 +10,8 @@
 
 /// Trait for linear operators on `X`.
 pub trait Linear<X> : Apply<X, Output=Self::Codomain>
-                      + for<'a> Apply<&'a X, Output=Self::Codomain> {
+                      + for<'a> Apply<&'a X, Output=Self::Codomain>
+                      + HasScalarField {
     type Codomain;
 }
 
@@ -18,32 +19,32 @@
 #[replace_float_literals(F::cast_from(literal))]
 pub trait AXPY<F : Num, X = Self> {
     /// Computes  `y = βy + αx`, where `y` is `Self`.
-    fn axpy(&mut self, α : F, x : &X, β : F);
+    fn axpy(&mut self, α : F, x : X, β : F);
 
     /// Copies `x` to `self`.
-    fn copy_from(&mut self, x : &X) {
+    fn copy_from(&mut self, x : X) {
         self.axpy(1.0, x, 0.0)
     }
 
     /// Computes  `y = αx`, where `y` is `Self`.
-    fn scale_from(&mut self, α : F, x : &X) {
+    fn scale_from(&mut self, α : F, x : X) {
         self.axpy(α, x, 0.0)
     }
 }
 
 /// Efficient in-place application for [`Linear`] operators.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> {
+pub trait GEMV<F : Num, X, Y = <Self as Apply<X>>::Output> : Apply<X> {
     /// Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv(&self, y : &mut Y, α : F, x : &X, β : F);
+    fn gemv(&self, y : &mut Y, α : F, x : X, β : F);
 
     /// Computes `y = Ax`, where `A` is `Self`
-    fn apply_mut(&self, y : &mut Y, x : &X){
+    fn apply_mut(&self, y : &mut Y, x : X){
         self.gemv(y, 1.0, x, 0.0)
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add(&self, y : &mut Y, x : &X){
+    fn apply_add(&self, y : &mut Y, x : X){
         self.gemv(y, 1.0, x, 1.0)
     }
 }
@@ -51,11 +52,10 @@
 
 /// Bounded linear operators
 pub trait BoundedLinear<X> : Linear<X> {
-    type FloatType : Float;
     /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
     /// This is not expected to be the norm, just any bound on it that can be
     /// reasonably implemented.
-    fn opnorm_bound(&self) -> Self::FloatType;
+    fn opnorm_bound(&self) -> Self::Field;
 }
 
 // Linear operator application into mutable target. The [`AsRef`] bound
@@ -70,7 +70,7 @@
 
 /// Trait for forming the adjoint operator of `Self`.
 pub trait Adjointable<X,Yʹ> : Linear<X> {
-    type AdjointCodomain;
+    type AdjointCodomain : HasScalarField<Field=Self::Field>;
     type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
 
     /// Form the adjoint operator of `self`.
@@ -89,7 +89,7 @@
 /// domain of the adjointed operator. `Self::Preadjoint` should be
 /// [`Adjointable`]`<'a,Ypre,X>`.
 pub trait Preadjointable<X,Ypre> : Linear<X> {
-    type PreadjointCodomain;
+    type PreadjointCodomain : HasScalarField<Field=Self::Field>;
     type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a;
 
     /// Form the preadjoint operator of `self`.
@@ -100,55 +100,49 @@
 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
 
-/// The identity operator
+/// The identity operator on `X` with scalar field `F`.
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct IdOp<X> (PhantomData<X>);
+pub struct IdOp<X, F : Num> (PhantomData<(X, F)>);
 
-impl<X> IdOp<X> {
-    fn new() -> IdOp<X> { IdOp(PhantomData) }
+impl<X, F : Num> HasScalarField for IdOp<X, F> {
+    type Field = F;
 }
 
-impl<X> Apply<X> for IdOp<X> {
+impl<X, F : Num> IdOp<X, F> {
+    fn new() -> IdOp<X, F> { IdOp(PhantomData) }
+}
+
+impl<X, F : Num, T : CloneIfNeeded<X>> Apply<T> for IdOp<X, F> {
     type Output = X;
 
-    fn apply(&self, x : X) -> X {
-        x
+    fn apply(&self, x : T) -> X {
+        x.clone_if_needed()
     }
 }
 
-impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone {
-    type Output = X;
-    
-    fn apply(&self, x : &'a X) -> X {
-        x.clone()
-    }
-}
-
-impl<X> Linear<X> for IdOp<X> where X : Clone {
+impl<X : Clone, F : Num> Linear<X> for IdOp<X, F> {
     type Codomain = X;
 }
 
-
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone {
+impl<F : Num, X : Clone, Y> GEMV<F, X, Y> for IdOp<X, F> where Y : AXPY<F, X> {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) {
+    fn gemv(&self, y : &mut Y, α : F, x : X, β : F) {
         y.axpy(α, x, β)
     }
 
-    fn apply_mut(&self, y : &mut Y, x : &X){
+    fn apply_mut(&self, y : &mut Y, x : X){
         y.copy_from(x);
     }
 }
 
-impl<X> BoundedLinear<X> for IdOp<X> where X : Clone {
-    type FloatType = float;
-    fn opnorm_bound(&self) -> float { 1.0 }
+impl<X, F : Num> BoundedLinear<X> for IdOp<X, F> where X : Clone {
+    fn opnorm_bound(&self) -> F { F::ONE }
 }
 
-impl<X> Adjointable<X,X> for IdOp<X> where X : Clone {
+impl<X, F : Num> Adjointable<X,X> for IdOp<X, F> where X : Clone + HasScalarField<Field=F> {
     type AdjointCodomain=X;
-    type Adjoint<'a> = IdOp<X> where X : 'a;
+    type Adjoint<'a> = IdOp<X, F> where X : 'a;
     fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
 }