|
5610
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2005 David Bateman |
|
|
4 |
|
7016
|
5 This file is part of Octave. |
|
|
6 |
|
5610
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
7016
|
9 Free Software Foundation; either version 3 of the License, or (at your |
|
|
10 option) any later version. |
|
5610
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
7016
|
18 along with Octave; see the file COPYING. If not, see |
|
|
19 <http://www.gnu.org/licenses/>. |
|
5610
|
20 |
|
|
21 */ |
|
|
22 |
|
|
23 #ifdef HAVE_CONFIG_H |
|
|
24 #include <config.h> |
|
|
25 #endif |
|
|
26 #include <vector> |
|
|
27 |
|
|
28 #include "lo-error.h" |
|
|
29 #include "SparseQR.h" |
|
|
30 |
|
|
31 SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order) |
|
|
32 { |
|
|
33 #ifdef HAVE_CXSPARSE |
|
5648
|
34 CXSPARSE_DNAME () A; |
|
5610
|
35 A.nzmax = a.nzmax (); |
|
|
36 A.m = a.rows (); |
|
|
37 A.n = a.cols (); |
|
|
38 nrows = A.m; |
|
|
39 // Cast away const on A, with full knowledge that CSparse won't touch it |
|
|
40 // Prevents the methods below making a copy of the data. |
|
|
41 A.p = const_cast<octave_idx_type *>(a.cidx ()); |
|
|
42 A.i = const_cast<octave_idx_type *>(a.ridx ()); |
|
|
43 A.x = const_cast<double *>(a.data ()); |
|
|
44 A.nz = -1; |
|
|
45 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
46 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
47 S = CXSPARSE_DNAME (_sqr) (order, &A, 1); |
|
|
48 #else |
|
|
49 S = CXSPARSE_DNAME (_sqr) (&A, order - 1, 1); |
|
|
50 #endif |
|
|
51 |
|
5648
|
52 N = CXSPARSE_DNAME (_qr) (&A, S); |
|
5610
|
53 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
54 if (!N) |
|
|
55 (*current_liboctave_error_handler) |
|
|
56 ("SparseQR: sparse matrix QR factorization filled"); |
|
|
57 count = 1; |
|
|
58 #else |
|
|
59 (*current_liboctave_error_handler) |
|
|
60 ("SparseQR: sparse matrix QR factorization not implemented"); |
|
|
61 #endif |
|
|
62 } |
|
|
63 |
|
|
64 SparseQR::SparseQR_rep::~SparseQR_rep (void) |
|
|
65 { |
|
|
66 #ifdef HAVE_CXSPARSE |
|
5648
|
67 CXSPARSE_DNAME (_sfree) (S); |
|
|
68 CXSPARSE_DNAME (_nfree) (N); |
|
5610
|
69 #endif |
|
|
70 } |
|
|
71 |
|
|
72 SparseMatrix |
|
|
73 SparseQR::SparseQR_rep::V (void) const |
|
|
74 { |
|
|
75 #ifdef HAVE_CXSPARSE |
|
|
76 // Drop zeros from V and sort |
|
5775
|
77 // FIXME Is the double transpose to sort necessary? |
|
5610
|
78 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
79 CXSPARSE_DNAME (_dropzeros) (N->L); |
|
|
80 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1); |
|
|
81 CXSPARSE_DNAME (_spfree) (N->L); |
|
|
82 N->L = CXSPARSE_DNAME (_transpose) (D, 1); |
|
|
83 CXSPARSE_DNAME (_spfree) (D); |
|
5610
|
84 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
85 |
|
|
86 octave_idx_type nc = N->L->n; |
|
|
87 octave_idx_type nz = N->L->nzmax; |
|
|
88 SparseMatrix ret (N->L->m, nc, nz); |
|
|
89 for (octave_idx_type j = 0; j < nc+1; j++) |
|
|
90 ret.xcidx (j) = N->L->p[j]; |
|
|
91 for (octave_idx_type j = 0; j < nz; j++) |
|
|
92 { |
|
|
93 ret.xridx (j) = N->L->i[j]; |
|
|
94 ret.xdata (j) = N->L->x[j]; |
|
|
95 } |
|
|
96 return ret; |
|
|
97 #else |
|
|
98 return SparseMatrix (); |
|
|
99 #endif |
|
|
100 } |
|
|
101 |
|
|
102 ColumnVector |
|
|
103 SparseQR::SparseQR_rep::Pinv (void) const |
|
|
104 { |
|
|
105 #ifdef HAVE_CXSPARSE |
|
|
106 ColumnVector ret(N->L->m); |
|
|
107 for (octave_idx_type i = 0; i < N->L->m; i++) |
|
5792
|
108 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
109 ret.xelem(i) = S->pinv[i]; |
|
|
110 #else |
|
5610
|
111 ret.xelem(i) = S->Pinv[i]; |
|
5792
|
112 #endif |
|
5610
|
113 return ret; |
|
|
114 #else |
|
|
115 return ColumnVector (); |
|
|
116 #endif |
|
|
117 } |
|
|
118 |
|
|
119 ColumnVector |
|
|
120 SparseQR::SparseQR_rep::P (void) const |
|
|
121 { |
|
|
122 #ifdef HAVE_CXSPARSE |
|
|
123 ColumnVector ret(N->L->m); |
|
|
124 for (octave_idx_type i = 0; i < N->L->m; i++) |
|
5792
|
125 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
126 ret.xelem(S->pinv[i]) = i; |
|
|
127 #else |
|
5610
|
128 ret.xelem(S->Pinv[i]) = i; |
|
5792
|
129 #endif |
|
5610
|
130 return ret; |
|
|
131 #else |
|
|
132 return ColumnVector (); |
|
|
133 #endif |
|
|
134 } |
|
|
135 |
|
|
136 SparseMatrix |
|
|
137 SparseQR::SparseQR_rep::R (const bool econ) const |
|
|
138 { |
|
|
139 #ifdef HAVE_CXSPARSE |
|
|
140 // Drop zeros from R and sort |
|
5775
|
141 // FIXME Is the double transpose to sort necessary? |
|
5610
|
142 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
143 CXSPARSE_DNAME (_dropzeros) (N->U); |
|
|
144 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1); |
|
|
145 CXSPARSE_DNAME (_spfree) (N->U); |
|
|
146 N->U = CXSPARSE_DNAME (_transpose) (D, 1); |
|
|
147 CXSPARSE_DNAME (_spfree) (D); |
|
5610
|
148 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
149 |
|
|
150 octave_idx_type nc = N->U->n; |
|
|
151 octave_idx_type nz = N->U->nzmax; |
|
|
152 |
|
|
153 SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); |
|
|
154 |
|
|
155 for (octave_idx_type j = 0; j < nc+1; j++) |
|
|
156 ret.xcidx (j) = N->U->p[j]; |
|
|
157 for (octave_idx_type j = 0; j < nz; j++) |
|
|
158 { |
|
|
159 ret.xridx (j) = N->U->i[j]; |
|
|
160 ret.xdata (j) = N->U->x[j]; |
|
|
161 } |
|
|
162 return ret; |
|
|
163 #else |
|
|
164 return SparseMatrix (); |
|
|
165 #endif |
|
|
166 } |
|
|
167 |
|
|
168 Matrix |
|
|
169 SparseQR::SparseQR_rep::C (const Matrix &b) const |
|
|
170 { |
|
|
171 #ifdef HAVE_CXSPARSE |
|
|
172 octave_idx_type b_nr = b.rows(); |
|
|
173 octave_idx_type b_nc = b.cols(); |
|
|
174 octave_idx_type nc = N->L->n; |
|
|
175 octave_idx_type nr = nrows; |
|
|
176 const double *bvec = b.fortran_vec(); |
|
6924
|
177 Matrix ret (b_nr, b_nc); |
|
5610
|
178 double *vec = ret.fortran_vec(); |
|
6924
|
179 if (nr < 0 || nc < 0 || nr != b_nr) |
|
5610
|
180 (*current_liboctave_error_handler) ("matrix dimension mismatch"); |
|
6924
|
181 else if (nr == 0 || nc == 0 || b_nc == 0) |
|
|
182 ret = Matrix (nc, b_nc, 0.0); |
|
5610
|
183 else |
|
|
184 { |
|
|
185 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); |
|
|
186 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) |
|
|
187 { |
|
|
188 OCTAVE_QUIT; |
|
5681
|
189 for (octave_idx_type i = nr; i < S->m2; i++) |
|
|
190 buf[i] = 0.; |
|
5610
|
191 volatile octave_idx_type nm = (nr < nc ? nr : nc); |
|
|
192 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
193 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
194 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr); |
|
|
195 #else |
|
5648
|
196 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf); |
|
5792
|
197 #endif |
|
5610
|
198 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
199 |
|
|
200 for (volatile octave_idx_type i = 0; i < nm; i++) |
|
|
201 { |
|
|
202 OCTAVE_QUIT; |
|
|
203 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
204 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); |
|
5610
|
205 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
206 } |
|
|
207 for (octave_idx_type i = 0; i < b_nr; i++) |
|
|
208 vec[i+idx] = buf[i]; |
|
|
209 } |
|
|
210 } |
|
|
211 return ret; |
|
|
212 #else |
|
|
213 return Matrix (); |
|
|
214 #endif |
|
|
215 } |
|
|
216 |
|
|
217 Matrix |
|
|
218 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info) |
|
|
219 { |
|
5797
|
220 info = -1; |
|
5610
|
221 #ifdef HAVE_CXSPARSE |
|
|
222 octave_idx_type nr = a.rows(); |
|
|
223 octave_idx_type nc = a.cols(); |
|
|
224 octave_idx_type b_nc = b.cols(); |
|
|
225 octave_idx_type b_nr = b.rows(); |
|
|
226 const double *bvec = b.fortran_vec(); |
|
|
227 Matrix x; |
|
|
228 |
|
6924
|
229 if (nr < 0 || nc < 0 || nr != b_nr) |
|
5610
|
230 (*current_liboctave_error_handler) |
|
|
231 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
6924
|
232 else if (nr == 0 || nc == 0 || b_nc == 0) |
|
|
233 x = Matrix (nc, b_nc, 0.0); |
|
5610
|
234 else if (nr >= nc) |
|
|
235 { |
|
5792
|
236 SparseQR q (a, 3); |
|
5610
|
237 if (! q.ok ()) |
|
5797
|
238 return Matrix(); |
|
5610
|
239 x.resize(nc, b_nc); |
|
|
240 double *vec = x.fortran_vec(); |
|
|
241 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
242 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; |
|
|
243 i++, idx+=nc, bidx+=b_nr) |
|
|
244 { |
|
|
245 OCTAVE_QUIT; |
|
5681
|
246 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
247 buf[j] = 0.; |
|
5610
|
248 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
249 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
250 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr); |
|
|
251 #else |
|
5648
|
252 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); |
|
5792
|
253 #endif |
|
5610
|
254 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
255 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
256 { |
|
|
257 OCTAVE_QUIT; |
|
|
258 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
259 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
260 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
261 } |
|
|
262 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
263 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
264 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
265 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, vec + idx, nc); |
|
|
266 #else |
|
5648
|
267 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); |
|
5792
|
268 #endif |
|
5610
|
269 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
270 } |
|
5797
|
271 info = 0; |
|
5610
|
272 } |
|
|
273 else |
|
|
274 { |
|
|
275 SparseMatrix at = a.hermitian(); |
|
5792
|
276 SparseQR q (at, 3); |
|
5610
|
277 if (! q.ok ()) |
|
5797
|
278 return Matrix(); |
|
5610
|
279 x.resize(nc, b_nc); |
|
|
280 double *vec = x.fortran_vec(); |
|
5681
|
281 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
|
282 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
283 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; |
|
|
284 i++, idx+=nc, bidx+=b_nr) |
|
|
285 { |
|
|
286 OCTAVE_QUIT; |
|
5681
|
287 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
288 buf[j] = 0.; |
|
5610
|
289 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
290 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
291 CXSPARSE_DNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr); |
|
|
292 #else |
|
5648
|
293 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); |
|
5792
|
294 #endif |
|
5648
|
295 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
296 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
297 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
298 { |
|
|
299 OCTAVE_QUIT; |
|
|
300 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
301 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
302 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
303 } |
|
|
304 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
305 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
306 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc); |
|
|
307 #else |
|
5648
|
308 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); |
|
5792
|
309 #endif |
|
5610
|
310 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
311 } |
|
5797
|
312 info = 0; |
|
5610
|
313 } |
|
|
314 |
|
|
315 return x; |
|
|
316 #else |
|
|
317 return Matrix (); |
|
|
318 #endif |
|
|
319 } |
|
|
320 |
|
|
321 SparseMatrix |
|
|
322 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info) |
|
|
323 { |
|
5797
|
324 info = -1; |
|
5610
|
325 #ifdef HAVE_CXSPARSE |
|
|
326 octave_idx_type nr = a.rows(); |
|
|
327 octave_idx_type nc = a.cols(); |
|
|
328 octave_idx_type b_nr = b.rows(); |
|
|
329 octave_idx_type b_nc = b.cols(); |
|
|
330 SparseMatrix x; |
|
|
331 volatile octave_idx_type ii, x_nz; |
|
|
332 |
|
6924
|
333 if (nr < 0 || nc < 0 || nr != b_nr) |
|
5610
|
334 (*current_liboctave_error_handler) |
|
|
335 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
6924
|
336 else if (nr == 0 || nc == 0 || b_nc == 0) |
|
|
337 x = SparseMatrix (nc, b_nc); |
|
5610
|
338 else if (nr >= nc) |
|
|
339 { |
|
5792
|
340 SparseQR q (a, 3); |
|
5610
|
341 if (! q.ok ()) |
|
5797
|
342 return SparseMatrix(); |
|
5610
|
343 x = SparseMatrix (nc, b_nc, b.nzmax()); |
|
|
344 x.xcidx(0) = 0; |
|
|
345 x_nz = b.nzmax(); |
|
|
346 ii = 0; |
|
|
347 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
348 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
349 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
350 { |
|
|
351 OCTAVE_QUIT; |
|
|
352 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
353 Xx[j] = b.xelem(j,i); |
|
5681
|
354 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
355 buf[j] = 0.; |
|
5610
|
356 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
357 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
358 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); |
|
|
359 #else |
|
5648
|
360 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); |
|
5792
|
361 #endif |
|
5610
|
362 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
363 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
364 { |
|
|
365 OCTAVE_QUIT; |
|
|
366 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
367 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
368 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
369 } |
|
|
370 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
371 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
372 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
373 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); |
|
|
374 #else |
|
5648
|
375 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); |
|
5792
|
376 #endif |
|
5610
|
377 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
378 |
|
|
379 for (octave_idx_type j = 0; j < nc; j++) |
|
|
380 { |
|
|
381 double tmp = Xx[j]; |
|
|
382 if (tmp != 0.0) |
|
|
383 { |
|
|
384 if (ii == x_nz) |
|
|
385 { |
|
|
386 // Resize the sparse matrix |
|
|
387 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
388 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
389 x.change_capacity (sz); |
|
|
390 x_nz = sz; |
|
|
391 } |
|
|
392 x.xdata(ii) = tmp; |
|
|
393 x.xridx(ii++) = j; |
|
|
394 } |
|
|
395 } |
|
|
396 x.xcidx(i+1) = ii; |
|
|
397 } |
|
5797
|
398 info = 0; |
|
5610
|
399 } |
|
|
400 else |
|
|
401 { |
|
|
402 SparseMatrix at = a.hermitian(); |
|
5792
|
403 SparseQR q (at, 3); |
|
5610
|
404 if (! q.ok ()) |
|
5797
|
405 return SparseMatrix(); |
|
5610
|
406 x = SparseMatrix (nc, b_nc, b.nzmax()); |
|
|
407 x.xcidx(0) = 0; |
|
|
408 x_nz = b.nzmax(); |
|
|
409 ii = 0; |
|
5681
|
410 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
5610
|
411 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
5681
|
412 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
413 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
414 { |
|
|
415 OCTAVE_QUIT; |
|
|
416 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
417 Xx[j] = b.xelem(j,i); |
|
5681
|
418 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
419 buf[j] = 0.; |
|
5610
|
420 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
421 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
422 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); |
|
|
423 #else |
|
5648
|
424 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); |
|
5792
|
425 #endif |
|
5648
|
426 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
427 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
428 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
429 { |
|
|
430 OCTAVE_QUIT; |
|
|
431 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
432 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
433 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
434 } |
|
|
435 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
436 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
437 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); |
|
|
438 #else |
|
5648
|
439 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); |
|
5792
|
440 #endif |
|
5610
|
441 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
442 |
|
|
443 for (octave_idx_type j = 0; j < nc; j++) |
|
|
444 { |
|
|
445 double tmp = Xx[j]; |
|
|
446 if (tmp != 0.0) |
|
|
447 { |
|
|
448 if (ii == x_nz) |
|
|
449 { |
|
|
450 // Resize the sparse matrix |
|
|
451 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
452 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
453 x.change_capacity (sz); |
|
|
454 x_nz = sz; |
|
|
455 } |
|
|
456 x.xdata(ii) = tmp; |
|
|
457 x.xridx(ii++) = j; |
|
|
458 } |
|
|
459 } |
|
|
460 x.xcidx(i+1) = ii; |
|
|
461 } |
|
5797
|
462 info = 0; |
|
5610
|
463 } |
|
|
464 |
|
|
465 x.maybe_compress (); |
|
|
466 return x; |
|
|
467 #else |
|
|
468 return SparseMatrix (); |
|
|
469 #endif |
|
|
470 } |
|
|
471 |
|
|
472 ComplexMatrix |
|
|
473 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info) |
|
|
474 { |
|
5797
|
475 info = -1; |
|
5610
|
476 #ifdef HAVE_CXSPARSE |
|
|
477 octave_idx_type nr = a.rows(); |
|
|
478 octave_idx_type nc = a.cols(); |
|
|
479 octave_idx_type b_nc = b.cols(); |
|
|
480 octave_idx_type b_nr = b.rows(); |
|
|
481 ComplexMatrix x; |
|
|
482 |
|
6924
|
483 if (nr < 0 || nc < 0 || nr != b_nr) |
|
5610
|
484 (*current_liboctave_error_handler) |
|
|
485 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
6924
|
486 else if (nr == 0 || nc == 0 || b_nc == 0) |
|
|
487 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); |
|
5610
|
488 else if (nr >= nc) |
|
|
489 { |
|
5792
|
490 SparseQR q (a, 3); |
|
5610
|
491 if (! q.ok ()) |
|
5797
|
492 return ComplexMatrix(); |
|
5610
|
493 x.resize(nc, b_nc); |
|
|
494 Complex *vec = x.fortran_vec(); |
|
|
495 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
496 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
|
497 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
498 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
499 { |
|
|
500 OCTAVE_QUIT; |
|
|
501 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
502 { |
|
|
503 Complex c = b.xelem (j,i); |
|
|
504 Xx[j] = std::real (c); |
|
|
505 Xz[j] = std::imag (c); |
|
|
506 } |
|
5681
|
507 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
508 buf[j] = 0.; |
|
5610
|
509 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
510 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
511 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); |
|
|
512 #else |
|
5648
|
513 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); |
|
5792
|
514 #endif |
|
5610
|
515 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
516 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
517 { |
|
|
518 OCTAVE_QUIT; |
|
|
519 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
520 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
521 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
522 } |
|
|
523 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
524 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
525 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
526 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); |
|
|
527 #else |
|
5648
|
528 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); |
|
5792
|
529 #endif |
|
5681
|
530 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
531 buf[j] = 0.; |
|
5792
|
532 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
533 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); |
|
|
534 #else |
|
5648
|
535 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); |
|
5792
|
536 #endif |
|
5610
|
537 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
538 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
539 { |
|
|
540 OCTAVE_QUIT; |
|
|
541 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
542 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
543 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
544 } |
|
|
545 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
546 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
547 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
548 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); |
|
|
549 #else |
|
5648
|
550 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); |
|
5792
|
551 #endif |
|
5610
|
552 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
553 for (octave_idx_type j = 0; j < nc; j++) |
|
|
554 vec[j+idx] = Complex (Xx[j], Xz[j]); |
|
|
555 } |
|
5797
|
556 info = 0; |
|
5610
|
557 } |
|
|
558 else |
|
|
559 { |
|
|
560 SparseMatrix at = a.hermitian(); |
|
5792
|
561 SparseQR q (at, 3); |
|
5610
|
562 if (! q.ok ()) |
|
5797
|
563 return ComplexMatrix(); |
|
5610
|
564 x.resize(nc, b_nc); |
|
|
565 Complex *vec = x.fortran_vec(); |
|
5681
|
566 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
5610
|
567 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
568 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
5681
|
569 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
570 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
571 { |
|
|
572 OCTAVE_QUIT; |
|
|
573 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
574 { |
|
|
575 Complex c = b.xelem (j,i); |
|
|
576 Xx[j] = std::real (c); |
|
|
577 Xz[j] = std::imag (c); |
|
|
578 } |
|
5681
|
579 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
580 buf[j] = 0.; |
|
5610
|
581 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
582 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
583 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); |
|
|
584 #else |
|
5648
|
585 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); |
|
5792
|
586 #endif |
|
5648
|
587 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
588 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
589 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
590 { |
|
|
591 OCTAVE_QUIT; |
|
|
592 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
593 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
594 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
595 } |
|
|
596 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
597 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
598 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); |
|
|
599 #else |
|
5648
|
600 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); |
|
5792
|
601 #endif |
|
5681
|
602 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
603 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
604 buf[j] = 0.; |
|
|
605 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
606 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
607 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); |
|
|
608 #else |
|
5648
|
609 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); |
|
5792
|
610 #endif |
|
5648
|
611 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
612 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
613 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
614 { |
|
|
615 OCTAVE_QUIT; |
|
|
616 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
617 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
618 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
619 } |
|
|
620 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
621 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
622 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); |
|
|
623 #else |
|
5648
|
624 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); |
|
5792
|
625 #endif |
|
5610
|
626 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
627 for (octave_idx_type j = 0; j < nc; j++) |
|
|
628 vec[j+idx] = Complex (Xx[j], Xz[j]); |
|
|
629 } |
|
5797
|
630 info = 0; |
|
5610
|
631 } |
|
|
632 |
|
|
633 return x; |
|
|
634 #else |
|
|
635 return ComplexMatrix (); |
|
|
636 #endif |
|
|
637 } |
|
|
638 |
|
|
639 SparseComplexMatrix |
|
|
640 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info) |
|
|
641 { |
|
5797
|
642 info = -1; |
|
5610
|
643 #ifdef HAVE_CXSPARSE |
|
|
644 octave_idx_type nr = a.rows(); |
|
|
645 octave_idx_type nc = a.cols(); |
|
|
646 octave_idx_type b_nr = b.rows(); |
|
|
647 octave_idx_type b_nc = b.cols(); |
|
|
648 SparseComplexMatrix x; |
|
|
649 volatile octave_idx_type ii, x_nz; |
|
|
650 |
|
6924
|
651 if (nr < 0 || nc < 0 || nr != b_nr) |
|
5610
|
652 (*current_liboctave_error_handler) |
|
|
653 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
6924
|
654 else if (nr == 0 || nc == 0 || b_nc == 0) |
|
|
655 x = SparseComplexMatrix (nc, b_nc); |
|
5610
|
656 else if (nr >= nc) |
|
|
657 { |
|
5792
|
658 SparseQR q (a, 3); |
|
5610
|
659 if (! q.ok ()) |
|
5797
|
660 return SparseComplexMatrix(); |
|
5610
|
661 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); |
|
|
662 x.xcidx(0) = 0; |
|
|
663 x_nz = b.nzmax(); |
|
|
664 ii = 0; |
|
|
665 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
666 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
|
667 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
668 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
669 { |
|
|
670 OCTAVE_QUIT; |
|
|
671 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
672 { |
|
|
673 Complex c = b.xelem (j,i); |
|
|
674 Xx[j] = std::real (c); |
|
|
675 Xz[j] = std::imag (c); |
|
|
676 } |
|
5681
|
677 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
678 buf[j] = 0.; |
|
5610
|
679 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
680 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
681 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); |
|
|
682 #else |
|
5648
|
683 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); |
|
5792
|
684 #endif |
|
5610
|
685 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
686 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
687 { |
|
|
688 OCTAVE_QUIT; |
|
|
689 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
690 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
691 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
692 } |
|
|
693 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
694 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
695 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
696 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); |
|
|
697 #else |
|
5648
|
698 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); |
|
5792
|
699 #endif |
|
5681
|
700 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
701 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
702 buf[j] = 0.; |
|
|
703 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
704 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
705 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); |
|
|
706 #else |
|
5648
|
707 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); |
|
5792
|
708 #endif |
|
5610
|
709 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
710 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
711 { |
|
|
712 OCTAVE_QUIT; |
|
|
713 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
714 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
715 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
716 } |
|
|
717 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
718 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
719 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
720 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); |
|
|
721 #else |
|
5648
|
722 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); |
|
5792
|
723 #endif |
|
5610
|
724 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
725 |
|
|
726 for (octave_idx_type j = 0; j < nc; j++) |
|
|
727 { |
|
|
728 Complex tmp = Complex (Xx[j], Xz[j]); |
|
|
729 if (tmp != 0.0) |
|
|
730 { |
|
|
731 if (ii == x_nz) |
|
|
732 { |
|
|
733 // Resize the sparse matrix |
|
|
734 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
735 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
736 x.change_capacity (sz); |
|
|
737 x_nz = sz; |
|
|
738 } |
|
|
739 x.xdata(ii) = tmp; |
|
|
740 x.xridx(ii++) = j; |
|
|
741 } |
|
|
742 } |
|
|
743 x.xcidx(i+1) = ii; |
|
|
744 } |
|
5797
|
745 info = 0; |
|
5610
|
746 } |
|
|
747 else |
|
|
748 { |
|
|
749 SparseMatrix at = a.hermitian(); |
|
5792
|
750 SparseQR q (at, 3); |
|
5610
|
751 if (! q.ok ()) |
|
5797
|
752 return SparseComplexMatrix(); |
|
5610
|
753 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); |
|
|
754 x.xcidx(0) = 0; |
|
|
755 x_nz = b.nzmax(); |
|
|
756 ii = 0; |
|
5681
|
757 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
5610
|
758 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
759 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
5681
|
760 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
761 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
762 { |
|
|
763 OCTAVE_QUIT; |
|
|
764 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
765 { |
|
|
766 Complex c = b.xelem (j,i); |
|
|
767 Xx[j] = std::real (c); |
|
|
768 Xz[j] = std::imag (c); |
|
|
769 } |
|
5681
|
770 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
771 buf[j] = 0.; |
|
5610
|
772 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
773 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
774 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); |
|
|
775 #else |
|
5648
|
776 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); |
|
5792
|
777 #endif |
|
5648
|
778 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
779 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
780 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
781 { |
|
|
782 OCTAVE_QUIT; |
|
|
783 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
784 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
785 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
786 } |
|
|
787 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
788 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
789 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); |
|
|
790 #else |
|
5648
|
791 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); |
|
5792
|
792 #endif |
|
5681
|
793 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
794 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
795 buf[j] = 0.; |
|
|
796 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
797 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
798 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); |
|
|
799 #else |
|
5648
|
800 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); |
|
5792
|
801 #endif |
|
5648
|
802 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
803 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
804 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
805 { |
|
|
806 OCTAVE_QUIT; |
|
|
807 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
808 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
809 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
810 } |
|
|
811 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
812 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
813 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); |
|
|
814 #else |
|
5648
|
815 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); |
|
5792
|
816 #endif |
|
5610
|
817 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
818 |
|
|
819 for (octave_idx_type j = 0; j < nc; j++) |
|
|
820 { |
|
|
821 Complex tmp = Complex (Xx[j], Xz[j]); |
|
|
822 if (tmp != 0.0) |
|
|
823 { |
|
|
824 if (ii == x_nz) |
|
|
825 { |
|
|
826 // Resize the sparse matrix |
|
|
827 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
828 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
829 x.change_capacity (sz); |
|
|
830 x_nz = sz; |
|
|
831 } |
|
|
832 x.xdata(ii) = tmp; |
|
|
833 x.xridx(ii++) = j; |
|
|
834 } |
|
|
835 } |
|
|
836 x.xcidx(i+1) = ii; |
|
|
837 } |
|
5797
|
838 info = 0; |
|
5610
|
839 } |
|
|
840 |
|
|
841 x.maybe_compress (); |
|
|
842 return x; |
|
|
843 #else |
|
|
844 return SparseComplexMatrix (); |
|
|
845 #endif |
|
|
846 } |
|
|
847 |
|
|
848 /* |
|
|
849 ;;; Local Variables: *** |
|
|
850 ;;; mode: C++ *** |
|
|
851 ;;; End: *** |
|
|
852 */ |
