|
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 { |
|
5797
|
217 info = -1; |
|
5610
|
218 #ifdef HAVE_CXSPARSE |
|
|
219 octave_idx_type nr = a.rows(); |
|
|
220 octave_idx_type nc = a.cols(); |
|
|
221 octave_idx_type b_nc = b.cols(); |
|
|
222 octave_idx_type b_nr = b.rows(); |
|
|
223 const double *bvec = b.fortran_vec(); |
|
|
224 Matrix x; |
|
|
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 ()) |
|
5797
|
233 return Matrix(); |
|
5610
|
234 x.resize(nc, b_nc); |
|
|
235 double *vec = x.fortran_vec(); |
|
|
236 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
237 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; |
|
|
238 i++, idx+=nc, bidx+=b_nr) |
|
|
239 { |
|
|
240 OCTAVE_QUIT; |
|
5681
|
241 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
242 buf[j] = 0.; |
|
5610
|
243 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
244 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
245 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr); |
|
|
246 #else |
|
5648
|
247 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); |
|
5792
|
248 #endif |
|
5610
|
249 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
250 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
251 { |
|
|
252 OCTAVE_QUIT; |
|
|
253 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
254 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
255 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
256 } |
|
|
257 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
258 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
259 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
260 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, vec + idx, nc); |
|
|
261 #else |
|
5648
|
262 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); |
|
5792
|
263 #endif |
|
5610
|
264 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
265 } |
|
5797
|
266 info = 0; |
|
5610
|
267 } |
|
|
268 else |
|
|
269 { |
|
|
270 SparseMatrix at = a.hermitian(); |
|
5792
|
271 SparseQR q (at, 3); |
|
5610
|
272 if (! q.ok ()) |
|
5797
|
273 return Matrix(); |
|
5610
|
274 x.resize(nc, b_nc); |
|
|
275 double *vec = x.fortran_vec(); |
|
5681
|
276 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
|
277 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
278 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; |
|
|
279 i++, idx+=nc, bidx+=b_nr) |
|
|
280 { |
|
|
281 OCTAVE_QUIT; |
|
5681
|
282 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
283 buf[j] = 0.; |
|
5610
|
284 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
285 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
286 CXSPARSE_DNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr); |
|
|
287 #else |
|
5648
|
288 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); |
|
5792
|
289 #endif |
|
5648
|
290 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
291 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
292 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
293 { |
|
|
294 OCTAVE_QUIT; |
|
|
295 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
296 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
297 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
298 } |
|
|
299 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
300 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
301 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc); |
|
|
302 #else |
|
5648
|
303 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); |
|
5792
|
304 #endif |
|
5610
|
305 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
306 } |
|
5797
|
307 info = 0; |
|
5610
|
308 } |
|
|
309 |
|
|
310 return x; |
|
|
311 #else |
|
|
312 return Matrix (); |
|
|
313 #endif |
|
|
314 } |
|
|
315 |
|
|
316 SparseMatrix |
|
|
317 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info) |
|
|
318 { |
|
5797
|
319 info = -1; |
|
5610
|
320 #ifdef HAVE_CXSPARSE |
|
|
321 octave_idx_type nr = a.rows(); |
|
|
322 octave_idx_type nc = a.cols(); |
|
|
323 octave_idx_type b_nr = b.rows(); |
|
|
324 octave_idx_type b_nc = b.cols(); |
|
|
325 SparseMatrix x; |
|
|
326 volatile octave_idx_type ii, x_nz; |
|
|
327 |
|
|
328 if (nr < 1 || nc < 1 || nr != b_nr) |
|
|
329 (*current_liboctave_error_handler) |
|
|
330 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
|
331 else if (nr >= nc) |
|
|
332 { |
|
5792
|
333 SparseQR q (a, 3); |
|
5610
|
334 if (! q.ok ()) |
|
5797
|
335 return SparseMatrix(); |
|
5610
|
336 x = SparseMatrix (nc, b_nc, b.nzmax()); |
|
|
337 x.xcidx(0) = 0; |
|
|
338 x_nz = b.nzmax(); |
|
|
339 ii = 0; |
|
|
340 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
341 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
342 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
343 { |
|
|
344 OCTAVE_QUIT; |
|
|
345 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
346 Xx[j] = b.xelem(j,i); |
|
5681
|
347 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
348 buf[j] = 0.; |
|
5610
|
349 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
350 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
351 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); |
|
|
352 #else |
|
5648
|
353 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); |
|
5792
|
354 #endif |
|
5610
|
355 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
356 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
357 { |
|
|
358 OCTAVE_QUIT; |
|
|
359 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
360 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
361 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
362 } |
|
|
363 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
364 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
365 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
366 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); |
|
|
367 #else |
|
5648
|
368 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); |
|
5792
|
369 #endif |
|
5610
|
370 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
371 |
|
|
372 for (octave_idx_type j = 0; j < nc; j++) |
|
|
373 { |
|
|
374 double tmp = Xx[j]; |
|
|
375 if (tmp != 0.0) |
|
|
376 { |
|
|
377 if (ii == x_nz) |
|
|
378 { |
|
|
379 // Resize the sparse matrix |
|
|
380 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
381 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
382 x.change_capacity (sz); |
|
|
383 x_nz = sz; |
|
|
384 } |
|
|
385 x.xdata(ii) = tmp; |
|
|
386 x.xridx(ii++) = j; |
|
|
387 } |
|
|
388 } |
|
|
389 x.xcidx(i+1) = ii; |
|
|
390 } |
|
5797
|
391 info = 0; |
|
5610
|
392 } |
|
|
393 else |
|
|
394 { |
|
|
395 SparseMatrix at = a.hermitian(); |
|
5792
|
396 SparseQR q (at, 3); |
|
5610
|
397 if (! q.ok ()) |
|
5797
|
398 return SparseMatrix(); |
|
5610
|
399 x = SparseMatrix (nc, b_nc, b.nzmax()); |
|
|
400 x.xcidx(0) = 0; |
|
|
401 x_nz = b.nzmax(); |
|
|
402 ii = 0; |
|
5681
|
403 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
5610
|
404 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
5681
|
405 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
406 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
407 { |
|
|
408 OCTAVE_QUIT; |
|
|
409 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
410 Xx[j] = b.xelem(j,i); |
|
5681
|
411 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
412 buf[j] = 0.; |
|
5610
|
413 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
414 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
415 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); |
|
|
416 #else |
|
5648
|
417 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); |
|
5792
|
418 #endif |
|
5648
|
419 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
420 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
421 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
422 { |
|
|
423 OCTAVE_QUIT; |
|
|
424 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
425 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
426 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
427 } |
|
|
428 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
429 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
430 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); |
|
|
431 #else |
|
5648
|
432 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); |
|
5792
|
433 #endif |
|
5610
|
434 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
435 |
|
|
436 for (octave_idx_type j = 0; j < nc; j++) |
|
|
437 { |
|
|
438 double tmp = Xx[j]; |
|
|
439 if (tmp != 0.0) |
|
|
440 { |
|
|
441 if (ii == x_nz) |
|
|
442 { |
|
|
443 // Resize the sparse matrix |
|
|
444 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
445 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
446 x.change_capacity (sz); |
|
|
447 x_nz = sz; |
|
|
448 } |
|
|
449 x.xdata(ii) = tmp; |
|
|
450 x.xridx(ii++) = j; |
|
|
451 } |
|
|
452 } |
|
|
453 x.xcidx(i+1) = ii; |
|
|
454 } |
|
5797
|
455 info = 0; |
|
5610
|
456 } |
|
|
457 |
|
|
458 x.maybe_compress (); |
|
|
459 return x; |
|
|
460 #else |
|
|
461 return SparseMatrix (); |
|
|
462 #endif |
|
|
463 } |
|
|
464 |
|
|
465 ComplexMatrix |
|
|
466 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info) |
|
|
467 { |
|
5797
|
468 info = -1; |
|
5610
|
469 #ifdef HAVE_CXSPARSE |
|
|
470 octave_idx_type nr = a.rows(); |
|
|
471 octave_idx_type nc = a.cols(); |
|
|
472 octave_idx_type b_nc = b.cols(); |
|
|
473 octave_idx_type b_nr = b.rows(); |
|
|
474 ComplexMatrix x; |
|
|
475 |
|
|
476 if (nr < 1 || nc < 1 || nr != b_nr) |
|
|
477 (*current_liboctave_error_handler) |
|
|
478 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
|
479 else if (nr >= nc) |
|
|
480 { |
|
5792
|
481 SparseQR q (a, 3); |
|
5610
|
482 if (! q.ok ()) |
|
5797
|
483 return ComplexMatrix(); |
|
5610
|
484 x.resize(nc, b_nc); |
|
|
485 Complex *vec = x.fortran_vec(); |
|
|
486 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
487 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
|
488 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
489 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
490 { |
|
|
491 OCTAVE_QUIT; |
|
|
492 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
493 { |
|
|
494 Complex c = b.xelem (j,i); |
|
|
495 Xx[j] = std::real (c); |
|
|
496 Xz[j] = std::imag (c); |
|
|
497 } |
|
5681
|
498 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
499 buf[j] = 0.; |
|
5610
|
500 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
501 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
502 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); |
|
|
503 #else |
|
5648
|
504 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); |
|
5792
|
505 #endif |
|
5610
|
506 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
507 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
508 { |
|
|
509 OCTAVE_QUIT; |
|
|
510 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
511 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
512 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
513 } |
|
|
514 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
515 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
516 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
517 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); |
|
|
518 #else |
|
5648
|
519 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); |
|
5792
|
520 #endif |
|
5681
|
521 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
522 buf[j] = 0.; |
|
5792
|
523 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
524 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); |
|
|
525 #else |
|
5648
|
526 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); |
|
5792
|
527 #endif |
|
5610
|
528 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
529 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
530 { |
|
|
531 OCTAVE_QUIT; |
|
|
532 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
533 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
534 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
535 } |
|
|
536 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
537 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
538 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
539 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); |
|
|
540 #else |
|
5648
|
541 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); |
|
5792
|
542 #endif |
|
5610
|
543 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
544 for (octave_idx_type j = 0; j < nc; j++) |
|
|
545 vec[j+idx] = Complex (Xx[j], Xz[j]); |
|
|
546 } |
|
5797
|
547 info = 0; |
|
5610
|
548 } |
|
|
549 else |
|
|
550 { |
|
|
551 SparseMatrix at = a.hermitian(); |
|
5792
|
552 SparseQR q (at, 3); |
|
5610
|
553 if (! q.ok ()) |
|
5797
|
554 return ComplexMatrix(); |
|
5610
|
555 x.resize(nc, b_nc); |
|
|
556 Complex *vec = x.fortran_vec(); |
|
5681
|
557 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
5610
|
558 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
559 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
5681
|
560 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
561 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
562 { |
|
|
563 OCTAVE_QUIT; |
|
|
564 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
565 { |
|
|
566 Complex c = b.xelem (j,i); |
|
|
567 Xx[j] = std::real (c); |
|
|
568 Xz[j] = std::imag (c); |
|
|
569 } |
|
5681
|
570 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
571 buf[j] = 0.; |
|
5610
|
572 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
573 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
574 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); |
|
|
575 #else |
|
5648
|
576 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); |
|
5792
|
577 #endif |
|
5648
|
578 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
579 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
580 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
581 { |
|
|
582 OCTAVE_QUIT; |
|
|
583 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
584 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
585 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
586 } |
|
|
587 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
588 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
589 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); |
|
|
590 #else |
|
5648
|
591 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); |
|
5792
|
592 #endif |
|
5681
|
593 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
594 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
595 buf[j] = 0.; |
|
|
596 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
597 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
598 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); |
|
|
599 #else |
|
5648
|
600 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); |
|
5792
|
601 #endif |
|
5648
|
602 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
603 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
604 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
605 { |
|
|
606 OCTAVE_QUIT; |
|
|
607 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
608 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
609 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
610 } |
|
|
611 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
612 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
613 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); |
|
|
614 #else |
|
5648
|
615 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); |
|
5792
|
616 #endif |
|
5610
|
617 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
618 for (octave_idx_type j = 0; j < nc; j++) |
|
|
619 vec[j+idx] = Complex (Xx[j], Xz[j]); |
|
|
620 } |
|
5797
|
621 info = 0; |
|
5610
|
622 } |
|
|
623 |
|
|
624 return x; |
|
|
625 #else |
|
|
626 return ComplexMatrix (); |
|
|
627 #endif |
|
|
628 } |
|
|
629 |
|
|
630 SparseComplexMatrix |
|
|
631 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info) |
|
|
632 { |
|
5797
|
633 info = -1; |
|
5610
|
634 #ifdef HAVE_CXSPARSE |
|
|
635 octave_idx_type nr = a.rows(); |
|
|
636 octave_idx_type nc = a.cols(); |
|
|
637 octave_idx_type b_nr = b.rows(); |
|
|
638 octave_idx_type b_nc = b.cols(); |
|
|
639 SparseComplexMatrix x; |
|
|
640 volatile octave_idx_type ii, x_nz; |
|
|
641 |
|
|
642 if (nr < 1 || nc < 1 || nr != b_nr) |
|
|
643 (*current_liboctave_error_handler) |
|
|
644 ("matrix dimension mismatch in solution of minimum norm problem"); |
|
|
645 else if (nr >= nc) |
|
|
646 { |
|
5792
|
647 SparseQR q (a, 3); |
|
5610
|
648 if (! q.ok ()) |
|
5797
|
649 return SparseComplexMatrix(); |
|
5610
|
650 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); |
|
|
651 x.xcidx(0) = 0; |
|
|
652 x_nz = b.nzmax(); |
|
|
653 ii = 0; |
|
|
654 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
655 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
|
656 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); |
|
|
657 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
658 { |
|
|
659 OCTAVE_QUIT; |
|
|
660 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
661 { |
|
|
662 Complex c = b.xelem (j,i); |
|
|
663 Xx[j] = std::real (c); |
|
|
664 Xz[j] = std::imag (c); |
|
|
665 } |
|
5681
|
666 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
667 buf[j] = 0.; |
|
5610
|
668 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
669 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
670 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); |
|
|
671 #else |
|
5648
|
672 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); |
|
5792
|
673 #endif |
|
5610
|
674 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
675 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
676 { |
|
|
677 OCTAVE_QUIT; |
|
|
678 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
679 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
680 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
681 } |
|
|
682 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
683 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
684 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
685 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); |
|
|
686 #else |
|
5648
|
687 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); |
|
5792
|
688 #endif |
|
5681
|
689 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
690 for (octave_idx_type j = nr; j < q.S()->m2; j++) |
|
|
691 buf[j] = 0.; |
|
|
692 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
693 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
694 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); |
|
|
695 #else |
|
5648
|
696 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); |
|
5792
|
697 #endif |
|
5610
|
698 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
699 for (volatile octave_idx_type j = 0; j < nc; j++) |
|
|
700 { |
|
|
701 OCTAVE_QUIT; |
|
|
702 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
703 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
704 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
705 } |
|
|
706 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
707 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); |
|
5792
|
708 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
709 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); |
|
|
710 #else |
|
5648
|
711 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); |
|
5792
|
712 #endif |
|
5610
|
713 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
714 |
|
|
715 for (octave_idx_type j = 0; j < nc; j++) |
|
|
716 { |
|
|
717 Complex tmp = Complex (Xx[j], Xz[j]); |
|
|
718 if (tmp != 0.0) |
|
|
719 { |
|
|
720 if (ii == x_nz) |
|
|
721 { |
|
|
722 // Resize the sparse matrix |
|
|
723 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
724 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
725 x.change_capacity (sz); |
|
|
726 x_nz = sz; |
|
|
727 } |
|
|
728 x.xdata(ii) = tmp; |
|
|
729 x.xridx(ii++) = j; |
|
|
730 } |
|
|
731 } |
|
|
732 x.xcidx(i+1) = ii; |
|
|
733 } |
|
5797
|
734 info = 0; |
|
5610
|
735 } |
|
|
736 else |
|
|
737 { |
|
|
738 SparseMatrix at = a.hermitian(); |
|
5792
|
739 SparseQR q (at, 3); |
|
5610
|
740 if (! q.ok ()) |
|
5797
|
741 return SparseComplexMatrix(); |
|
5610
|
742 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); |
|
|
743 x.xcidx(0) = 0; |
|
|
744 x_nz = b.nzmax(); |
|
|
745 ii = 0; |
|
5681
|
746 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); |
|
5610
|
747 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); |
|
|
748 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); |
|
5681
|
749 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); |
|
5610
|
750 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) |
|
|
751 { |
|
|
752 OCTAVE_QUIT; |
|
|
753 for (octave_idx_type j = 0; j < b_nr; j++) |
|
|
754 { |
|
|
755 Complex c = b.xelem (j,i); |
|
|
756 Xx[j] = std::real (c); |
|
|
757 Xz[j] = std::imag (c); |
|
|
758 } |
|
5681
|
759 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
760 buf[j] = 0.; |
|
5610
|
761 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
762 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
763 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); |
|
|
764 #else |
|
5648
|
765 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); |
|
5792
|
766 #endif |
|
5648
|
767 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
768 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
769 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
770 { |
|
|
771 OCTAVE_QUIT; |
|
|
772 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
773 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
774 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
775 } |
|
|
776 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
777 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
778 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); |
|
|
779 #else |
|
5648
|
780 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); |
|
5792
|
781 #endif |
|
5681
|
782 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
783 for (octave_idx_type j = nr; j < nbuf; j++) |
|
|
784 buf[j] = 0.; |
|
|
785 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
786 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
787 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); |
|
|
788 #else |
|
5648
|
789 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); |
|
5792
|
790 #endif |
|
5648
|
791 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); |
|
5610
|
792 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
793 for (volatile octave_idx_type j = nr-1; j >= 0; j--) |
|
|
794 { |
|
|
795 OCTAVE_QUIT; |
|
|
796 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5648
|
797 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); |
|
5610
|
798 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
799 } |
|
|
800 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
5792
|
801 #if defined(CS_VER) && (CS_VER >= 2) |
|
|
802 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); |
|
|
803 #else |
|
5648
|
804 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); |
|
5792
|
805 #endif |
|
5610
|
806 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; |
|
|
807 |
|
|
808 for (octave_idx_type j = 0; j < nc; j++) |
|
|
809 { |
|
|
810 Complex tmp = Complex (Xx[j], Xz[j]); |
|
|
811 if (tmp != 0.0) |
|
|
812 { |
|
|
813 if (ii == x_nz) |
|
|
814 { |
|
|
815 // Resize the sparse matrix |
|
|
816 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; |
|
|
817 sz = (sz > 10 ? sz : 10) + x_nz; |
|
|
818 x.change_capacity (sz); |
|
|
819 x_nz = sz; |
|
|
820 } |
|
|
821 x.xdata(ii) = tmp; |
|
|
822 x.xridx(ii++) = j; |
|
|
823 } |
|
|
824 } |
|
|
825 x.xcidx(i+1) = ii; |
|
|
826 } |
|
5797
|
827 info = 0; |
|
5610
|
828 } |
|
|
829 |
|
|
830 x.maybe_compress (); |
|
|
831 return x; |
|
|
832 #else |
|
|
833 return SparseComplexMatrix (); |
|
|
834 #endif |
|
|
835 } |
|
|
836 |
|
|
837 /* |
|
|
838 ;;; Local Variables: *** |
|
|
839 ;;; mode: C++ *** |
|
|
840 ;;; End: *** |
|
|
841 */ |
