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