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