Mercurial > hg > octave-nkf

comparison liboctave/SparseCmplxQR.cc @ 7505:f5005d9510f4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove dispatched sparse functions and treat in the generic versions of the functions
author David Bateman <dbateman@free.fr>
date Wed, 20 Feb 2008 15:52:11 -0500
parents a1dbe9d80eee
children 82be108cc558
comparison
equal deleted inserted replaced
7504:ddcf233d765b 7505:f5005d9510f4
37 #define OCTAVE_C99_COMPLEX(buf, n) \ 37 #define OCTAVE_C99_COMPLEX(buf, n) \
38 OCTAVE_LOCAL_BUFFER (double, buf ## tmp, (2 * (n))); \ 38 OCTAVE_LOCAL_BUFFER (double, buf ## tmp, (2 * (n))); \
39 cs_complex_t *buf = reinterpret_cast<cs_complex_t *> (buf ## tmp); 39 cs_complex_t *buf = reinterpret_cast<cs_complex_t *> (buf ## tmp);
40 40
41 #define OCTAVE_C99_ZERO (0. + 0.iF) 41 #define OCTAVE_C99_ZERO (0. + 0.iF)
42 #define OCTAVE_C99_ONE (1. + 0.iF)
42 #else 43 #else
43 #define OCTAVE_C99_COMPLEX(buf, n) \ 44 #define OCTAVE_C99_COMPLEX(buf, n) \
44 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, (n)); 45 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, (n));
45 #define OCTAVE_C99_ZERO cs_complex_t(0., 0.); 46 #define OCTAVE_C99_ZERO cs_complex_t(0., 0.);
47 #define OCTAVE_C99_ONE cs_complex_t(1., 0.);
46 #endif 48 #endif
47 49
48 SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep 50 SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep
49 (GCC_ATTR_UNUSED const SparseComplexMatrix& a, GCC_ATTR_UNUSED int order) 51 (GCC_ATTR_UNUSED const SparseComplexMatrix& a, GCC_ATTR_UNUSED int order)
50 { 52 {
224 for (octave_idx_type i = 0; i < b_nr; i++) 226 for (octave_idx_type i = 0; i < b_nr; i++)
225 vec[i+idx] = buf[i]; 227 vec[i+idx] = buf[i];
226 } 228 }
227 } 229 }
228 return ret; 230 return ret;
231 #else
232 return ComplexMatrix ();
233 #endif
234 }
235
236 ComplexMatrix
237 SparseComplexQR::SparseComplexQR_rep::Q (void) const
238 {
239 #ifdef HAVE_CXSPARSE
240 octave_idx_type nc = N->L->n;
241 octave_idx_type nr = nrows;
242 ComplexMatrix ret(nr, nr);
243 Complex *vec = ret.fortran_vec();
244 if (nr < 0 || nc < 0)
245 (*current_liboctave_error_handler) ("matrix dimension mismatch");
246 else if (nr == 0 || nc == 0)
247 ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
248 else
249 {
250 OCTAVE_C99_COMPLEX (bvec, nr);
251 for (octave_idx_type i = 0; i < nr; i++)
252 bvec[i] = OCTAVE_C99_ZERO;
253 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
254 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
255 {
256 OCTAVE_QUIT;
257 bvec[j] = OCTAVE_C99_ONE;
258 volatile octave_idx_type nm = (nr < nc ? nr : nc);
259 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
260 #if defined(CS_VER) && (CS_VER >= 2)
261 CXSPARSE_ZNAME (_ipvec)
262 (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr);
263 #else
264 CXSPARSE_ZNAME (_ipvec)
265 (nr, S->Pinv, bvec, reinterpret_cast<cs_complex_t *>(buf));
266 #endif
267 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
268 for (volatile octave_idx_type i = 0; i < nm; i++)
269 {
270 OCTAVE_QUIT;
271 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
272 CXSPARSE_ZNAME (_happly)
273 (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
274 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
275 }
276 for (octave_idx_type i = 0; i < nr; i++)
277 vec[i+idx] = buf[i];
278 bvec[j] = OCTAVE_C99_ZERO;
279 }
280 }
281 return ret.hermitian ();
229 #else 282 #else
230 return ComplexMatrix (); 283 return ComplexMatrix ();
231 #endif 284 #endif
232 } 285 }
233 286