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