Extend abstract operator

Cargo.toml

@@ -23,9 +23,10 @@ rand_chacha = "0.3"
 rand_distr = "0.4"
 itertools = "0.14"
 blas = "0.22"
-lapack = "0.19"
+lapack = "0.20.0"
+lapack-sys = { git = "https://github.com/ignacia-fp/lapack-sys", rev = "35ac4ab6811252772d28185b5c3dbfef82188c6b" }
 thiserror = "2"
-serde = "1"
+serde = { version = "1.0", features = ["derive"] }
 bempp-distributed-tools = { git = "https://github.com/bempp/distributed_tools.git", optional = true }
 log = { version = "0.4", optional = true}
 
@@ -34,6 +35,8 @@ pulp = { version = "0.21" }
 bytemuck = "1.16.0"
 
 mpi = { version = "0.8.0", optional = true }
+arpack-ng-sys = "0.2.1"
+csv = "1.1"
 
 
 [dev-dependencies]

examples/X0.csv

@@ -0,0 +1,80 @@
+-3.73162828553711,0.1183907761085212,0.9972123936280015
+0.4741392065971786,-2.541872013614432,2.288913411054258
+3.940068239531302,-1.251132980315346,-1.729893078633118
+-3.287807699592993,-1.128344249213524,-0.7000761269532703
+-3.514151978910112,1.237208598365918,3.58000836673225
+-0.657820584471593,-1.621271592386599,-3.019443025261761
+-1.37250785190388,-2.479524669849622,0.4755592803658795
+3.328494185303328,2.74264332172786,-2.215096324116073
+-2.13184746233437,-3.263924952924406,2.677434362140795
+2.105521797095058,1.845357062105455,1.536662697646119
+3.75348782618712,2.76122803402387,3.413389654147215
+-2.310181737890268,-2.72383332155594,0.8400742233740655
+1.231288275223015,-1.761023838728049,0.5687170718989085
+-0.9862518003133482,-2.348499166190018,-1.462370325782947
+-0.9265255608289218,0.4266346859524193,0.1632971652997064
+2.8115628903337,3.433808918089146,0.795905075764086
+2.535274820720547,0.8380555978924038,0.1763041521793314
+3.084696242490499,3.618117636768116,1.749916756659586
+3.883049756104727,1.71093196848618,1.043796500252091
+0.1059876751619333,3.696043239868189,-0.8039499727774793
+-3.719212487329726,-0.05371873555315343,2.351456155241632
+3.483501752182511,3.298657608461658,-0.1509146092060227
+-0.328705484409042,-3.462605174718437,-3.437381480218432
+0.2728115279335235,-1.595098316574369,0.3406603453289341
+-0.8117075964884117,1.464111786745515,-1.346484093333903
+3.288569360074701,-3.560408560087333,-1.278745533722304
+-0.632021489498932,3.470179165304919,-2.41076134417167
+1.843315969444414,1.940568872898417,-1.80440303491904
+-0.4514604844497248,3.416002828548363,-1.285656261405941
+-1.157073873426262,0.5895659357604153,3.170313416307478
+0.1517090290678524,2.38482581965371,1.066929997860353
+1.158086820456133,0.8268269184468933,-2.352483122099351
+2.006783029110328,-1.638879292395481,3.478458111519607
+-3.468123984338773,-3.098826733801999,2.879858738884896
+0.3728958441880517,-0.004142358084177289,-3.242881351432046
+-0.2412896552927224,-2.548252391669281,-0.7041397292482525
+3.845603469522265,1.109963898001604,-3.841350099378027
+-0.06481540623833748,-2.869489384000321,-2.028642103251894
+3.953758100658055,-0.1202734046980671,-3.965423601249017
+2.849852579235237,-1.6193238151773,1.611044008640642
+0.2173845972207493,2.60248065277665,-1.488452111665494
+0.9812204738192287,-0.6760944845235253,1.355306055375801
+3.566598410830857,1.585783537607742,-3.567760822270386
+-0.7002280650485551,-2.502811154871331,1.6497493647021
+2.68797421006431,-0.2747259827292243,2.443469374825868
+-0.4586418995317452,0.003301566454140925,-2.514274880264145
+2.846034227184734,-2.186601635866459,1.76772272154781
+-3.498456677569798,-0.3922453804275987,-0.272269290303341
+-2.167839772158199,1.329935900048681,-0.8176297070980691
+-3.655464978795243,2.952537116523236,1.28292605219653
+2.870177884668004,-0.9960479238676898,1.326121096970781
+1.556858333698411,-3.402769753034622,1.431510540364358
+3.683283006176307,-2.738710853222885,1.85996961087249
+1.785730441666463,3.44200171746023,-2.077378781275213
+-2.45210004092483,1.397718359766978,-0.7138378967364458
+2.581255717773968,3.986349006400111,-0.3859604668388696
+2.763435882313683,-1.397577587415256,-2.593108132198631
+-3.637749245619159,-1.106240933395653,1.942795135267155
+2.763588799287656,-2.063270198828206,2.903333552445571
+-2.671114353209854,2.813808565406188,3.772967434254785
+3.67149754982554,2.263730463986914,3.701168378315072
+0.910531773313374,0.5378209715222635,-3.925409796480833
+1.20326158429046,-1.893798847461403,-3.790517521423674
+1.266337657996369,2.343047978180254,-1.332798072982384
+-2.531704453713843,-2.582695405167952,-3.925283391058958
+0.1925207349461742,0.9503912963635708,-0.1252200219290147
+-2.263122316179229,-2.994286234093984,0.5235568339384757
+0.2332789975292391,-3.680197315218904,-0.1222738568504367
+-2.538752054865134,-0.2999017615226105,3.448025235887442
+-3.052540789244851,-3.368538874658002,-2.076646325437413
+0.4832288185628828,0.6237050700811455,1.558707817707959
+-2.248508782339848,-1.482160721883134,-1.58464357170877
+-2.517929282061579,-2.775270769084503,-3.27602157472457
+1.634568674859941,-2.710142286577208,2.243207255140116
+3.008986605544752,3.257049395779147,-0.03501022355988503
+2.515025114040092,1.917473895667771,-1.225134426980443
+0.5115216512046832,0.7179358363111685,-2.487666607985341
+3.044260743910304,-2.862594876583988,-1.46201137215194
+-0.2383931340630649,1.786059765469997,1.351269014921268
+-2.635369291339539,3.090088683267917,-0.5134386557105834

