31:T825,ja = new Int32Array([5, 6, 11, 13, 21, 24, 26, 28, 30, 32, 34, 36, 37, 1, 2, 3, 6, 7, 16, 20, 22, 26, 29, 31, 38, 39, 0, 1, 5, 8, 13, 16, 22, 24, 25, 30, 31, 37, 39, 0, 2, 12, 15, 17, 21, 23, 25, 31, 32, 38, 0, 1, 8, 9, 10, 14, 16, 25, 26, 27, 30, 32, 36, 39, 0, 2, 6, 10, 12, 14, 15, 16, 24, 28, 34, 36, 37, 5, 8, 9, 10, 11, 14, 15, 18, 21, 25, 30, 34, 36, 0, 2, 3, 6, 9, 10, 11, 12, 13, 14, 15, 16, 19, 22, 24, 27, 28, 31, 34, 35, 37, 38, 2, 4, 8, 9, 11, 15, 17, 18, 19, 20, 23, 26, 28, 30, 39, 0, 7, 8, 11, 12, 13, 17, 20, 23, 24, 25, 26, 27, 29, 30, 34, 35, 39, 0, 4, 5, 6, 7, 9, 13, 18, 19, 20, 21, 29, 32, 36, 37, 1, 4, 6, 7, 8, 9, 10, 12, 14, 19, 23, 27, 29, 30, 31, 33, 35, 36, 39, 4, 5, 8, 9, 12, 14, 18, 20, 22, 23, 26, 28, 29, 33, 38, 2, 4, 5, 6, 10, 15, 18, 19, 25, 26, 29, 32, 33, 34, 38, 0, 5, 6, 7, 12, 17, 19, 20, 32, 34, 35, 36, 37, 2, 5, 10, 11, 12, 13, 14, 15, 17, 19, 21, 24, 28, 31, 33, 34, 37, 4, 5, 6, 7, 8, 9, 14, 16, 18, 19, 22, 24, 25, 27, 29, 36, 37, 38, 6, 7, 8, 11, 14, 15, 20, 21, 26, 27, 35, 37, 39, 8, 10, 12, 13, 17, 18, 19, 21, 26, 28, 34, 39, 1, 3, 4, 5, 9, 16, 19, 27, 28, 30, 31, 32, 36, 37, 0, 5, 9, 11, 13, 16, 21, 25, 28, 29, 32, 33, 34, 37, 39, 1, 11, 13, 17, 21, 27, 28, 30, 34, 35, 37, 6, 12, 13, 15, 16, 19, 20, 21, 22, 23, 27, 38, 39, 3, 6, 7, 8, 10, 18, 26, 27, 28, 31, 34, 0, 2, 4, 7, 8, 13, 16, 17, 18, 20, 24, 29, 31, 33, 1, 5, 7, 10, 11, 15, 18, 25, 36, 3, 6, 9, 11, 15, 16, 18, 19, 21, 23, 28, 29, 34, 35, 37, 2, 4, 18, 22, 23, 24, 30, 32, 33, 36, 37, 39, 1, 3, 5, 10, 14, 18, 22, 25, 28, 31, 38, 2, 3, 4, 8, 9, 13, 15, 16, 18, 19, 24, 25, 28, 30, 33, 35, 36, 38, 0, 4, 5, 11, 12, 17, 20, 22, 24, 26, 27, 34, 35, 0, 1, 3, 13, 17, 21, 22, 26, 27, 29, 35, 0, 2, 3, 8, 11, 14, 27, 30, 31, 33, 38, 1, 2, 10, 14, 16, 19, 22, 24, 28, 29, 31, 33, 38, 2, 4, 6, 10, 13, 14, 16, 22, 23, 29, 30, 32, 35, 38, 4, 7, 17, 18, 21, 22, 32, 33, 36, 38, 39, 3, 9, 12, 15, 21, 25, 27, 32, 33, 34, 35, 1, 2, 6, 7, 9, 10, 14, 16, 20, 31, 32, 35, 36, 39, 1, 2, 8, 9, 13, 15, 17, 20, 22, 23, 24, 31, 38, 39, 1, 4, 7, 11, 15, 23, 24, 25, 26, 31, 33, 35, 38]);32:T26ba,a = new Float64Array([0.66409721489255, 0.84198102161381, 0.086985417216544, 0.73676829493454, 0.8640918027908, 0.38778369007063, 0.16104673881133, 0.23472161089756, 0.78371407454075, 0.42437513518351, 0.024136640142713, 0.097865603444104, 0.84387106115179, 0.40016761161395, 0.69923549084889, 0.66702987796954, 0.85354897931849, 0.53823151650756, 0.65741016280717, 0.2866589153589, 0.96450714811893, 0.48811498865863, 0.36257565597192, 0.85061888203519, 0.013656753587377, 0.52713099705387, 0.065797721066418, 0.17739803026309, 0.91163121345995, 0.088972918730682, 0.9539713808121, 0.32966474784988, 0.98276154183911, 0.46334126194163, 0.66517087522204, 0.15530571767842, 0.59740136684729, 0.94834261618012, 0.24604026984705, 0.74229067738275, 0.45494995613347, 0.27949460841692, 0.38072374108281, 0.59645310770555, 0.73620055184523, 0.22060646918631, 0.6646953083876, 0.27446741437282, 0.63267413928764, 0.85027011476935, 0.076471245417591, 0.69235237720066, 0.63454246457412, 0.10859496384328, 0.63197515561803, 0.7907377494456, 0.81536307596386, 0.81703976579804, 0.0042168763485816, 0.49465305474338, 0.43954942349324, 0.26972704905538, 0.75249431037926, 0.048900618240657, 0.97404789644016, 0.045381852912429, 0.83543993059333, 0.16537250027311, 0.41784518836897, 0.10470877546105, 0.68321012877077, 0.80950425370108, 0.5383067436229, 0.52140686359322, 0.51983720740296, 0.40603951337097, 0.826644857799, 0.90021519777375, 0.25824479025707, 0.60824775305029, 0.29895941740785, 0.85713173582085, 0.086828700772873, 0.58318777828626, 0.69904044768728, 0.95569524166905, 0.2209895561547, 0.65667160631003, 0.18771165245572, 0.56381605219274, 0.64488580992673, 0.1191987163011, 0.96312462536764, 0.34366194500758, 0.70841572839227, 0.092106863899206, 0.89669200028139, 0.23998477227985, 0.18544674999334, 0.46504025555451, 0.6157540216184, 0.59274066639726, 0.012456462724347, 0.2183824406091, 0.41556357146034, 0.32293527122724, 0.83634807161817, 0.839323311038, 0.1180727817668, 0.067804348220958, 0.0168316830028, 0.98710941569279, 0.067667180703798, 0.99725481355435, 0.96254051195576, 0.16354285607279, 0.95137052608252, 0.34901871129359, 0.34084752404171, 0.43883430466002, 0.96618846802329, 0.26318011352009, 0.22334194938808, 0.30256319758602, 0.40367053281687, 0.14676253131906, 0.063912935584743, 0.10605741576574, 0.25270096922885, 0.50208168546766, 0.58305751233504, 0.45182946438521, 0.4616837308098, 0.090379560874021, 0.32335861135431, 0.57247322591603, 0.52984879935619, 0.42524448755441, 0.88823030371602, 0.46542373414404, 0.41597074569015, 0.083536626344331, 0.26222169225208, 0.90931407823661, 0.70966936494674, 0.87580041535003, 0.66117594142499, 0.47075983158814, 0.87090980162421, 0.36050090443366, 0.86851316172095, 0.6799574078433, 0.7745788059079, 0.26199168025609, 0.36368794337087, 0.23506123211005, 0.9521097256579, 0.4349337995215, 0.86129106621318, 0.47955187665278, 0.98045169374927, 0.69774883971445, 0.46728924543936, 0.33193904689138, 0.087280241347514, 0.32964839382546, 0.45209510226366, 0.19578087245849, 0.054554409372878, 0.9325557583629, 0.28983588530209, 0.16614358879912, 0.12937159795751, 0.79695714814447, 0.63796531578431, 0.69116359655334, 0.6566041087064, 0.3164737016505, 0.99781568860533, 0.95958142911996, 0.67137671060459, 0.45421183083868, 0.010211845864641, 0.3420369244842, 0.29203939591164, 0.76199786121119, 0.13811758167023, 0.028742313398394, 0.30381065993747, 0.033628078658892, 0.037584267574169, 0.49243469605801, 0.69312154114857, 0.16271650798745, 0.98056168760199, 0.38932682452226, 0.53628787656142, 0.65222502809587, 0.56253893233954, 0.37130878091385, 0.79117077160216, 0.62218569248085, 0.77720587410834, 0.76728066837754, 0.34080476283133, 0.38598243165109, 0.20025661550474, 0.88695712149467, 0.99892438249612, 0.55910784916911, 0.89780516917715, 0.83775788165525, 0.47000354922842, 0.31529232781161, 0.084250014314544, 0.70394091014934, 0.13394123927408, 0.82585938592714, 0.4215777904827, 0.80308818156043, 0.0012260861700522, 0.1680285437815, 0.7902342238418, 0.55499938901281, 0.31224083169933, 0.92806459028649, 0.2420692142295, 0.39609093563449, 0.21914706296248, 0.5283689221965, 0.67927487505566, 0.2301632246143, 0.72912379807286, 0.49473904329107, 0.41151909223363, 0.45780567520196, 0.34025149156351, 0.3164872640355, 0.93844740229587, 0.3651033003652, 0.6165725503194, 0.84540067978455, 0.37591202388327, 0.78814887524962, 0.3644102720378, 0.73768145718504, 0.064922351420355, 0.96139853259614, 0.17633116439745, 0.41544446415056, 0.47244308910912, 0.066939796352265, 0.21064433372144, 0.948309211502, 0.40649708798458, 0.52105481294964, 0.94304930835173, 0.86209596826793, 0.62740343605513, 0.53512420716468, 0.87615022010922, 0.76083926659117, 0.57008802498229, 0.51158007071893, 0.96231187878284, 0.29276684685273, 0.34704392559223, 0.1821283447473, 0.27888971021347, 0.93254981512788, 0.63901270070999, 0.83034742056874, 0.65817396094006, 0.37613910314447, 0.92425559364457, 0.90825014184613, 0.24631454946767, 0.77816218593072, 0.0060010682819416, 0.7341695026188, 0.5955489471534, 0.49783113528873, 0.72781118551633, 0.66225712590956, 0.7389398006438, 0.47505937538811, 0.28376496037644, 0.36458299512257, 0.49072650656604, 0.51917632972783, 0.71484031608088, 0.10374841098848, 0.32517009988668, 0.71911738473881, 0.56227632358776, 0.26972431841759, 0.17937273773336, 0.36611479444714, 0.67785077853028, 0.74924887518829, 0.87669388059373, 0.53247514810994, 0.72014530595399, 0.24130796419518, 0.89353328239803, 0.91474564928317, 0.8568825315949, 0.78219005827894, 0.079410784449154, 0.34469563995707, 0.74882430245579, 0.49253279412749, 0.73667793149905, 0.40242958692947, 0.89357904712371, 0.76445568388535, 0.49091491684826, 0.29588948297123, 0.82801906151139, 0.33152289704025, 0.28011863831436, 0.42024107855756, 0.86541421518913, 0.32689955752664, 0.43324178011773, 0.58557211448698, 0.77960841719974, 0.85975404403161, 0.57419783183103, 0.07184490983926, 0.2184185628865, 0.45143160384681, 0.60698629152355, 0.39187620225916, 0.40578644136236, 0.97716984384561, 0.16895345746025, 0.77094133932653, 0.32650160990958, 0.83766613753404, 0.95224665522214, 0.90949811083707, 0.0092934235042396, 0.59570208033347, 0.30175896049559, 0.26463790855586, 0.83370270246347, 0.55775214152306, 0.82235520417912, 0.7795484954396, 0.68825480699924, 0.79787023868313, 0.76447208307892, 0.866067035527, 0.59562939898839, 0.5910492556128, 0.63083996280601, 0.70544138630174, 0.99306430993279, 0.75679666630775, 0.17000235578511, 0.66414927116788, 0.30998595725279, 0.60282881399748, 0.14451311721677, 0.63912770274986, 0.76688914222964, 0.92015726441525, 0.28946118535915, 0.36297261079912, 0.33081542855632, 0.077584788705029, 0.82335714196011, 0.65077656817193, 0.39485370991512, 0.83880657974575, 0.84802095631511, 0.60272482950367, 0.85445700811895, 0.89812707756559, 0.61188216209965, 0.52371965419674, 0.50468107243287, 0.88919592317622, 0.47637113345618, 0.42407930522416, 0.92980666874433, 0.29910925044637, 0.47354941697491, 0.50776368962031, 0.44764549818246, 0.60804483136537, 0.3715119116807, 0.86661882832023, 0.38576697808959, 0.13372335822029, 0.73497184819308, 0.24546756280841, 0.37115529988481, 0.52904775530521, 0.070124424095324, 0.50181817752394, 0.67260104262857, 0.83276055745443, 0.14863509458892, 0.43772465197264, 0.99723941041028, 0.53125292925688, 0.18484625508303, 0.96057264644679, 0.36008066980172, 0.75400907627959, 0.02332577110423, 0.76534378936763, 0.24261500418308, 0.51724102325609, 0.38662863727036, 0.92349223137065, 0.66844386126308, 0.29406603439435, 0.88001435104758, 0.24085031135047, 0.31939553018631, 0.41957975943553, 0.46050797843398, 0.88158148521631, 0.71018992537176, 0.9050517156278, 0.32363275313919, 0.91786842789402, 0.7728007057555, 0.96934607530448, 0.66778048904044, 0.37203299131805, 0.61298084566974, 0.70801004241594, 0.64868363721701, 0.10132475341732, 0.73439186379984, 0.48581840772453, 0.15738899873448, 0.032589597177035, 0.89255448053244, 0.39578449465138, 0.32574296525016, 0.51726393891371, 0.87528906989623, 0.30485229906852, 0.81157375025171, 0.16405088462124, 0.86235585196985, 0.064627548709804, 0.53317611409965, 0.2220269484548, 0.88269043708345, 0.082297921219048, 0.96366148486904, 0.83976589461778, 0.087819044519132, 0.97993288095106, 0.80399745833315, 0.36103780770723, 0.85339976141853, 0.17507189986067, 0.7354778357481, 0.22558387053459, 0.49166348878837, 0.74962741497421, 0.15058505309307, 0.2824946624611, 0.386168069386, 0.89072686289005, 0.064191195212394, 0.25423533434711, 0.51557263430002, 0.61019803425772, 0.99825095245533, 0.84769435173259, 0.75066340330553, 0.80584430173312, 0.14365524665623, 0.69047203552466, 0.4211538962187, 0.0023073316562489, 0.85397449035848, 0.49806651635937, 0.85461601654748, 0.45061970290291, 0.34775738248032, 0.97727156382858, 0.14744104777809, 0.31511626919504, 0.7035026432497, 0.59201451604814, 0.19412905545632, 0.20894501274868, 0.69308936348795, 0.87096057686534, 0.12921391526666, 0.1444018656129, 0.093096737327565, 0.022832828584515, 0.85579704952231, 0.11565840948171, 0.15990576202046, 0.75305123988215, 0.77118766390308, 0.81809826140203, 0.35093810518782, 0.77150895389333, 0.90893232212818, 0.61360928584524, 0.49531520646779, 0.35935842728212, 0.19411258687922, 0.94410850337898, 0.80554843172689, 0.16489376088832, 0.63628597354343, 0.74645545601214, 0.12300008773943, 0.3709410589053, 0.34068128109941, 0.67384117593702, 0.87835343781782, 0.52485860210138, 0.34031006989084, 0.031038194443583, 0.78729296931405, 0.30877335896193, 0.24698937556101, 0.3872920877241, 0.33903670792423, 0.0057043447185794, 0.0087801571976301, 0.67155921211073, 0.21207364844721, 0.65415653989378, 0.66649673258257, 0.57024202196404, 0.18992019686379, 0.49792923987747, 0.52621778264931, 0.75353723799509, 0.66908486218615, 0.4471997508999, 0.013198933570273, 0.24060259211837, 0.85191674011383, 0.85572038956719]);33:T5f8,class Permutation { static Apply(p, b, x, n) { if (p == null) { x.set(b); } else { for (let k = 0; k < n; k++) { x[k] = b[p[k]]; } } } static ApplyInverse(p, b, x, n) { if (p == null) { x.set(b); } else { for (let k = 0; k < n; k++) { x[p[k]] = b[k]; } } } static Create(n, seed = 0) { let i, j, tmp; let p = new Int32Array(n); if (seed == 0) { // Identity for (i = 0; i < n; i++) { p[i] = i; } } else { // Inverse for (i = 0; i < n; i++) { p[i] = n - i - 1; } } if (seed > 0) { // Randomize permutation. for (i = 0; i < n; i++) { let rand = Math.random() * 2147483647; // j = rand integer in range k to n-1 j = i + (rand.Next() % (n - i)); // swap p[k] and p[j] tmp = p[j]; p[j] = p[i]; p[i] = tmp; } } return p; } static Invert(p) { let k, n = p.length; let pinv = new Int32Array(n); for (k = 0; k < n; k++) { // Invert the permutation. pinv[p[k]] = k; } return pinv; } static IsValid(p) { let length = p.length; let check = new Int8Array(length); for (let i = 0; i < length; i++) { if (p[i] >= length || p[i] < 0) { return false; } check[p[i]] = true; } for (let i = 0; i < length; i++) { if (check[i] == false) { return false; } } return true; } }34:T1801,class SparseLU { static Create(A, p = null, tol = 1.0, progress = null) { should.exist(A); let n = A.ColumnCount; if (p == null) { p = Permutation.Create(n); } // FIXME: Check.SquareMatrix(A, "A"); // FIXME: Check.Permutation(p, n, "p"); // FIXME: Check.NotNaN(tol, "tol"); // Ensure tol is in range. tol = Math.min(Math.max(tol, 0.0), 1.0); let C = new SparseLU(n); // Ordering and symbolic analysis C.SymbolicAnalysis(A, p); // Numeric Cholesky factorization C.Factorize(A, tol, progress); return C; } constructor(n) { this.n = n; this.temp = new Float64Array(n); // workspace this.S = null; // SymbolicFactorization this.L = null; // CompressedColumnStorage this.U = null; // CompressedColumnStorage this.pinv = new Int32Array(0); // partial pivoting } get NonZerosCount() { return (this.L.NonZerosCount + this.U.NonZerosCount - this.n); } Solve(input, result) { if (input == null) throw "input"; if (result == null) throw "result"; let x = this.temp; Permutation.ApplyInverse(this.pinv, input, x, this.n); // x = b(p) SolverHelper.SolveLower(this.L, x); // x = L\x. SolverHelper.SolveUpper(this.U, x); // x = U\x. Permutation.ApplyInverse(this.S.q, x, result, this.n); // b(q) = x } SolveTranspose(input, result) { if (input == null) throw "input"; if (result == null) throw "result"; let x = this.temp; Permutation.Apply(this.S.q, input, x, this.n); // x = Q'*b SolverHelper.SolveUpperTranspose(this.U, x); // x = U'\x. SolverHelper.SolveLowerTranspose(this.L, x); // x = L'\x. Permutation.Apply(this.pinv, x, result, this.n); // b = P'*x } Factorize(A, tol, progress) { let q = this.S.q; let i; let lnz = this.S.lnz; let unz = this.S.unz; this.L = CompressedColumnStorage.Create(this.n, this.n, lnz); this.U = CompressedColumnStorage.Create(this.n, this.n, unz); this.pinv = new Float32Array(this.n); // Workspace let x = this.temp; let xi = new Float32Array(2 * this.n); for (i = 0; i < this.n; i++) { // No rows pivotal yet. this.pinv[i] = -1; } lnz = unz = 0; let ipiv, top, p, col; let pivot; let a, t; let li, ui; let lp = this.L.ColumnPointers; let up = this.U.ColumnPointers; let lx, ux; let current = 0.0; let step = this.n / 100.0; // Now compute L(:,k) and U(:,k) for (let k = 0; k < this.n; k++) { // Progress reporting. if (k >= current) { current += step; if (progress != null) { progress(k / this.n); } } // Triangular solve lp[k] = lnz; // L(:,k) starts here up[k] = unz; // U(:,k) starts here if (lnz + this.n > this.L.Values.length) { this.L.Resize(2 * this.L.Values.length + this.n); } if (unz + this.n > this.U.Values.length) { this.U.Resize(2 * this.U.Values.length + this.n); } li = this.L.RowIndices; ui = this.U.RowIndices; lx = this.L.Values; ux = this.U.Values; col = q != null ? (q[k]) : k; top = this.SolveSp(this.L, A, col, xi, x, this.pinv, true); // x = L\A(:,col) // Find pivot ipiv = -1; a = -1; for (p = top; p < this.n; p++) { i = xi[p]; // x(i) is nonzero if (this.pinv[i] < 0) /* Row i is not yet pivotal */ { if ((t = Math.abs(x[i])) > a) { a = t; // Largest pivot candidate so far ipiv = i; } } else /* x(i) is the entry U(pinv[i],k) */ { ui[unz] = this.pinv[i]; ux[unz++] = x[i]; } } if (ipiv == -1 || a <= 0.0) { throw "No pivot element found."; } if (this.pinv[col] < 0 && Math.abs(x[col]) >= a * tol) { ipiv = col; } // Divide by pivot pivot = x[ipiv]; // the chosen pivot ui[unz] = k; // last entry in U(:,k) is U(k,k) ux[unz++] = pivot; this.pinv[ipiv] = k; // ipiv is the kth pivot row li[lnz] = ipiv; // first entry in L(:,k) is L(k,k) = 1 lx[lnz++] = 1.0; for (p = top; p < this.n; p++) /* L(k+1:n,k) = x / pivot */ { i = xi[p]; if (this.pinv[i] < 0) /* x(i) is an entry in L(:,k) */ { li[lnz] = i; // save unpermuted row in L lx[lnz++] = x[i] / pivot; // scale pivot column } x[i] = 0.0; // x [0..n-1] = 0 for next k } } // Finalize L and U lp[this.n] = lnz; up[this.n] = unz; li = this.L.RowIndices; // fix row indices of L for final pinv for (p = 0; p < lnz; p++) { li[p] = this.pinv[li[p]]; } // Remove extra space from L and U this.L.Resize(0); this.U.Resize(0); } SymbolicAnalysis(A, p) { let sym = this.S = new SymbolicFactorization(); // Fill-reducing ordering sym.q = p; // Guess nnz(L) and nnz(U) sym.unz = sym.lnz = 4 * (A.ColumnPointers[this.n]) + this.n; } SolveSp(G, B, k, xi, x, pinv, lo) { if (xi == null || x == null) return -1; let gp = G.ColumnPointers; let gi = G.RowIndices; let gx = G.Values; let bp = B.ColumnPointers; let bi = B.RowIndices; let bx = B.Values; let n = G.ColumnCount; // xi[top..n-1]=Reach(B(:,k)) let top = GraphHelper.Reach(gp, gi, bp, bi, n, k, xi, pinv); let j, J, p, q, px; for (p = top; p < n; p++) { x[xi[p]] = 0; // clear x } for (p = bp[k]; p < bp[k + 1]; p++) { x[bi[p]] = bx[p]; // scatter B } for (px = top; px < n; px++) { j = xi[px]; // x(j) is nonzero J = pinv != null ? (pinv[j]) : j; // j maps to col J of G if (J < 0) continue; // column J is empty x[j] /= gx[lo ? (gp[J]) : (gp[J + 1] - 1)]; // x(j) /= G(j,j) p = lo ? (gp[J] + 1) : (gp[J]); // lo: L(j,j) 1st entry q = lo ? (gp[J + 1]) : (gp[J + 1] - 1); // up: U(j,j) last entry for (; p < q; p++) { x[gi[p]] -= gx[p] * x[j]; // x(i) -= G(i,j) * x(j) } } // Return top of stack. return top; } }35:T719,class CompressedColumnStorage { get NonZerosCount() { return this.ColumnPointers[this.columns]; } constructor(rowCount, columnCount, values, rowIndices, columnPointers) { this.rows = rowCount this.columns = columnCount; if (columnPointers.length !== columnCount + 1) { throw "columnPointers.Length must equal (columnCount + 1)"; } if (values.length !== rowIndices.length) { throw "rowIndices.Length must equal values.Length"; } this.ColumnPointers = columnPointers; this.RowIndices = rowIndices; this.Values = values; } static Create(rowCount, columnCount, valueCount = 0) { return new CompressedColumnStorage(rowCount, columnCount, new Float64Array(valueCount), new Int32Array(valueCount), new Int32Array(columnCount + 1)); } Resize(size) { if (size <= 0) { size = this.ColumnPointers[this.columns]; } // TODO: check available memory // and throw OutOfMemoryException before trying to resize. let RowIndices = new Int32Array(size); for (let i = 0; i < Math.min(RowIndices.length, this.RowIndices.length); i++) { RowIndices[i] = this.RowIndices[i]; } this.RowIndices = RowIndices; let Values = new Float64Array(size); for (let i = 0; i < Math.min(Values.length, this.Values.length); i++) { Values[i] = this.Values[i]; } this.Values = Values; return true; } get RowCount() { return this.rows; } get ColumnCount() { return this.columns; } *EnumerateIndexed() { let n = this.ColumnCount; let ax = this.Values; let ap = this.ColumnPointers; let ai = this.RowIndices; for (let i = 0; i < n; i++) { let end = ap[i + 1]; for (let j = ap[i]; j < end; j++) { yield [ai[j], i, ax[j]]; } } } }36:T210f,class GraphHelper { static DepthFirstSearch(j, Gp, Gi, top, xi, pstack, offset, pinv) { let i, p, p2, jnew, head = 0; let done; if (xi == null || pstack == null) return -1; xi[0] = j; // initialize the recursion stack while (head >= 0) { j = xi[head]; // get j from the top of the recursion stack jnew = pinv != null ? (pinv[j]) : j; if (!(Gp[j] < 0)) { //CS_MARK(Gp, j); Gp[j] = -(Gp[j]) - 2; // mark node j as visited pstack[offset + head] = (jnew < 0) ? 0 : // CS_UNFLIP(Gp[jnew]); ((Gp[jnew] < 0) ? -Gp[jnew] - 2 : Gp[jnew]); } done = true; // node j done if no unvisited neighbors p2 = (jnew < 0) ? 0 : // CS_UNFLIP(Gp[jnew + 1]); ((Gp[jnew + 1] < 0) ? -Gp[jnew + 1] - 2 : Gp[jnew + 1]); for (p = pstack[offset + head]; p < p2; p++) // examine all neighbors of j { i = Gi[p]; // consider neighbor node i if (Gp[i] < 0) continue; // skip visited node i pstack[offset + head] = p; // pause depth-first search of node j xi[++head] = i; // start dfs at node i done = false; // node j is not done break; // break, to start dfs (i) } if (done) // depth-first search at node j is done { head--; // remove j from the recursion stack xi[--top] = j; // and place in the output stack } } return top; } static EliminationTree(m, n, colptr, rowind, ata) { let i, k, p, inext; let parent = new Int32Array(n); // allocate result let ancestor = new Int32Array(n); let prev = null; if (ata) { prev = new Int32Array(m); for (i = 0; i < m; i++) prev[i] = -1; } for (k = 0; k < n; k++) { parent[k] = -1; // node k has no parent yet ancestor[k] = -1; // nor does k have an ancestor for (p = colptr[k]; p < colptr[k + 1]; p++) { i = ata ? (prev[rowind[p]]) : (rowind[p]); for (; i != -1 && i < k; i = inext) // traverse from i to k { inext = ancestor[i]; // inext = ancestor of i ancestor[i] = k; // path compression if (inext == -1) parent[i] = k; // no anc., parent is k } if (ata) prev[rowind[p]] = k; } } return parent; } static TreePostorder(parent, n) { let j, k = 0; if (parent == null) return null; // check inputs let post = new Int32Array(n); // allocate result let w = new Int32Array(n); // get workspace let next = new Int32Array(n); let stack = new Int32Array(n); for (j = 0; j < n; j++) w[j] = -1; // empty linked lists for (j = n - 1; j >= 0; j--) // traverse nodes in reverse order*/ { if (parent[j] == -1) continue; // j is a root next[j] = w[parent[j]]; // add j to list of its parent w[parent[j]] = j; } for (j = 0; j < n; j++) { if (parent[j] != -1) continue; // skip j if it is not a root k = this.TreeDepthFirstSearch(j, k, w, next, post, stack); } return post; // success; free w, return post } static ColumnCounts(A, parent, post, ata) { let i, j, k, J, p, q, jleaf = 0; let ATp, ATi, colcount, delta, head = null, next = null; if (parent == null || post == null) return (null); // check inputs let m = A.RowCount; let n = A.ColumnCount; delta = colcount = new Float32Array(n); // allocate result let AT = A.Transpose(); // AT = A' // w is ancestor let w = new Int32Array(n); // get workspace let maxfirst = new Int32Array(n); let prevleaf = new Int32Array(n); let first = new Int32Array(n); for (k = 0; k < n; k++) { w[k] = -1; // clear workspace w [0..s-1] } maxfirst.set(w); prevleaf.set(w); first.set(w); for (k = 0; k < n; k++) // find first [j] { j = post[k]; delta[j] = (first[j] == -1) ? 1 : 0; // delta[j]=1 if j is a leaf for (; j != -1 && first[j] == -1; j = parent[j]) { first[j] = k; } } ATp = AT.ColumnPointers; ATi = AT.RowIndices; if (ata) // Init ata { head = new Int32Array(n + 1); next = new Int32Array(m); head.set(w); head[n] = -1; for (k = 0; k < n; k++) { w[post[k]] = k; // invert post } for (i = 0; i < m; i++) { for (k = n, p = ATp[i]; p < ATp[i + 1]; p++) { k = Math.min(k, w[ATi[p]]); } next[i] = head[k]; // place row i in linked list k head[k] = i; } } for (i = 0; i < n; i++) w[i] = i; // each node in its own set for (k = 0; k < n; k++) { j = post[k]; // j is the kth node in postordered etree if (parent[j] != -1) delta[parent[j]]--; // j is not a root //int HEAD(k,j) (ata ? head [k] : j) //int NEXT(J) (ata ? next [J] : -1) for (J = (ata ? head[k] : j); J != -1; J = (ata ? next[J] : -1)) // J=j for LL'=A case { for (p = ATp[J]; p < ATp[J + 1]; p++) { i = ATi[p]; q, jleaf = this.IsLeaf(i, j, first, maxfirst, prevleaf, w, jleaf); if (jleaf >= 1) delta[j]++; // A(i,j) is in skeleton if (jleaf == 2) delta[q]--; // account for overlap in q } } if (parent[j] != -1) w[j] = parent[j]; } for (j = 0; j < n; j++) // sum up delta's of each child { if (parent[j] != -1) colcount[parent[j]] += colcount[j]; } return colcount; // success: free workspace } static IsLeaf(i, j, first, maxfirst, prevleaf, ancestor, jleaf) { let q, s, sparent, jprev; if (first == null || maxfirst == null || prevleaf == null || ancestor == null) { return (-1, jleaf); } jleaf = 0; if (i <= j || first[j] <= maxfirst[i]) { return (-1, jleaf); // j not a leaf } maxfirst[i] = first[j]; // update max first[j] seen so far jprev = prevleaf[i]; // jprev = previous leaf of ith subtree prevleaf[i] = j; jleaf = (jprev == -1) ? 1 : 2; // j is first or subsequent leaf if (jleaf == 1) return (i, jleaf); // if 1st leaf, q = root of ith subtree for (q = jprev; q != ancestor[q]; q = ancestor[q]); for (s = jprev; s != q; s = sparent) { sparent = ancestor[s]; // path compression ancestor[s] = q; } return (q, jleaf); } static TreeDepthFirstSearch(j, k, head, next, post, stack) { let i, p, top = 0; if (head == null || next == null || post == null || stack == null) return (-1); // check inputs stack[0] = j; // place j on the stack while (top >= 0) // while (stack is not empty) { p = stack[top]; // p = top of stack i = head[p]; // i = youngest child of p if (i == -1) { top--; // p has no unordered children left post[k++] = p; // node p is the kth postordered node } else { head[p] = next[i]; // remove i from children of p top++; stack[top] = i; // start dfs on child node i } } return (k); } static Reach(Gp, Gi, Bp, Bi, n, k, xi, pinv) { if (xi == null) return (-1); // check inputs let p, top = n; for (p = Bp[k]; p < Bp[k + 1]; p++) { //if (!CS_MARKED(Gp, Bi[p])) if (!(Gp[Bi[p]] < 0)) // start a dfs at unmarked node i { top = this.DepthFirstSearch(Bi[p], Gp, Gi, top, xi, xi, n, pinv); } } for (p = top; p < n; p++) { //CS_MARK(Gp, xi[p]); Gp[xi[p]] = -(Gp[xi[p]]) - 2; // restore G } return top; } static EtreeReach(A, k, parent, s, w) { let i, p, n, len; if (parent == null || s == null || w == null) return -1; // check inputs let top = n = A.ColumnCount; let Ap = A.ColumnPointers; let Ai = A.RowIndices; //CS_MARK(w, k); w[k] = -w[k] - 2; // mark node k as visited for (p = Ap[k]; p < Ap[k + 1]; p++) { i = Ai[p]; // A(i,k) is nonzero if (i > k) { continue; // only use upper triangular part of A } for (len = 0; !(w[i] < 0); i = parent[i]) // traverse up etree { s[len++] = i; // L(k,i) is nonzero //CS_MARK(w, i); w[i] = -w[i] - 2; // mark i as visited } while (len > 0) { s[--top] = s[--len]; // push path onto stack } } for (p = top; p < n; p++) { //CS_MARK(w, s[p]); w[s[p]] = -w[s[p]] - 2; // unmark all nodes } //CS_MARK(w, k); w[k] = -w[k] - 2; // unmark node k return top; // s [top..n-1] contains pattern of L(k,:) } }37:T539,class SolverHelper { static SolveLower(L, x) { let p, j, k, n = L.ColumnCount; let lp = L.ColumnPointers; let li = L.RowIndices; let lx = L.Values; for (j = 0; j < n; j++) { x[j] /= lx[lp[j]]; k = lp[j + 1]; for (p = lp[j] + 1; p < k; p++) { x[li[p]] -= lx[p] * x[j]; } } } static SolveLowerTranspose(L, x) { let p, j, k, n = L.ColumnCount; let lp = L.ColumnPointers; let li = L.RowIndices; let lx = L.Values; for (j = n - 1; j >= 0; j--) { k = lp[j + 1]; for (p = lp[j] + 1; p < k; p++) { x[j] -= lx[p] * x[li[p]]; } x[j] /= lx[lp[j]]; } } static SolveUpper(U, x) { let p, j, k, n = U.ColumnCount; let up = U.ColumnPointers; let ui = U.RowIndices; let ux = U.Values; for (j = n - 1; j >= 0; j--) { x[j] /= ux[up[j + 1] - 1]; k = up[j + 1] - 1; for (p = up[j]; p < k; p++) { x[ui[p]] -= ux[p] * x[j]; } } } static SolveUpperTranspose(U, x) { let p, j, k, n = U.ColumnCount; let up = U.ColumnPointers; let ui = U.RowIndices; let ux = U.Values; for (j = 0; j < n; j++) { k = up[j + 1] - 1; for (p = up[j]; p < k; p++) { x[j] -= ux[p] * x[ui[p]]; } x[j] /= ux[up[j + 1] - 1]; } } }38:T3cd2,// TODO: Check this class class AMD { static Generate(A, order) { let Cp, Ci, P, W, nv, next, head, elen, degree, w, hhead; let d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1, k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n; let ok; let h; n = A.ColumnCount; if (order == ColumnOrdering.Natural) { // TODO: return null here? return Permutation.Create(n); } let C = this.ConstructMatrix(SymbolicColumnStorage.Create(A), order); Cp = C.ColumnPointers; cnz = Cp[n]; // Find dense threshold dense = Math.max(16, 10 * Math.trunc(Math.sqrt(n))); dense = Math.min(n - 2, dense); // add elbow room to C if (!C.Resize(cnz + cnz / 5 + 2 * n)) return null; P = new Int32Array(n + 1); // allocate result W = new Int32Array(n + 1); // get workspace w = new Int32Array(n + 1); degree = new Int32Array(n + 1); elen = new Int32Array(n + 1); // Initialized to 0's // Initialize quotient graph for (k = 0; k < n; k++) { W[k] = Cp[k + 1] - Cp[k]; } W[n] = 0; nzmax = C.RowIndices.Length; Ci = C.RowIndices; for (i = 0; i <= n; i++) { P[i] = -1; w[i] = 1; // node i is alive degree[i] = W[i]; // degree of node i } next = new Int32Array(P); hhead = new Int32Array(P); head = new Int32Array(P); nv = new Int32Array(w); mark = this.Clear(0, 0, w, n); // clear w elen[n] = -2; // n is a dead element Cp[n] = -1; // n is a root of assembly tree w[n] = 0; // n is a dead element // Initialize degree lists for (i = 0; i < n; i++) { d = degree[i]; if (d == 0) // node i is empty { elen[i] = -2; // element i is dead nel++; Cp[i] = -1; // i is a root of assembly tree w[i] = 0; } else if (d > dense) // node i is dense { nv[i] = 0; // absorb i into element n elen[i] = -1; // node i is dead nel++; Cp[i] = -(n + 2); // FLIP(n) nv[n]++; } else { if (head[d] != -1) P[head[d]] = i; next[i] = head[d]; // put node i in degree list d head[d] = i; } } while (nel < n) // while (selecting pivots) do { // Select node of minimum approximate degree for (k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) ; if (next[k] != -1) P[next[k]] = -1; head[mindeg] = next[k]; // remove k from degree list elenk = elen[k]; // elenk = |Ek| nvk = nv[k]; // # of nodes k represents nel += nvk; // nv[k] nodes of A eliminated // Garbage collection if (elenk > 0 && cnz + mindeg >= nzmax) { for (j = 0; j < n; j++) { if ((p = Cp[j]) >= 0) // j is a live node or element { Cp[j] = Ci[p]; // save first entry of object Ci[p] = -(j + 2); // first entry is now CS_FLIP(j) } } for (q = 0, p = 0; p < cnz;) // scan all of memory { if ((j = this.FLIP(Ci[p++])) >= 0) // found object j { Ci[q] = Cp[j]; // restore first entry of object Cp[j] = q++; // new pointer to object j for (k3 = 0; k3 < W[j] - 1; k3++) { Ci[q++] = Ci[p++]; } } } cnz = q; // Ci [cnz...nzmax-1] now free } // Construct new element dk = 0; nv[k] = -nvk; // flag k as in Lk p = Cp[k]; pk1 = (elenk == 0) ? p : cnz; // do in place if elen[k] == 0 pk2 = pk1; for (k1 = 1; k1 <= elenk + 1; k1++) { if (k1 > elenk) { e = k; // search the nodes in k pj = p; // list of nodes starts at Ci[pj]*/ ln = W[k] - elenk; // length of list of nodes in k } else { e = Ci[p++]; // search the nodes in e pj = Cp[e]; ln = W[e]; // length of list of nodes in e } for (k2 = 1; k2 <= ln; k2++) { i = Ci[pj++]; if ((nvi = nv[i]) <= 0) continue; // node i dead, or seen dk += nvi; // degree[Lk] += size of node i nv[i] = -nvi; // negate nv[i] to denote i in Lk Ci[pk2++] = i; // place i in Lk if (next[i] != -1) P[next[i]] = P[i]; if (P[i] != -1) // remove i from degree list { next[P[i]] = next[i]; } else { head[degree[i]] = next[i]; } } if (e != k) { Cp[e] = -(k + 2); // absorb e into k // FLIP(k) w[e] = 0; // e is now a dead element } } if (elenk != 0) cnz = pk2; // Ci [cnz...nzmax] is free degree[k] = dk; // external degree of k - |Lk\i| Cp[k] = pk1; // element k is in Ci[pk1..pk2-1] W[k] = pk2 - pk1; elen[k] = -2; // k is now an element // Find set differences mark = this.Clear(mark, lemax, w, n); // clear w if necessary for (pk = pk1; pk < pk2; pk++) // scan 1: find |Le\Lk| { i = Ci[pk]; if ((eln = elen[i]) <= 0) continue; // skip if elen[i] empty nvi = -nv[i]; // nv [i] was negated wnvi = mark - nvi; for (p = Cp[i]; p <= Cp[i] + eln - 1; p++) // scan Ei { e = Ci[p]; if (w[e] >= mark) { w[e] -= nvi; // decrement |Le\Lk| } else if (w[e] != 0) // ensure e is a live element { w[e] = degree[e] + wnvi; // 1st time e seen in scan 1 } } } // Degree update for (pk = pk1; pk < pk2; pk++) // scan2: degree update { i = Ci[pk]; // consider node i in Lk p1 = Cp[i]; p2 = p1 + elen[i] - 1; pn = p1; for (h = 0, d = 0, p = p1; p <= p2; p++) // scan Ei { e = Ci[p]; if (w[e] != 0) // e is an unabsorbed element { dext = w[e] - mark; // dext = |Le\Lk| if (dext > 0) { d += dext; // sum up the set differences Ci[pn++] = e; // keep e in Ei h += e; // compute the hash of node i } else { Cp[e] = -(k + 2); // aggressive absorb. e.k // FLIP(k) w[e] = 0; // e is a dead element } } } elen[i] = pn - p1 + 1; // elen[i] = |Ei| p3 = pn; p4 = p1 + W[i]; for (p = p2 + 1; p < p4; p++) // prune edges in Ai { j = Ci[p]; if ((nvj = nv[j]) <= 0) continue; // node j dead or in Lk d += nvj; // degree(i) += |j| Ci[pn++] = j; // place j in node list of i h += j; // compute hash for node i } if (d == 0) // check for mass elimination { Cp[i] = -(k + 2); // absorb i into k // FLIP(k) nvi = -nv[i]; dk -= nvi; // |Lk| -= |i| nvk += nvi; // |k| += nv[i] nel += nvi; nv[i] = 0; elen[i] = -1; // node i is dead } else { degree[i] = Math.min(degree[i], d); // update degree(i) Ci[pn] = Ci[p3]; // move first node to end Ci[p3] = Ci[p1]; // move 1st el. to end of Ei Ci[p1] = k; // add k as 1st element in of Ei W[i] = pn - p1 + 1; // new len of adj. list of node i h = ((h < 0) ? (-h) : h) % n; // finalize hash of i next[i] = hhead[h]; // place i in hash bucket hhead[h] = i; P[i] = h; // save hash of i in last[i] } } // scan2 is done degree[k] = dk; // finalize |Lk| lemax = Math.max(lemax, dk); mark = this.Clear(mark + lemax, lemax, w, n); // clear w // Supernode detection for (pk = pk1; pk < pk2; pk++) { i = Ci[pk]; if (nv[i] >= 0) continue; // skip if i is dead h = P[i]; // scan hash bucket of node i i = hhead[h]; hhead[h] = -1; // hash bucket will be empty for (; i != -1 && next[i] != -1; i = next[i], mark++) { ln = W[i]; eln = elen[i]; for (p = Cp[i] + 1; p <= Cp[i] + ln - 1; p++) w[Ci[p]] = mark; jlast = i; for (j = next[i]; j != -1;) // compare i with all j { ok = (W[j] == ln) && (elen[j] == eln); for (p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) { if (w[Ci[p]] != mark) ok = false; // compare i and j } if (ok) // i and j are identical { Cp[j] = -(i + 2); // absorb j into i // FLIP(i) nv[i] += nv[j]; nv[j] = 0; elen[j] = -1; // node j is dead j = next[j]; // delete j from hash bucket next[jlast] = j; } else { jlast = j; // j and i are different j = next[j]; } } } } // Finalize new element for (p = pk1, pk = pk1; pk < pk2; pk++) // finalize Lk { i = Ci[pk]; if ((nvi = -nv[i]) <= 0) continue; // skip if i is dead nv[i] = nvi; // restore nv[i] d = degree[i] + dk - nvi; // compute external degree(i) d = Math.min(d, n - nel - nvi); if (head[d] != -1) P[head[d]] = i; next[i] = head[d]; // put i back in degree list P[i] = -1; head[d] = i; mindeg = Math.min(mindeg, d); // find new minimum degree degree[i] = d; Ci[p++] = i; // place i in Lk } nv[k] = nvk; // # nodes absorbed into k if ((W[k] = p - pk1) == 0) // length of adj list of element k { Cp[k] = -1; // k is a root of the tree w[k] = 0; // k is now a dead element } if (elenk != 0) cnz = p; // free unused space in Lk } // Postordering for (i = 0; i < n; i++) Cp[i] = -(Cp[i] + 2); // fix assembly tree // FLIP(Cp[i]) for (j = 0; j <= n; j++) head[j] = -1; for (j = n; j >= 0; j--) // place unordered nodes in lists { if (nv[j] > 0) continue; // skip if j is an element next[j] = head[Cp[j]]; // place j in list of its parent head[Cp[j]] = j; } for (e = n; e >= 0; e--) // place elements in lists { if (nv[e] <= 0) continue; // skip unless e is an element if (Cp[e] != -1) { next[e] = head[Cp[e]]; // place e in list of its parent head[Cp[e]] = e; } } for (k = 0, i = 0; i <= n; i++) // postorder the assembly tree { if (Cp[i] == -1) { k = GraphHelper.TreeDepthFirstSearch(i, k, head, next, P, w); } } return P; } static FLIP(i) { return -i - 2; } // clear w static Clear(mark, lemax, w, n) { let k; if (mark < 2 || (mark + lemax < 0)) { for (k = 0; k < n; k++) { if (w[k] != 0) { w[k] = 1; } } mark = 2; } return mark; // at this point, w [0..n-1] < mark holds } // keep off-diagonal entries; drop diagonal entries static KeepOffDiag(i, j) { return (i != j); } // Construct matrix C static ConstructMatrix(A, order) { let result = null; // Compute A' let AT = A.Transpose(); let m = A.RowCount; let n = A.ColumnCount; if (order == ColumnOrdering.MinimumDegreeAtPlusA) { if (n != m) { throw "Resources.MatrixSquare A"; } // Return A+A' result = A.Add(AT); } else if (order == ColumnOrdering.MinimumDegreeStS) { // Drop dense columns from AT let dense, p, p2 = 0; // Find dense threshold dense = Math.max(16, 10 * Math.trunc(Math.sqrt(n))); dense = Math.min(n - 2, dense); let colptr = AT.ColumnPointers; let rowind = AT.RowIndices; for (let j = 0; j < m; j++) { // Column j of AT starts here. p = colptr[j]; // New column j starts here. colptr[j] = p2; if (colptr[j + 1] - p > dense) { // Skip dense column j continue; } for (; p < colptr[j + 1]; p++) { rowind[p2++] = rowind[p]; } } colptr[m] = p2; // Return A'*A with no dense rows result = AT.Multiply(AT.Transpose()); } else { // Return A'*A result = AT.Multiply(A); } // Drop diagonal entries. result.Keep(this.KeepOffDiag); return result; } }39:[] 2f:{"title":"Sparse Linear Algebra","cells":[{"id":0,"value":"md`# Sparse Linear Algebra`","pinned":false,"hidden":false,"mode":"ojs"},{"id":394,"value":"md`Sparse linear algebra in pure JavaScript.\n\n*References:*\n- [CSparse.NET](https://github.com/wo80/CSparse.NET)\n- *Direct Methods for Sparse Linear Systems, SIAM, Philadelphia, PA, 2006 by Tim Davis*`","pinned":false,"hidden":false,"mode":"ojs"},{"id":412,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":240,"value":"md`## Examples`","pinned":false,"hidden":false,"mode":"ojs"},{"id":447,"value":"md`### Create a sparse matrix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":488,"value":"A = {\n // size of the squared matrix\n const n = 40;\n \n // create csc matrix from column pointers (ia), row indices (ja), and values (a)\n return new CompressedColumnStorage(n, n, a, ja, ia);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":491,"value":"md`### Plot sparse structure`","pinned":false,"hidden":false,"mode":"ojs"},{"id":438,"value":"{ \n // plot sparse structure\n return sparseStructurePlot(A);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":464,"value":"sparseStructurePlot = (A, unit = 10) => {\n const ctx = DOM.context2d(A.RowCount * unit, A.ColumnCount * unit);\n \n for (let [row, col, _] of A.EnumerateIndexed()) {\n ctx.fillRect(row * unit, col * unit, unit, unit);\n }\n \n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":441,"value":"md`### SparseLU Solver`","pinned":false,"hidden":false,"mode":"ojs"},{"id":504,"value":"md`Initialize Solver`","pinned":false,"hidden":false,"mode":"ojs"},{"id":499,"value":"sparseLU = {\n // size of the squared matrix\n let n = 40;\n \n // create csc matrix from column pointers (ia), row indices (ja), and values (a)\n let A = new CompressedColumnStorage(n, n, a, ja, ia);\n \n // initialize solver\n return SparseLU.Create(A);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":508,"value":"md`Lower`","pinned":false,"hidden":false,"mode":"ojs"},{"id":510,"value":"sparseStructurePlot(sparseLU.U); // check original implementation","pinned":false,"hidden":false,"mode":"ojs"},{"id":518,"value":"md`Upper`","pinned":false,"hidden":false,"mode":"ojs"},{"id":513,"value":"sparseStructurePlot(sparseLU.L); // check original implementation","pinned":false,"hidden":false,"mode":"ojs"},{"id":522,"value":"md`Check results`","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"{\n // size of the squared matrix\n let n = 40;\n \n // initialize solution\n let x = new Float64Array(n);\n \n // solve A.x = b\n sparseLU.Solve(b, x);\n \n // check result\n return test(\"Solution is correct\", () => {\n for (let i = 0; i < n; i++) {\n if (assertAlmostEqual(x[i] - xExpected[i], 1e-10)) {\n throw `invalid result at position ${i}:\\n actual = ${x[i]}\\n expected = ${xExpected[i]}\\n error = ${x[i] - xExpected[i]}`;\n }\n }\n });\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":456,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"ia = new Int32Array([0, 13, 26, 39, 50, 64, 77, 90, 112, 127, 145, 160, 179, 194, 209, 222, 239, 257, 270, 282, 296, 311, 322, 335, 346, 360, 369, 384, 396, 407, 425, 438, 449, 460, 473, 487, 498, 509, 523, 537, 550]);","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"b = new Float64Array([9.597773156848987, 12.497482264478945, 12.733251045250235, 6.356058789165047, 14.445331778142895, 12.292862649056847, 13.708320393688217, 9.57041653616123, 10.654345998869106, 8.436660298839161, 11.364646202254415, 11.773693047012246, 6.349697409694839, 10.138160327263032, 9.822266966833693, 12.023820335999808, 12.942922901669931, 11.777190222778742, 15.792966796717185, 9.995897423354668, 6.940160733691441, 12.92014229238994, 9.948963320682866, 7.9259990545075105, 11.024778963018479, 10.256768029111687, 7.97352619251819, 8.280025929486266, 10.347912485383768, 10.08187832110433, 8.969518363736242, 11.082687545164479, 12.605874448437115, 10.122010951440961, 8.998042400806032, 10.409428786688997, 8.727692433661325, 9.762796828042688, 15.952252144580681, 9.878602133374466]);","pinned":false,"hidden":false,"mode":"ojs"},{"id":150,"value":"xExpected = new Float64Array([1, 1.0256410256410255, 1.0512820512820513, 1.0769230769230769, 1.1025641025641026, 1.1282051282051282, 1.1538461538461537, 1.1794871794871795, 1.205128205128205, 1.2307692307692308, 1.2564102564102564, 1.282051282051282, 1.3076923076923077, 1.3333333333333333, 1.358974358974359, 1.3846153846153846, 1.4102564102564101, 1.435897435897436, 1.4615384615384617, 1.4871794871794872, 1.5128205128205128, 1.5384615384615383, 1.564102564102564, 1.5897435897435899, 1.6153846153846154, 1.641025641025641, 1.6666666666666665, 1.6923076923076923, 1.717948717948718, 1.7435897435897436, 1.7692307692307692, 1.7948717948717947, 1.8205128205128205, 1.8461538461538463, 1.8717948717948718, 1.8974358974358974, 1.9230769230769231, 1.9487179487179487, 1.9743589743589745, 2]);","pinned":false,"hidden":false,"mode":"ojs"},{"id":68,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":53,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":70,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":74,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":90,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":203,"value":"ColumnOrdering = {\n return {\n Natural: 0,\n MinimumDegreeAtPlusA: 1,\n MinimumDegreeStS: 2,\n MinimumDegreeAtA: 3\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":98,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":115,"value":"class SymbolicFactorization\n{\n constructor()\n {\n this.pinv = new Int32Array(0); // inverse row perm. for QR, fill red. perm for Chol\n this.q = new Int32Array(0); // fill-reducing column permutation for LU and QR\n this.parent = new Int32Array(0); // elimination tree for Cholesky and QR\n this.cp = new Int32Array(0); // column pointers for Cholesky, row counts for QR\n this.leftmost = new Int32Array(0); // leftmost[i] = min(find(A(i,:))), for QR\n this.m2 = 0; // # of rows for QR, after adding fictitious rows\n this.lnz = 0; // # entries in L for LU or Cholesky; in V for QR\n this.unz = 0; // # entries in U for LU; in R for QR\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":95,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":206,"value":"class SymbolicColumnStorage\n{\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":361,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":214,"value":"md`**ToDo:**`","pinned":false,"hidden":false,"mode":"ojs"},{"id":243,"value":"md`\n- Add Multiply, Transpose, etc.\n- JS conform naming (lowerCase)\n- Better error handling\n- Implement AMD and SymbolicColumnStorage\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":272,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":251,"value":"test = (title, testFunction) => {\n testFunction();\n return html `✔️ ${title||'Test passing '}`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":249,"value":"assertAlmostEqual = (a, b, epsilon) => {\n return Math.abs(b, a) < epsilon\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":284,"value":"assert = chai.assert;","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"should = chai.should();","pinned":false,"hidden":false,"mode":"ojs"},{"id":281,"value":"chai = require('https://bundle.run/chai@4.2.0')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W39","stdlib":"1","titling":"auto"}