Mercurial > hg > octave-lyh

comparison liboctave/SparseQR.cc @ 10314:07ebe522dac2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave C++ sources
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:23:32 -0500
parents 4c0cdbe0acca
children 12884915a8e4
comparison
equal deleted inserted replaced
10313:f3b65e1ae355 10314:07ebe522dac2
183 ret = Matrix (nc, b_nc, 0.0); 183 ret = Matrix (nc, b_nc, 0.0);
184 else 184 else
185 { 185 {
186 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 186 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
187 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) 187 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
188 { 188 {
189 octave_quit (); 189 octave_quit ();
190 for (octave_idx_type i = nr; i < S->m2; i++) 190 for (octave_idx_type i = nr; i < S->m2; i++)
191 buf[i] = 0.; 191 buf[i] = 0.;
192 volatile octave_idx_type nm = (nr < nc ? nr : nc); 192 volatile octave_idx_type nm = (nr < nc ? nr : nc);
193 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 193 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
194 #if defined(CS_VER) && (CS_VER >= 2) 194 #if defined(CS_VER) && (CS_VER >= 2)
195 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr); 195 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
196 #else 196 #else
197 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf); 197 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
198 #endif 198 #endif
199 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 199 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
200 200
201 for (volatile octave_idx_type i = 0; i < nm; i++) 201 for (volatile octave_idx_type i = 0; i < nm; i++)
202 { 202 {
203 octave_quit (); 203 octave_quit ();
204 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 204 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
205 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); 205 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
206 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 206 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
207 } 207 }
208 for (octave_idx_type i = 0; i < b_nr; i++) 208 for (octave_idx_type i = 0; i < b_nr; i++)
209 vec[i+idx] = buf[i]; 209 vec[i+idx] = buf[i];
210 } 210 }
211 } 211 }
212 return ret; 212 return ret;
213 #else 213 #else
214 return Matrix (); 214 return Matrix ();
215 #endif 215 #endif
229 ret = Matrix (nc, nr, 0.0); 229 ret = Matrix (nc, nr, 0.0);
230 else 230 else
231 { 231 {
232 OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1); 232 OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
233 for (octave_idx_type i = 0; i < nr; i++) 233 for (octave_idx_type i = 0; i < nr; i++)
234 bvec[i] = 0.; 234 bvec[i] = 0.;
235 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 235 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
236 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) 236 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
237 { 237 {
238 octave_quit (); 238 octave_quit ();
239 bvec[j] = 1.0; 239 bvec[j] = 1.0;
240 for (octave_idx_type i = nr; i < S->m2; i++) 240 for (octave_idx_type i = nr; i < S->m2; i++)
241 buf[i] = 0.; 241 buf[i] = 0.;
242 volatile octave_idx_type nm = (nr < nc ? nr : nc); 242 volatile octave_idx_type nm = (nr < nc ? nr : nc);
243 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 243 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
244 #if defined(CS_VER) && (CS_VER >= 2) 244 #if defined(CS_VER) && (CS_VER >= 2)
245 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr); 245 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
246 #else 246 #else
247 CXSPARSE_DNAME (_ipvec) (nr, S->Pinv, bvec, buf); 247 CXSPARSE_DNAME (_ipvec) (nr, S->Pinv, bvec, buf);
248 #endif 248 #endif
249 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 249 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
250 250
251 for (volatile octave_idx_type i = 0; i < nm; i++) 251 for (volatile octave_idx_type i = 0; i < nm; i++)
252 { 252 {
253 octave_quit (); 253 octave_quit ();
254 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 254 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
255 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); 255 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
256 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 256 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
257 } 257 }
258 for (octave_idx_type i = 0; i < nr; i++) 258 for (octave_idx_type i = 0; i < nr; i++)
259 vec[i+idx] = buf[i]; 259 vec[i+idx] = buf[i];
260 bvec[j] = 0.0; 260 bvec[j] = 0.0;
261 } 261 }
262 } 262 }
263 return ret.transpose (); 263 return ret.transpose ();
264 #else 264 #else
265 return Matrix (); 265 return Matrix ();
266 #endif 266 #endif
285 x = Matrix (nc, b_nc, 0.0); 285 x = Matrix (nc, b_nc, 0.0);
286 else if (nr >= nc) 286 else if (nr >= nc)
287 { 287 {
288 SparseQR q (a, 3); 288 SparseQR q (a, 3);
289 if (! q.ok ()) 289 if (! q.ok ())
290 return Matrix(); 290 return Matrix();
291 x.resize(nc, b_nc); 291 x.resize(nc, b_nc);
292 double *vec = x.fortran_vec(); 292 double *vec = x.fortran_vec();
293 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 293 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
294 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 294 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
295 i++, idx+=nc, bidx+=b_nr) 295 i++, idx+=nc, bidx+=b_nr)
296 { 296 {
297 octave_quit (); 297 octave_quit ();
298 for (octave_idx_type j = nr; j < q.S()->m2; j++) 298 for (octave_idx_type j = nr; j < q.S()->m2; j++)
299 buf[j] = 0.; 299 buf[j] = 0.;
300 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 300 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
301 #if defined(CS_VER) && (CS_VER >= 2) 301 #if defined(CS_VER) && (CS_VER >= 2)
302 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr); 302 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr);
303 #else 303 #else
304 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); 304 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
305 #endif 305 #endif
306 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 306 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
307 for (volatile octave_idx_type j = 0; j < nc; j++) 307 for (volatile octave_idx_type j = 0; j < nc; j++)
308 { 308 {
309 octave_quit (); 309 octave_quit ();
310 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 310 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
311 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 311 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
312 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 312 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
313 } 313 }
314 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 314 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
315 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 315 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
316 #if defined(CS_VER) && (CS_VER >= 2) 316 #if defined(CS_VER) && (CS_VER >= 2)
317 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, vec + idx, nc); 317 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, vec + idx, nc);
318 #else 318 #else
319 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 319 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
320 #endif 320 #endif
321 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 321 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
322 } 322 }
323 info = 0; 323 info = 0;
324 } 324 }
325 else 325 else
326 { 326 {
327 SparseMatrix at = a.hermitian(); 327 SparseMatrix at = a.hermitian();
328 SparseQR q (at, 3); 328 SparseQR q (at, 3);
329 if (! q.ok ()) 329 if (! q.ok ())
330 return Matrix(); 330 return Matrix();
331 x.resize(nc, b_nc); 331 x.resize(nc, b_nc);
332 double *vec = x.fortran_vec(); 332 double *vec = x.fortran_vec();
333 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 333 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
334 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 334 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
335 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 335 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
336 i++, idx+=nc, bidx+=b_nr) 336 i++, idx+=nc, bidx+=b_nr)
337 { 337 {
338 octave_quit (); 338 octave_quit ();
339 for (octave_idx_type j = nr; j < nbuf; j++) 339 for (octave_idx_type j = nr; j < nbuf; j++)
340 buf[j] = 0.; 340 buf[j] = 0.;
341 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 341 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
342 #if defined(CS_VER) && (CS_VER >= 2) 342 #if defined(CS_VER) && (CS_VER >= 2)
343 CXSPARSE_DNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr); 343 CXSPARSE_DNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr);
344 #else 344 #else
345 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); 345 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf);
346 #endif 346 #endif
347 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 347 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
348 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 348 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
349 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 349 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
350 { 350 {
351 octave_quit (); 351 octave_quit ();
352 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 352 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
353 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 353 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
354 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 354 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
355 } 355 }
356 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 356 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
357 #if defined(CS_VER) && (CS_VER >= 2) 357 #if defined(CS_VER) && (CS_VER >= 2)
358 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc); 358 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc);
359 #else 359 #else
360 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 360 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
361 #endif 361 #endif
362 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 362 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
363 } 363 }
364 info = 0; 364 info = 0;
365 } 365 }
366 366
367 return x; 367 return x;
368 #else 368 #else
389 x = SparseMatrix (nc, b_nc); 389 x = SparseMatrix (nc, b_nc);
390 else if (nr >= nc) 390 else if (nr >= nc)
391 { 391 {
392 SparseQR q (a, 3); 392 SparseQR q (a, 3);
393 if (! q.ok ()) 393 if (! q.ok ())
394 return SparseMatrix(); 394 return SparseMatrix();
395 x = SparseMatrix (nc, b_nc, b.nzmax()); 395 x = SparseMatrix (nc, b_nc, b.nzmax());
396 x.xcidx(0) = 0; 396 x.xcidx(0) = 0;
397 x_nz = b.nzmax(); 397 x_nz = b.nzmax();
398 ii = 0; 398 ii = 0;
399 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 399 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
400 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 400 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
401 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 401 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
402 { 402 {
403 octave_quit (); 403 octave_quit ();
404 for (octave_idx_type j = 0; j < b_nr; j++) 404 for (octave_idx_type j = 0; j < b_nr; j++)
405 Xx[j] = b.xelem(j,i); 405 Xx[j] = b.xelem(j,i);
406 for (octave_idx_type j = nr; j < q.S()->m2; j++) 406 for (octave_idx_type j = nr; j < q.S()->m2; j++)
407 buf[j] = 0.; 407 buf[j] = 0.;
408 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 408 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
409 #if defined(CS_VER) && (CS_VER >= 2) 409 #if defined(CS_VER) && (CS_VER >= 2)
410 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); 410 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr);
411 #else 411 #else
412 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 412 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
413 #endif 413 #endif
414 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 414 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
415 for (volatile octave_idx_type j = 0; j < nc; j++) 415 for (volatile octave_idx_type j = 0; j < nc; j++)
416 { 416 {
417 octave_quit (); 417 octave_quit ();
418 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 418 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
419 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 419 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
420 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 420 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
421 } 421 }
422 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 422 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
423 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 423 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
424 #if defined(CS_VER) && (CS_VER >= 2) 424 #if defined(CS_VER) && (CS_VER >= 2)
425 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); 425 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc);
426 #else 426 #else
427 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 427 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
428 #endif 428 #endif
429 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 429 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
430 430
431 for (octave_idx_type j = 0; j < nc; j++) 431 for (octave_idx_type j = 0; j < nc; j++)
432 { 432 {
433 double tmp = Xx[j]; 433 double tmp = Xx[j];
434 if (tmp != 0.0) 434 if (tmp != 0.0)
435 { 435 {
436 if (ii == x_nz) 436 if (ii == x_nz)
437 { 437 {
438 // Resize the sparse matrix 438 // Resize the sparse matrix
439 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 439 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
440 sz = (sz > 10 ? sz : 10) + x_nz; 440 sz = (sz > 10 ? sz : 10) + x_nz;
441 x.change_capacity (sz); 441 x.change_capacity (sz);
442 x_nz = sz; 442 x_nz = sz;
443 } 443 }
444 x.xdata(ii) = tmp; 444 x.xdata(ii) = tmp;
445 x.xridx(ii++) = j; 445 x.xridx(ii++) = j;
446 } 446 }
447 } 447 }
448 x.xcidx(i+1) = ii; 448 x.xcidx(i+1) = ii;
449 } 449 }
450 info = 0; 450 info = 0;
451 } 451 }
452 else 452 else
453 { 453 {
454 SparseMatrix at = a.hermitian(); 454 SparseMatrix at = a.hermitian();
455 SparseQR q (at, 3); 455 SparseQR q (at, 3);
456 if (! q.ok ()) 456 if (! q.ok ())
457 return SparseMatrix(); 457 return SparseMatrix();
458 x = SparseMatrix (nc, b_nc, b.nzmax()); 458 x = SparseMatrix (nc, b_nc, b.nzmax());
459 x.xcidx(0) = 0; 459 x.xcidx(0) = 0;
460 x_nz = b.nzmax(); 460 x_nz = b.nzmax();
461 ii = 0; 461 ii = 0;
462 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 462 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
463 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 463 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
464 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 464 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
465 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 465 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
466 { 466 {
467 octave_quit (); 467 octave_quit ();
468 for (octave_idx_type j = 0; j < b_nr; j++) 468 for (octave_idx_type j = 0; j < b_nr; j++)
469 Xx[j] = b.xelem(j,i); 469 Xx[j] = b.xelem(j,i);
470 for (octave_idx_type j = nr; j < nbuf; j++) 470 for (octave_idx_type j = nr; j < nbuf; j++)
471 buf[j] = 0.; 471 buf[j] = 0.;
472 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 472 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
473 #if defined(CS_VER) && (CS_VER >= 2) 473 #if defined(CS_VER) && (CS_VER >= 2)
474 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); 474 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr);
475 #else 475 #else
476 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 476 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
477 #endif 477 #endif
478 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 478 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
479 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 479 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
480 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 480 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
481 { 481 {
482 octave_quit (); 482 octave_quit ();
483 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 483 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
484 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 484 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
485 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 485 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
486 } 486 }
487 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 487 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
488 #if defined(CS_VER) && (CS_VER >= 2) 488 #if defined(CS_VER) && (CS_VER >= 2)
489 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); 489 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc);
490 #else 490 #else
491 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 491 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
492 #endif 492 #endif
493 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 493 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
494 494
495 for (octave_idx_type j = 0; j < nc; j++) 495 for (octave_idx_type j = 0; j < nc; j++)
496 { 496 {
497 double tmp = Xx[j]; 497 double tmp = Xx[j];
498 if (tmp != 0.0) 498 if (tmp != 0.0)
499 { 499 {
500 if (ii == x_nz) 500 if (ii == x_nz)
501 { 501 {
502 // Resize the sparse matrix 502 // Resize the sparse matrix
503 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 503 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
504 sz = (sz > 10 ? sz : 10) + x_nz; 504 sz = (sz > 10 ? sz : 10) + x_nz;
505 x.change_capacity (sz); 505 x.change_capacity (sz);
506 x_nz = sz; 506 x_nz = sz;
507 } 507 }
508 x.xdata(ii) = tmp; 508 x.xdata(ii) = tmp;
509 x.xridx(ii++) = j; 509 x.xridx(ii++) = j;
510 } 510 }
511 } 511 }
512 x.xcidx(i+1) = ii; 512 x.xcidx(i+1) = ii;
513 } 513 }
514 info = 0; 514 info = 0;
515 } 515 }
516 516
517 x.maybe_compress (); 517 x.maybe_compress ();
518 return x; 518 return x;
539 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); 539 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
540 else if (nr >= nc) 540 else if (nr >= nc)
541 { 541 {
542 SparseQR q (a, 3); 542 SparseQR q (a, 3);
543 if (! q.ok ()) 543 if (! q.ok ())
544 return ComplexMatrix(); 544 return ComplexMatrix();
545 x.resize(nc, b_nc); 545 x.resize(nc, b_nc);
546 Complex *vec = x.fortran_vec(); 546 Complex *vec = x.fortran_vec();
547 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 547 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
548 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 548 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
549 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 549 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
550 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 550 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
551 { 551 {
552 octave_quit (); 552 octave_quit ();
553 for (octave_idx_type j = 0; j < b_nr; j++) 553 for (octave_idx_type j = 0; j < b_nr; j++)
554 { 554 {
555 Complex c = b.xelem (j,i); 555 Complex c = b.xelem (j,i);
556 Xx[j] = std::real (c); 556 Xx[j] = std::real (c);
557 Xz[j] = std::imag (c); 557 Xz[j] = std::imag (c);
558 } 558 }
559 for (octave_idx_type j = nr; j < q.S()->m2; j++) 559 for (octave_idx_type j = nr; j < q.S()->m2; j++)
560 buf[j] = 0.; 560 buf[j] = 0.;
561 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 561 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
562 #if defined(CS_VER) && (CS_VER >= 2) 562 #if defined(CS_VER) && (CS_VER >= 2)
563 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); 563 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr);
564 #else 564 #else
565 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 565 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
566 #endif 566 #endif
567 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 567 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
568 for (volatile octave_idx_type j = 0; j < nc; j++) 568 for (volatile octave_idx_type j = 0; j < nc; j++)
569 { 569 {
570 octave_quit (); 570 octave_quit ();
571 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 571 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
572 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 572 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
573 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 573 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
574 } 574 }
575 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 575 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
576 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 576 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
577 #if defined(CS_VER) && (CS_VER >= 2) 577 #if defined(CS_VER) && (CS_VER >= 2)
578 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); 578 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc);
579 #else 579 #else
580 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 580 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
581 #endif 581 #endif
582 for (octave_idx_type j = nr; j < q.S()->m2; j++) 582 for (octave_idx_type j = nr; j < q.S()->m2; j++)
583 buf[j] = 0.; 583 buf[j] = 0.;
584 #if defined(CS_VER) && (CS_VER >= 2) 584 #if defined(CS_VER) && (CS_VER >= 2)
585 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); 585 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr);
586 #else 586 #else
587 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); 587 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
588 #endif 588 #endif
589 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 589 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
590 for (volatile octave_idx_type j = 0; j < nc; j++) 590 for (volatile octave_idx_type j = 0; j < nc; j++)
591 { 591 {
592 octave_quit (); 592 octave_quit ();
593 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 593 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
594 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 594 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
595 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 595 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
596 } 596 }
597 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 597 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
598 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 598 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
599 #if defined(CS_VER) && (CS_VER >= 2) 599 #if defined(CS_VER) && (CS_VER >= 2)
600 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); 600 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc);
601 #else 601 #else
602 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); 602 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz);
603 #endif 603 #endif
604 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 604 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
605 for (octave_idx_type j = 0; j < nc; j++) 605 for (octave_idx_type j = 0; j < nc; j++)
606 vec[j+idx] = Complex (Xx[j], Xz[j]); 606 vec[j+idx] = Complex (Xx[j], Xz[j]);
607 } 607 }
608 info = 0; 608 info = 0;
609 } 609 }
610 else 610 else
611 { 611 {
612 SparseMatrix at = a.hermitian(); 612 SparseMatrix at = a.hermitian();
613 SparseQR q (at, 3); 613 SparseQR q (at, 3);
614 if (! q.ok ()) 614 if (! q.ok ())
615 return ComplexMatrix(); 615 return ComplexMatrix();
616 x.resize(nc, b_nc); 616 x.resize(nc, b_nc);
617 Complex *vec = x.fortran_vec(); 617 Complex *vec = x.fortran_vec();
618 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 618 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
619 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 619 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
620 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 620 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
621 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 621 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
622 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 622 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
623 { 623 {
624 octave_quit (); 624 octave_quit ();
625 for (octave_idx_type j = 0; j < b_nr; j++) 625 for (octave_idx_type j = 0; j < b_nr; j++)
626 { 626 {
627 Complex c = b.xelem (j,i); 627 Complex c = b.xelem (j,i);
628 Xx[j] = std::real (c); 628 Xx[j] = std::real (c);
629 Xz[j] = std::imag (c); 629 Xz[j] = std::imag (c);
630 } 630 }
631 for (octave_idx_type j = nr; j < nbuf; j++) 631 for (octave_idx_type j = nr; j < nbuf; j++)
632 buf[j] = 0.; 632 buf[j] = 0.;
633 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 633 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
634 #if defined(CS_VER) && (CS_VER >= 2) 634 #if defined(CS_VER) && (CS_VER >= 2)
635 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); 635 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr);
636 #else 636 #else
637 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 637 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
638 #endif 638 #endif
639 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 639 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
640 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 640 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
641 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 641 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
642 { 642 {
643 octave_quit (); 643 octave_quit ();
644 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 644 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
645 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 645 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
646 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 646 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
647 } 647 }
648 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 648 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
649 #if defined(CS_VER) && (CS_VER >= 2) 649 #if defined(CS_VER) && (CS_VER >= 2)
650 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); 650 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc);
651 #else 651 #else
652 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 652 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
653 #endif 653 #endif
654 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 654 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
655 for (octave_idx_type j = nr; j < nbuf; j++) 655 for (octave_idx_type j = nr; j < nbuf; j++)
656 buf[j] = 0.; 656 buf[j] = 0.;
657 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 657 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
658 #if defined(CS_VER) && (CS_VER >= 2) 658 #if defined(CS_VER) && (CS_VER >= 2)
659 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); 659 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr);
660 #else 660 #else
661 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); 661 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
662 #endif 662 #endif
663 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 663 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
664 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 664 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
665 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 665 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
666 { 666 {
667 octave_quit (); 667 octave_quit ();
668 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 668 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
669 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 669 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
670 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 670 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
671 } 671 }
672 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 672 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
673 #if defined(CS_VER) && (CS_VER >= 2) 673 #if defined(CS_VER) && (CS_VER >= 2)
674 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); 674 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc);
675 #else 675 #else
676 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); 676 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz);
677 #endif 677 #endif
678 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 678 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
679 for (octave_idx_type j = 0; j < nc; j++) 679 for (octave_idx_type j = 0; j < nc; j++)
680 vec[j+idx] = Complex (Xx[j], Xz[j]); 680 vec[j+idx] = Complex (Xx[j], Xz[j]);
681 } 681 }
682 info = 0; 682 info = 0;
683 } 683 }
684 684
685 return x; 685 return x;
686 #else 686 #else
707 x = SparseComplexMatrix (nc, b_nc); 707 x = SparseComplexMatrix (nc, b_nc);
708 else if (nr >= nc) 708 else if (nr >= nc)
709 { 709 {
710 SparseQR q (a, 3); 710 SparseQR q (a, 3);
711 if (! q.ok ()) 711 if (! q.ok ())
712 return SparseComplexMatrix(); 712 return SparseComplexMatrix();
713 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 713 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
714 x.xcidx(0) = 0; 714 x.xcidx(0) = 0;
715 x_nz = b.nzmax(); 715 x_nz = b.nzmax();
716 ii = 0; 716 ii = 0;
717 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 717 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
718 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 718 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
719 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 719 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
720 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 720 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
721 { 721 {
722 octave_quit (); 722 octave_quit ();
723 for (octave_idx_type j = 0; j < b_nr; j++) 723 for (octave_idx_type j = 0; j < b_nr; j++)
724 { 724 {
725 Complex c = b.xelem (j,i); 725 Complex c = b.xelem (j,i);
726 Xx[j] = std::real (c); 726 Xx[j] = std::real (c);
727 Xz[j] = std::imag (c); 727 Xz[j] = std::imag (c);
728 } 728 }
729 for (octave_idx_type j = nr; j < q.S()->m2; j++) 729 for (octave_idx_type j = nr; j < q.S()->m2; j++)
730 buf[j] = 0.; 730 buf[j] = 0.;
731 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 731 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
732 #if defined(CS_VER) && (CS_VER >= 2) 732 #if defined(CS_VER) && (CS_VER >= 2)
733 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); 733 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr);
734 #else 734 #else
735 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 735 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
736 #endif 736 #endif
737 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 737 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
738 for (volatile octave_idx_type j = 0; j < nc; j++) 738 for (volatile octave_idx_type j = 0; j < nc; j++)
739 { 739 {
740 octave_quit (); 740 octave_quit ();
741 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 741 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
742 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 742 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
743 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 743 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
744 } 744 }
745 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 745 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
746 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 746 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
747 #if defined(CS_VER) && (CS_VER >= 2) 747 #if defined(CS_VER) && (CS_VER >= 2)
748 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); 748 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc);
749 #else 749 #else
750 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 750 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
751 #endif 751 #endif
752 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 752 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
753 for (octave_idx_type j = nr; j < q.S()->m2; j++) 753 for (octave_idx_type j = nr; j < q.S()->m2; j++)
754 buf[j] = 0.; 754 buf[j] = 0.;
755 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 755 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
756 #if defined(CS_VER) && (CS_VER >= 2) 756 #if defined(CS_VER) && (CS_VER >= 2)
757 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); 757 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr);
758 #else 758 #else
759 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); 759 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
760 #endif 760 #endif
761 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 761 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
762 for (volatile octave_idx_type j = 0; j < nc; j++) 762 for (volatile octave_idx_type j = 0; j < nc; j++)
763 { 763 {
764 octave_quit (); 764 octave_quit ();
765 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 765 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
766 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 766 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
767 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 767 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
768 } 768 }
769 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 769 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
770 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 770 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
771 #if defined(CS_VER) && (CS_VER >= 2) 771 #if defined(CS_VER) && (CS_VER >= 2)
772 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); 772 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc);
773 #else 773 #else
774 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); 774 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz);
775 #endif 775 #endif
776 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 776 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
777 777
778 for (octave_idx_type j = 0; j < nc; j++) 778 for (octave_idx_type j = 0; j < nc; j++)
779 { 779 {
780 Complex tmp = Complex (Xx[j], Xz[j]); 780 Complex tmp = Complex (Xx[j], Xz[j]);
781 if (tmp != 0.0) 781 if (tmp != 0.0)
782 { 782 {
783 if (ii == x_nz) 783 if (ii == x_nz)
784 { 784 {
785 // Resize the sparse matrix 785 // Resize the sparse matrix
786 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 786 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
787 sz = (sz > 10 ? sz : 10) + x_nz; 787 sz = (sz > 10 ? sz : 10) + x_nz;
788 x.change_capacity (sz); 788 x.change_capacity (sz);
789 x_nz = sz; 789 x_nz = sz;
790 } 790 }
791 x.xdata(ii) = tmp; 791 x.xdata(ii) = tmp;
792 x.xridx(ii++) = j; 792 x.xridx(ii++) = j;
793 } 793 }
794 } 794 }
795 x.xcidx(i+1) = ii; 795 x.xcidx(i+1) = ii;
796 } 796 }
797 info = 0; 797 info = 0;
798 } 798 }
799 else 799 else
800 { 800 {
801 SparseMatrix at = a.hermitian(); 801 SparseMatrix at = a.hermitian();
802 SparseQR q (at, 3); 802 SparseQR q (at, 3);
803 if (! q.ok ()) 803 if (! q.ok ())
804 return SparseComplexMatrix(); 804 return SparseComplexMatrix();
805 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 805 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
806 x.xcidx(0) = 0; 806 x.xcidx(0) = 0;
807 x_nz = b.nzmax(); 807 x_nz = b.nzmax();
808 ii = 0; 808 ii = 0;
809 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 809 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
810 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 810 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
811 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 811 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
812 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 812 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
813 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 813 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
814 { 814 {
815 octave_quit (); 815 octave_quit ();
816 for (octave_idx_type j = 0; j < b_nr; j++) 816 for (octave_idx_type j = 0; j < b_nr; j++)
817 { 817 {
818 Complex c = b.xelem (j,i); 818 Complex c = b.xelem (j,i);
819 Xx[j] = std::real (c); 819 Xx[j] = std::real (c);
820 Xz[j] = std::imag (c); 820 Xz[j] = std::imag (c);
821 } 821 }
822 for (octave_idx_type j = nr; j < nbuf; j++) 822 for (octave_idx_type j = nr; j < nbuf; j++)
823 buf[j] = 0.; 823 buf[j] = 0.;
824 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 824 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
825 #if defined(CS_VER) && (CS_VER >= 2) 825 #if defined(CS_VER) && (CS_VER >= 2)
826 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); 826 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr);
827 #else 827 #else
828 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 828 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
829 #endif 829 #endif
830 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 830 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
831 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 831 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
832 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 832 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
833 { 833 {
834 octave_quit (); 834 octave_quit ();
835 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 835 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
836 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 836 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
837 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 837 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
838 } 838 }
839 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 839 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
840 #if defined(CS_VER) && (CS_VER >= 2) 840 #if defined(CS_VER) && (CS_VER >= 2)
841 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); 841 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc);
842 #else 842 #else
843 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 843 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
844 #endif 844 #endif
845 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 845 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
846 for (octave_idx_type j = nr; j < nbuf; j++) 846 for (octave_idx_type j = nr; j < nbuf; j++)
847 buf[j] = 0.; 847 buf[j] = 0.;
848 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 848 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
849 #if defined(CS_VER) && (CS_VER >= 2) 849 #if defined(CS_VER) && (CS_VER >= 2)
850 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); 850 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr);
851 #else 851 #else
852 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); 852 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
853 #endif 853 #endif
854 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 854 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
855 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 855 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
856 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 856 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
857 { 857 {
858 octave_quit (); 858 octave_quit ();
859 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 859 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
860 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 860 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
861 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 861 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
862 } 862 }
863 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 863 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
864 #if defined(CS_VER) && (CS_VER >= 2) 864 #if defined(CS_VER) && (CS_VER >= 2)
865 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); 865 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc);
866 #else 866 #else
867 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); 867 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz);
868 #endif 868 #endif
869 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 869 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
870 870
871 for (octave_idx_type j = 0; j < nc; j++) 871 for (octave_idx_type j = 0; j < nc; j++)
872 { 872 {
873 Complex tmp = Complex (Xx[j], Xz[j]); 873 Complex tmp = Complex (Xx[j], Xz[j]);
874 if (tmp != 0.0) 874 if (tmp != 0.0)
875 { 875 {
876 if (ii == x_nz) 876 if (ii == x_nz)
877 { 877 {
878 // Resize the sparse matrix 878 // Resize the sparse matrix
879 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 879 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
880 sz = (sz > 10 ? sz : 10) + x_nz; 880 sz = (sz > 10 ? sz : 10) + x_nz;
881 x.change_capacity (sz); 881 x.change_capacity (sz);
882 x_nz = sz; 882 x_nz = sz;
883 } 883 }
884 x.xdata(ii) = tmp; 884 x.xdata(ii) = tmp;
885 x.xridx(ii++) = j; 885 x.xridx(ii++) = j;
886 } 886 }
887 } 887 }
888 x.xcidx(i+1) = ii; 888 x.xcidx(i+1) = ii;
889 } 889 }
890 info = 0; 890 info = 0;
891 } 891 }
892 892
893 x.maybe_compress (); 893 x.maybe_compress ();
894 return x; 894 return x;
897 #endif 897 #endif
898 } 898 }
899 899
900 Matrix 900 Matrix
901 qrsolve(const SparseMatrix &a, const MArray2<double> &b, 901 qrsolve(const SparseMatrix &a, const MArray2<double> &b,
902 octave_idx_type &info) 902 octave_idx_type &info)
903 { 903 {
904 return qrsolve (a, Matrix (b), info); 904 return qrsolve (a, Matrix (b), info);
905 } 905 }
906 906
907 ComplexMatrix 907 ComplexMatrix
908 qrsolve(const SparseMatrix &a, const MArray2<Complex> &b, 908 qrsolve(const SparseMatrix &a, const MArray2<Complex> &b,
909 octave_idx_type &info) 909 octave_idx_type &info)
910 { 910 {
911 return qrsolve (a, ComplexMatrix (b), info); 911 return qrsolve (a, ComplexMatrix (b), info);
912 } 912 }