examples/Y0.csv

@@ -0,0 +1,80 @@
+14.69972565901391,-1.141468833465243,0.4397628058490257
+12.00467177123156,0.3450566535153889,-0.01473908746398855
+8.47961540962646,-3.287777076316099,1.308826713376965
+8.398538411157281,-3.01468837280452,1.849933218363242
+14.88049054873884,-3.737191079416643,-3.111293341819557
+11.74148973091656,-0.3484659609597651,0.8657812321388532
+11.04110550473725,0.7311205721480922,-1.484468077090113
+14.83574476650973,2.466228041426272,0.0516046693835106
+8.408072302382465,1.812773714425157,2.018471489371509
+14.05371411510352,-2.94161729377125,1.710374703538183
+14.68599916016895,2.606027045615759,-3.761562476391337
+11.94854696068537,-3.567674502239928,-3.231224919011964
+11.95997114829134,-0.6908649095010624,2.738962946269394
+9.812413013389474,-2.776469076240756,2.139536346417309
+9.174581477712177,-0.2409673622942501,-3.671616348336108
+13.38797467178571,-2.038087522910148,-3.371974526774131
+10.3341815471754,-0.697658984727064,1.169218562354129
+11.50622066461022,-2.624588875782207,-0.3508002284925729
+11.71935018240397,-0.1146741774843543,0.5592192460310068
+11.20640672181586,-1.793663187647437,0.4458471522208551
+12.42929792160473,-2.428047879697353,-0.3990291598062505
+13.02637016363813,1.015763663159283,-2.239240292858172
+8.562214531355789,3.179329412369906,-1.745619707387359
+9.896807560889087,-1.400019856040225,-1.739532832091071
+10.59304383943116,-2.802233958335766,3.433480512118521
+8.909511205874308,-2.672426019910846,2.044803212501556
+11.16766198802251,2.983659737938954,-1.171675823994299
+15.4327835253383,-2.275324074587767,-1.709757976852235
+14.17661591629654,0.09515702147714311,0.1163884174999934
+9.575837719548321,-0.3708408807348693,-2.236078964630921
+13.23050940395724,2.945290636556644,-3.653251006453718
+8.483041255398991,-2.157211974614906,-1.539945203228259
+10.48758162444119,2.088218212612035,-2.659188360759875
+12.96170600713591,-1.588018898419284,1.381801190081896
+12.05771261589465,1.317536569107611,-2.539715678653259
+15.58041947902209,-3.652895170068064,-0.1469802529782429
+12.41539485418864,-1.692018838059387,-2.445186058249809
+15.3585766656192,1.160183286340198,-1.458187941565711
+13.65468761763591,-0.38438861691135,1.578091871295611
+8.565667906981417,2.493226372866149,-2.160723270542434
+14.81684801065015,0.2166042521152436,-2.397459318886075
+10.58466120033887,-3.44907222239941,-0.008520673548414237
+10.79368870007821,-3.432172564606917,2.857572248764802
+10.4554581651538,1.443655613314474,0.996416483624766
+9.644239364278235,0.4673524547632306,2.583488909748998
+9.085950867078303,-0.1989639465749384,2.904416148321347
+13.14713045409713,-0.3133568651937866,1.588952654970732
+13.02870753914112,-3.916736332368862,0.09204855318969685
+14.60760455878589,3.313021655317749,-2.409720202437378
+10.88460409535118,2.420319189083149,-0.808699854228232
+14.30183423904584,-0.8616365729135129,1.863785632368431
+14.42840386421026,1.3128573195545,-3.335021363504484
+8.079261505296593,0.977957113296207,-3.59415742343417
+11.03662660434383,-1.559305186719115,-2.878542825148496
+10.15482365241914,2.403240366901536,0.2239871454969888
+11.98515534577828,3.764944543948487,-2.975058395592636
+14.93999661956639,-3.766644636495617,2.659933026675608
+8.718543216304262,1.471706294886816,-2.335908617928363
+10.79002876620167,3.201997039005546,2.157719526149253
+15.60392547020433,3.554480623796912,2.001995101335606
+8.688094639785552,-1.437166289112112,3.465469880899794
+12.25561819305717,0.6080516384134951,-2.604455837975218
+11.21792210763533,1.280354123835233,-2.526414253090011
+13.02566151092011,-3.837472176096592,-2.644067464772272
+12.82676622818703,-3.436770286186816,-0.4911104604117247
+9.77273992491547,-1.71817103485717,2.841030122100281
+14.12663905112325,-3.220269124523254,0.8474705229515056
+15.93003911380492,2.669892288699828,2.257057535583721
+13.49800511699696,-3.77942975727085,-0.4835008098131581
+8.117019220233026,3.929185790473598,2.830587153373334
+15.03621175621608,-0.77279672588885,2.303531983934888
+11.20075535563193,-3.371797876578698,-1.794896345037752
+9.079383956279774,-3.278999616500193,-0.944884026835445
+11.72809049063392,-1.630581815041499,-1.436742007842677
+11.25147817317208,-3.408120683233418,-0.147489366208716
+11.70242916441565,-0.8750391154790824,1.22973168598537
+13.42210547370169,2.751327702937805,2.469533643655149
+13.09715862046082,-3.366097859028909,2.210820684261365
+8.587144239372094,2.445616685742658,-0.4268829236380336
+10.30285545022036,0.3764137264011644,3.858750112222492

