csSqr.js (4764B)
1 import { csPermute } from './csPermute.js'; 2 import { csPost } from './csPost.js'; 3 import { csEtree } from './csEtree.js'; 4 import { createCsAmd } from './csAmd.js'; 5 import { createCsCounts } from './csCounts.js'; 6 import { factory } from '../../../utils/factory.js'; 7 var name = 'csSqr'; 8 var dependencies = ['add', 'multiply', 'transpose']; 9 export var createCsSqr = /* #__PURE__ */factory(name, dependencies, _ref => { 10 var { 11 add, 12 multiply, 13 transpose 14 } = _ref; 15 var csAmd = createCsAmd({ 16 add, 17 multiply, 18 transpose 19 }); 20 var csCounts = createCsCounts({ 21 transpose 22 }); 23 /** 24 * Symbolic ordering and analysis for QR and LU decompositions. 25 * 26 * @param {Number} order The ordering strategy (see csAmd for more details) 27 * @param {Matrix} a The A matrix 28 * @param {boolean} qr Symbolic ordering and analysis for QR decomposition (true) or 29 * symbolic ordering and analysis for LU decomposition (false) 30 * 31 * @return {Object} The Symbolic ordering and analysis for matrix A 32 * 33 * Reference: http://faculty.cse.tamu.edu/davis/publications.html 34 */ 35 36 return function csSqr(order, a, qr) { 37 // a arrays 38 var aptr = a._ptr; 39 var asize = a._size; // columns 40 41 var n = asize[1]; // vars 42 43 var k; // symbolic analysis result 44 45 var s = {}; // fill-reducing ordering 46 47 s.q = csAmd(order, a); // validate results 48 49 if (order && !s.q) { 50 return null; 51 } // QR symbolic analysis 52 53 54 if (qr) { 55 // apply permutations if needed 56 var c = order ? csPermute(a, null, s.q, 0) : a; // etree of C'*C, where C=A(:,q) 57 58 s.parent = csEtree(c, 1); // post order elimination tree 59 60 var post = csPost(s.parent, n); // col counts chol(C'*C) 61 62 s.cp = csCounts(c, s.parent, post, 1); // check we have everything needed to calculate number of nonzero elements 63 64 if (c && s.parent && s.cp && _vcount(c, s)) { 65 // calculate number of nonzero elements 66 for (s.unz = 0, k = 0; k < n; k++) { 67 s.unz += s.cp[k]; 68 } 69 } 70 } else { 71 // for LU factorization only, guess nnz(L) and nnz(U) 72 s.unz = 4 * aptr[n] + n; 73 s.lnz = s.unz; 74 } // return result S 75 76 77 return s; 78 }; 79 /** 80 * Compute nnz(V) = s.lnz, s.pinv, s.leftmost, s.m2 from A and s.parent 81 */ 82 83 function _vcount(a, s) { 84 // a arrays 85 var aptr = a._ptr; 86 var aindex = a._index; 87 var asize = a._size; // rows & columns 88 89 var m = asize[0]; 90 var n = asize[1]; // initialize s arrays 91 92 s.pinv = []; // (m + n) 93 94 s.leftmost = []; // (m) 95 // vars 96 97 var parent = s.parent; 98 var pinv = s.pinv; 99 var leftmost = s.leftmost; // workspace, next: first m entries, head: next n entries, tail: next n entries, nque: next n entries 100 101 var w = []; // (m + 3 * n) 102 103 var next = 0; 104 var head = m; 105 var tail = m + n; 106 var nque = m + 2 * n; // vars 107 108 var i, k, p, p0, p1; // initialize w 109 110 for (k = 0; k < n; k++) { 111 // queue k is empty 112 w[head + k] = -1; 113 w[tail + k] = -1; 114 w[nque + k] = 0; 115 } // initialize row arrays 116 117 118 for (i = 0; i < m; i++) { 119 leftmost[i] = -1; 120 } // loop columns backwards 121 122 123 for (k = n - 1; k >= 0; k--) { 124 // values & index for column k 125 for (p0 = aptr[k], p1 = aptr[k + 1], p = p0; p < p1; p++) { 126 // leftmost[i] = min(find(A(i,:))) 127 leftmost[aindex[p]] = k; 128 } 129 } // scan rows in reverse order 130 131 132 for (i = m - 1; i >= 0; i--) { 133 // row i is not yet ordered 134 pinv[i] = -1; 135 k = leftmost[i]; // check row i is empty 136 137 if (k === -1) { 138 continue; 139 } // first row in queue k 140 141 142 if (w[nque + k]++ === 0) { 143 w[tail + k] = i; 144 } // put i at head of queue k 145 146 147 w[next + i] = w[head + k]; 148 w[head + k] = i; 149 } 150 151 s.lnz = 0; 152 s.m2 = m; // find row permutation and nnz(V) 153 154 for (k = 0; k < n; k++) { 155 // remove row i from queue k 156 i = w[head + k]; // count V(k,k) as nonzero 157 158 s.lnz++; // add a fictitious row 159 160 if (i < 0) { 161 i = s.m2++; 162 } // associate row i with V(:,k) 163 164 165 pinv[i] = k; // skip if V(k+1:m,k) is empty 166 167 if (--nque[k] <= 0) { 168 continue; 169 } // nque[k] is nnz (V(k+1:m,k)) 170 171 172 s.lnz += w[nque + k]; // move all rows to parent of k 173 174 var pa = parent[k]; 175 176 if (pa !== -1) { 177 if (w[nque + pa] === 0) { 178 w[tail + pa] = w[tail + k]; 179 } 180 181 w[next + w[tail + k]] = w[head + pa]; 182 w[head + pa] = w[next + i]; 183 w[nque + pa] += w[nque + k]; 184 } 185 } 186 187 for (i = 0; i < m; i++) { 188 if (pinv[i] < 0) { 189 pinv[i] = k++; 190 } 191 } 192 193 return true; 194 } 195 });