examples/eigs.rs

@@ -0,0 +1,67 @@
+//! Example of GMRES iterations
+use num::Complex;
+use rlst::operator::operations::eigenvalues::eigs::Which;
+use rlst::operator::Operator;
+use rlst::{c64, rlst_dynamic_array2, Eigs};
+
+pub fn main() {
+    let mut mat = rlst_dynamic_array2!(f64, [5, 5]);
+    mat.r_mut()[[0, 0]] = 4.0;
+    mat.r_mut()[[0, 1]] = 1.0;
+
+    mat.r_mut()[[1, 0]] = 1.0;
+    mat.r_mut()[[1, 1]] = 4.0;
+    mat.r_mut()[[1, 2]] = 1.0;
+
+    mat.r_mut()[[2, 1]] = 1.0;
+    mat.r_mut()[[2, 2]] = 4.0;
+    mat.r_mut()[[2, 3]] = 1.0;
+
+    mat.r_mut()[[3, 2]] = 1.0;
+    mat.r_mut()[[3, 3]] = 4.0;
+    mat.r_mut()[[3, 4]] = 1.0;
+
+    mat.r_mut()[[4, 3]] = 1.0;
+    mat.r_mut()[[4, 4]] = 4.0;
+
+    // We can now wrap the matrix into an operator.
+    let op = Operator::from(mat);
+
+    let k = 3;
+    let mut eigs = Eigs::new(op, 1e-10, None, None, Some(Which::LM));
+
+    let (res1, _res2) = eigs.run(None, k, None, false);
+
+    println!("Eigenvalues: {:?}", res1);
+
+    // We can do the same for a complex case
+    let mut mat = rlst_dynamic_array2!(c64, [4, 4]);
+    mat.r_mut()[[0, 0]] = Complex::new(1.0, 2.0);
+    mat.r_mut()[[0, 1]] = Complex::new(2.0, -1.0);
+    mat.r_mut()[[0, 2]] = Complex::new(0.0, 0.0);
+    mat.r_mut()[[0, 3]] = Complex::new(1.0, 0.0);
+
+    mat.r_mut()[[1, 0]] = Complex::new(0.0, 0.0);
+    mat.r_mut()[[1, 1]] = Complex::new(3.0, 0.0);
+    mat.r_mut()[[1, 2]] = Complex::new(1.0, 1.0);
+    mat.r_mut()[[1, 3]] = Complex::new(0.0, 0.0);
+
+    mat.r_mut()[[2, 0]] = Complex::new(4.0, 0.0);
+    mat.r_mut()[[2, 1]] = Complex::new(0.0, 0.0);
+    mat.r_mut()[[2, 2]] = Complex::new(2.0, -1.0);
+    mat.r_mut()[[2, 3]] = Complex::new(1.0, 0.0);
+
+    mat.r_mut()[[3, 0]] = Complex::new(0.0, 1.0);
+    mat.r_mut()[[3, 1]] = Complex::new(0.0, 0.0);
+    mat.r_mut()[[3, 2]] = Complex::new(0.0, 0.0);
+    mat.r_mut()[[3, 3]] = Complex::new(1.0, 3.0);
+
+    let op = Operator::from(mat);
+
+    let k = 2;
+    let mut eigs = Eigs::new(op, 1e-10, None, None, Some(Which::LM));
+
+    let (res1, _res2) = eigs.run(None, k, None, false);
+
+    println!("Eigenvalues: {:?}", res1);
+}

examples/gmres.rs

@@ -0,0 +1,100 @@
+//! Example of GMRES iterations
+use num::Complex;
+use rlst::operator::operations::gmres::IdOperator;
+use rlst::operator::Operator;
+use rlst::{c64, rlst_dynamic_array2, zero_element, GmresIteration, OperatorBase, RawAccessMut};
+
+pub fn main() {
+    let mut mat = rlst_dynamic_array2!(f64, [5, 5]);
+    mat.r_mut()[[0, 0]] = 4.0;
+    mat.r_mut()[[0, 1]] = 1.0;
+
+    mat.r_mut()[[1, 0]] = 1.0;
+    mat.r_mut()[[1, 1]] = 4.0;
+    mat.r_mut()[[1, 2]] = 1.0;
+
+    mat.r_mut()[[2, 1]] = 1.0;
+    mat.r_mut()[[2, 2]] = 4.0;
+    mat.r_mut()[[2, 3]] = 1.0;
+
+    mat.r_mut()[[3, 2]] = 1.0;
+    mat.r_mut()[[3, 3]] = 4.0;
+    mat.r_mut()[[3, 4]] = 1.0;
+
+    mat.r_mut()[[4, 3]] = 1.0;
+    mat.r_mut()[[4, 4]] = 4.0;
+
+    // Wrap matrix into operator
+    let op = Operator::from(mat);
+
+    let id_operator = IdOperator::new(op.r().domain(), op.r().range());
+    // Right-hand side
+    let mut rhs = zero_element(op.range());
+    rhs.view_mut().data_mut()[0] = 1.0;
+    rhs.view_mut().data_mut()[1] = 2.0;
+    rhs.view_mut().data_mut()[2] = 3.0;
+    rhs.view_mut().data_mut()[3] = 4.0;
+    rhs.view_mut().data_mut()[4] = 5.0;
+
+    let mut residuals = Vec::<f64>::new();
+    let tol = 1e-3;
+
+    // GMRES solve
+    let gmres = (GmresIteration::new(op.r(), rhs.r(), 5))
+        .set_callable(|_, res| {
+            residuals.push(res);
+        })
+        .set_tol(tol)
+        .set_restart(5)
+        .set_preconditioner(id_operator.r());
+
+    let (_sol, _res) = gmres.run();
+
+    println!("Residuals: {:?}", residuals);
+    // As an example, we set up a complex matrix and
+    // hen wrap it in an operator to run GMRES.
+    let mut mat = rlst_dynamic_array2!(c64, [5, 5]);
+    mat.r_mut()[[0, 0]] = Complex::new(4.0, 1.0);
+    mat.r_mut()[[0, 1]] = Complex::new(1.0, -1.0);
+
+    mat.r_mut()[[1, 0]] = Complex::new(1.0, 1.0);
+    mat.r_mut()[[1, 1]] = Complex::new(4.0, 2.0);
+    mat.r_mut()[[1, 2]] = Complex::new(1.0, -1.0);
+
+    mat.r_mut()[[2, 1]] = Complex::new(1.0, 1.0);
+    mat.r_mut()[[2, 2]] = Complex::new(4.0, 3.0);
+    mat.r_mut()[[2, 3]] = Complex::new(1.0, -1.0);
+
+    mat.r_mut()[[3, 2]] = Complex::new(1.0, 1.0);
+    mat.r_mut()[[3, 3]] = Complex::new(4.0, 4.0);
+    mat.r_mut()[[3, 4]] = Complex::new(1.0, -1.0);
+
+    mat.r_mut()[[4, 3]] = Complex::new(1.0, 1.0);
+    mat.r_mut()[[4, 4]] = Complex::new(4.0, 5.0);
+
+    // We can now wrap the matrix into an operator.
+    let op = Operator::from(mat);
+    let id_operator = IdOperator::new(op.r().domain(), op.r().range());
+    // Let's create a right-hand side.
+    let mut rhs = zero_element(op.range());
+    rhs.view_mut().data_mut()[0] = Complex::new(1.0, 2.0);
+    rhs.view_mut().data_mut()[1] = Complex::new(2.0, 3.0);
+    rhs.view_mut().data_mut()[2] = Complex::new(3.0, 4.0);
+    rhs.view_mut().data_mut()[3] = Complex::new(4.0, 5.0);
+    rhs.view_mut().data_mut()[4] = Complex::new(5.0, 6.0);
+
+    let mut residuals = Vec::<f64>::new();
+    let tol = 1E-3;
+
+    // We can now run the GMRES iteration. We must specify the dimension in this case
+    let gmres = (GmresIteration::new(op.r(), rhs.r(), 5))
+        .set_callable(|_, res| {
+            residuals.push(res);
+        })
+        .set_tol(tol)
+        .set_restart(5)
+        .set_preconditioner(id_operator.r());
+    let (_sol, _res) = gmres.run();
+
+    println!("Residuals: {:?}", residuals)
+}