Mercurial > hg > octave-lyh

comparison liboctave/dSparse.cc @ 15018:3d8ace26c5b4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding conventions for cuddled parentheses in liboctave/. * Array-util.cc, Array.cc, Array.h, CMatrix.cc, CNDArray.cc, CSparse.cc, CmplxQR.cc, CollocWt.cc, DASPK.cc, DASRT.cc, DASSL.cc, EIG.cc, LSODE.cc, MSparse.cc, MatrixType.cc, Sparse-op-defs.h, Sparse-perm-op-defs.h, Sparse.cc, Sparse.h, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, boolSparse.cc, cmd-hist.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dSparse.cc, data-conv.cc, dbleQR.cc, dbleSVD.cc, dim-vector.cc, eigs-base.cc, f2c-main.c, fCMatrix.cc, fCNDArray.cc, fCmplxQR.cc, fEIG.cc, fMatrix.cc, fNDArray.cc, floatQR.cc, floatSVD.cc, idx-vector.cc, kpse.cc, lo-specfun.cc, mx-inlines.cc, mx-op-defs.h, oct-alloc.cc, oct-binmap.h, oct-fftw.cc, oct-group.h, oct-inttypes.cc, oct-inttypes.h, oct-locbuf.cc, oct-md5.cc, oct-rand.cc, oct-sort.cc, oct-syscalls.cc, randgamma.c, randmtzig.c, randpoisson.c, sparse-base-chol.cc, sparse-base-lu.cc, sparse-dmsolve.cc, str-vec.cc, str-vec.h, tempnam.c, tempname.c: Use Octave coding conventions for cuddled parentheses in liboctave/.
author Rik <rik@octave.org>
date Thu, 26 Jul 2012 08:13:22 -0700
parents 460a3c6d8bf1
children 4bbd3bbb8912
comparison
equal deleted inserted replaced
15017:dd4ad69e4ab9 15018:3d8ace26c5b4
175 ridx (j) = i; 175 ridx (j) = i;
176 j++; 176 j++;
177 } 177 }
178 } 178 }
179 for (octave_idx_type i = l; i <= a.cols (); i++) 179 for (octave_idx_type i = l; i <= a.cols (); i++)
180 cidx(i) = j; 180 cidx (i) = j;
181 } 181 }
182 182
183 bool 183 bool
184 SparseMatrix::operator == (const SparseMatrix& a) const 184 SparseMatrix::operator == (const SparseMatrix& a) const
185 { 185 {
192 192
193 if (nr != nr_a || nc != nc_a || nz != nz_a) 193 if (nr != nr_a || nc != nc_a || nz != nz_a)
194 return false; 194 return false;
195 195
196 for (octave_idx_type i = 0; i < nc + 1; i++) 196 for (octave_idx_type i = 0; i < nc + 1; i++)
197 if (cidx(i) != a.cidx(i)) 197 if (cidx (i) != a.cidx (i))
198 return false; 198 return false;
199 199
200 for (octave_idx_type i = 0; i < nz; i++) 200 for (octave_idx_type i = 0; i < nz; i++)
201 if (data(i) != a.data(i) || ridx(i) != a.ridx(i)) 201 if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
202 return false; 202 return false;
203 203
204 return true; 204 return true;
205 } 205 }
206 206
218 218
219 if (nr == nc && nr > 0) 219 if (nr == nc && nr > 0)
220 { 220 {
221 for (octave_idx_type j = 0; j < nc; j++) 221 for (octave_idx_type j = 0; j < nc; j++)
222 { 222 {
223 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 223 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
224 { 224 {
225 octave_idx_type ri = ridx(i); 225 octave_idx_type ri = ridx (i);
226 226
227 if (ri != j) 227 if (ri != j)
228 { 228 {
229 bool found = false; 229 bool found = false;
230 230
231 for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++) 231 for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
232 { 232 {
233 if (ridx(k) == j) 233 if (ridx (k) == j)
234 { 234 {
235 if (data(i) == data(k)) 235 if (data (i) == data (k))
236 found = true; 236 found = true;
237 break; 237 break;
238 } 238 }
239 } 239 }
240 240
292 octave_idx_type nel = 0; 292 octave_idx_type nel = 0;
293 for (octave_idx_type j = 0; j < nc; j++) 293 for (octave_idx_type j = 0; j < nc; j++)
294 { 294 {
295 double tmp_max = octave_NaN; 295 double tmp_max = octave_NaN;
296 octave_idx_type idx_j = 0; 296 octave_idx_type idx_j = 0;
297 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 297 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
298 { 298 {
299 if (ridx(i) != idx_j) 299 if (ridx (i) != idx_j)
300 break; 300 break;
301 else 301 else
302 idx_j++; 302 idx_j++;
303 } 303 }
304 304
305 if (idx_j != nr) 305 if (idx_j != nr)
306 tmp_max = 0.; 306 tmp_max = 0.;
307 307
308 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 308 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
309 { 309 {
310 double tmp = data (i); 310 double tmp = data (i);
311 311
312 if (xisnan (tmp)) 312 if (xisnan (tmp))
313 continue; 313 continue;
342 } 342 }
343 else 343 else
344 { 344 {
345 idx_arg.resize (dim_vector (nr, 1), 0); 345 idx_arg.resize (dim_vector (nr, 1), 0);
346 346
347 for (octave_idx_type i = cidx(0); i < cidx(1); i++) 347 for (octave_idx_type i = cidx (0); i < cidx (1); i++)
348 idx_arg.elem(ridx(i)) = -1; 348 idx_arg.elem (ridx (i)) = -1;
349 349
350 for (octave_idx_type j = 0; j < nc; j++) 350 for (octave_idx_type j = 0; j < nc; j++)
351 for (octave_idx_type i = 0; i < nr; i++) 351 for (octave_idx_type i = 0; i < nr; i++)
352 { 352 {
353 if (idx_arg.elem(i) != -1) 353 if (idx_arg.elem (i) != -1)
354 continue; 354 continue;
355 bool found = false; 355 bool found = false;
356 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) 356 for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
357 if (ridx(k) == i) 357 if (ridx (k) == i)
358 { 358 {
359 found = true; 359 found = true;
360 break; 360 break;
361 } 361 }
362 362
363 if (!found) 363 if (!found)
364 idx_arg.elem(i) = j; 364 idx_arg.elem (i) = j;
365 365
366 } 366 }
367 367
368 for (octave_idx_type j = 0; j < nc; j++) 368 for (octave_idx_type j = 0; j < nc; j++)
369 { 369 {
370 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 370 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
371 { 371 {
372 octave_idx_type ir = ridx (i); 372 octave_idx_type ir = ridx (i);
373 octave_idx_type ix = idx_arg.elem (ir); 373 octave_idx_type ix = idx_arg.elem (ir);
374 double tmp = data (i); 374 double tmp = data (i);
375 375
380 } 380 }
381 } 381 }
382 382
383 octave_idx_type nel = 0; 383 octave_idx_type nel = 0;
384 for (octave_idx_type j = 0; j < nr; j++) 384 for (octave_idx_type j = 0; j < nr; j++)
385 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) 385 if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
386 nel++; 386 nel++;
387 387
388 result = SparseMatrix (nr, 1, nel); 388 result = SparseMatrix (nr, 1, nel);
389 389
390 octave_idx_type ii = 0; 390 octave_idx_type ii = 0;
441 octave_idx_type nel = 0; 441 octave_idx_type nel = 0;
442 for (octave_idx_type j = 0; j < nc; j++) 442 for (octave_idx_type j = 0; j < nc; j++)
443 { 443 {
444 double tmp_min = octave_NaN; 444 double tmp_min = octave_NaN;
445 octave_idx_type idx_j = 0; 445 octave_idx_type idx_j = 0;
446 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 446 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
447 { 447 {
448 if (ridx(i) != idx_j) 448 if (ridx (i) != idx_j)
449 break; 449 break;
450 else 450 else
451 idx_j++; 451 idx_j++;
452 } 452 }
453 453
454 if (idx_j != nr) 454 if (idx_j != nr)
455 tmp_min = 0.; 455 tmp_min = 0.;
456 456
457 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 457 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
458 { 458 {
459 double tmp = data (i); 459 double tmp = data (i);
460 460
461 if (xisnan (tmp)) 461 if (xisnan (tmp))
462 continue; 462 continue;
491 } 491 }
492 else 492 else
493 { 493 {
494 idx_arg.resize (dim_vector (nr, 1), 0); 494 idx_arg.resize (dim_vector (nr, 1), 0);
495 495
496 for (octave_idx_type i = cidx(0); i < cidx(1); i++) 496 for (octave_idx_type i = cidx (0); i < cidx (1); i++)
497 idx_arg.elem(ridx(i)) = -1; 497 idx_arg.elem (ridx (i)) = -1;
498 498
499 for (octave_idx_type j = 0; j < nc; j++) 499 for (octave_idx_type j = 0; j < nc; j++)
500 for (octave_idx_type i = 0; i < nr; i++) 500 for (octave_idx_type i = 0; i < nr; i++)
501 { 501 {
502 if (idx_arg.elem(i) != -1) 502 if (idx_arg.elem (i) != -1)
503 continue; 503 continue;
504 bool found = false; 504 bool found = false;
505 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) 505 for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
506 if (ridx(k) == i) 506 if (ridx (k) == i)
507 { 507 {
508 found = true; 508 found = true;
509 break; 509 break;
510 } 510 }
511 511
512 if (!found) 512 if (!found)
513 idx_arg.elem(i) = j; 513 idx_arg.elem (i) = j;
514 514
515 } 515 }
516 516
517 for (octave_idx_type j = 0; j < nc; j++) 517 for (octave_idx_type j = 0; j < nc; j++)
518 { 518 {
519 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 519 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
520 { 520 {
521 octave_idx_type ir = ridx (i); 521 octave_idx_type ir = ridx (i);
522 octave_idx_type ix = idx_arg.elem (ir); 522 octave_idx_type ix = idx_arg.elem (ir);
523 double tmp = data (i); 523 double tmp = data (i);
524 524
529 } 529 }
530 } 530 }
531 531
532 octave_idx_type nel = 0; 532 octave_idx_type nel = 0;
533 for (octave_idx_type j = 0; j < nr; j++) 533 for (octave_idx_type j = 0; j < nr; j++)
534 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) 534 if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
535 nel++; 535 nel++;
536 536
537 result = SparseMatrix (nr, 1, nel); 537 result = SparseMatrix (nr, 1, nel);
538 538
539 octave_idx_type ii = 0; 539 octave_idx_type ii = 0;
618 octave_idx_type nc = a.cols (); 618 octave_idx_type nc = a.cols ();
619 octave_idx_type nz = a.nnz (); 619 octave_idx_type nz = a.nnz ();
620 SparseMatrix r (nr, nc, nz); 620 SparseMatrix r (nr, nc, nz);
621 621
622 for (octave_idx_type i = 0; i < nc +1; i++) 622 for (octave_idx_type i = 0; i < nc +1; i++)
623 r.cidx(i) = a.cidx(i); 623 r.cidx (i) = a.cidx (i);
624 624
625 for (octave_idx_type i = 0; i < nz; i++) 625 for (octave_idx_type i = 0; i < nz; i++)
626 { 626 {
627 r.data(i) = std::real (a.data(i)); 627 r.data (i) = std::real (a.data (i));
628 r.ridx(i) = a.ridx(i); 628 r.ridx (i) = a.ridx (i);
629 } 629 }
630 630
631 return r; 631 return r;
632 } 632 }
633 633
638 octave_idx_type nc = a.cols (); 638 octave_idx_type nc = a.cols ();
639 octave_idx_type nz = a.nnz (); 639 octave_idx_type nz = a.nnz ();
640 SparseMatrix r (nr, nc, nz); 640 SparseMatrix r (nr, nc, nz);
641 641
642 for (octave_idx_type i = 0; i < nc +1; i++) 642 for (octave_idx_type i = 0; i < nc +1; i++)
643 r.cidx(i) = a.cidx(i); 643 r.cidx (i) = a.cidx (i);
644 644
645 for (octave_idx_type i = 0; i < nz; i++) 645 for (octave_idx_type i = 0; i < nz; i++)
646 { 646 {
647 r.data(i) = std::imag (a.data(i)); 647 r.data (i) = std::imag (a.data (i));
648 r.ridx(i) = a.ridx(i); 648 r.ridx (i) = a.ridx (i);
649 } 649 }
650 650
651 return r; 651 return r;
652 } 652 }
653 653
665 // best way to handle it. 665 // best way to handle it.
666 Matrix tmp (nr, nc, atan2 (x, 0.)); 666 Matrix tmp (nr, nc, atan2 (x, 0.));
667 667
668 for (octave_idx_type j = 0; j < nc; j++) 668 for (octave_idx_type j = 0; j < nc; j++)
669 for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++) 669 for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++)
670 tmp.elem (y.ridx(i), j) = atan2 (x, y.data(i)); 670 tmp.elem (y.ridx (i), j) = atan2 (x, y.data (i));
671 671
672 return SparseMatrix (tmp); 672 return SparseMatrix (tmp);
673 } 673 }
674 } 674 }
675 675
681 octave_idx_type nz = x.nnz (); 681 octave_idx_type nz = x.nnz ();
682 682
683 SparseMatrix retval (nr, nc, nz); 683 SparseMatrix retval (nr, nc, nz);
684 684
685 octave_idx_type ii = 0; 685 octave_idx_type ii = 0;
686 retval.xcidx(0) = 0; 686 retval.xcidx (0) = 0;
687 for (octave_idx_type i = 0; i < nc; i++) 687 for (octave_idx_type i = 0; i < nc; i++)
688 { 688 {
689 for (octave_idx_type j = x.cidx(i); j < x.cidx(i+1); j++) 689 for (octave_idx_type j = x.cidx (i); j < x.cidx (i+1); j++)
690 { 690 {
691 double tmp = atan2 (x.data(j), y); 691 double tmp = atan2 (x.data (j), y);
692 if (tmp != 0.) 692 if (tmp != 0.)
693 { 693 {
694 retval.xdata (ii) = tmp; 694 retval.xdata (ii) = tmp;
695 retval.xridx (ii++) = x.ridx (j); 695 retval.xridx (ii++) = x.ridx (j);
696 } 696 }
735 735
736 octave_idx_type jx = 0; 736 octave_idx_type jx = 0;
737 r.cidx (0) = 0; 737 r.cidx (0) = 0;
738 for (octave_idx_type i = 0 ; i < x_nc ; i++) 738 for (octave_idx_type i = 0 ; i < x_nc ; i++)
739 { 739 {
740 octave_idx_type ja = x.cidx(i); 740 octave_idx_type ja = x.cidx (i);
741 octave_idx_type ja_max = x.cidx(i+1); 741 octave_idx_type ja_max = x.cidx (i+1);
742 bool ja_lt_max= ja < ja_max; 742 bool ja_lt_max= ja < ja_max;
743 743
744 octave_idx_type jb = y.cidx(i); 744 octave_idx_type jb = y.cidx (i);
745 octave_idx_type jb_max = y.cidx(i+1); 745 octave_idx_type jb_max = y.cidx (i+1);
746 bool jb_lt_max = jb < jb_max; 746 bool jb_lt_max = jb < jb_max;
747 747
748 while (ja_lt_max || jb_lt_max ) 748 while (ja_lt_max || jb_lt_max )
749 { 749 {
750 octave_quit (); 750 octave_quit ();
751 if ((! jb_lt_max) || 751 if ((! jb_lt_max) ||
752 (ja_lt_max && (x.ridx(ja) < y.ridx(jb)))) 752 (ja_lt_max && (x.ridx (ja) < y.ridx (jb))))
753 { 753 {
754 r.ridx(jx) = x.ridx(ja); 754 r.ridx (jx) = x.ridx (ja);
755 r.data(jx) = atan2 (x.data(ja), 0.); 755 r.data (jx) = atan2 (x.data (ja), 0.);
756 jx++; 756 jx++;
757 ja++; 757 ja++;
758 ja_lt_max= ja < ja_max; 758 ja_lt_max= ja < ja_max;
759 } 759 }
760 else if (( !ja_lt_max ) || 760 else if (( !ja_lt_max ) ||
761 (jb_lt_max && (y.ridx(jb) < x.ridx(ja)) ) ) 761 (jb_lt_max && (y.ridx (jb) < x.ridx (ja)) ) )
762 { 762 {
763 jb++; 763 jb++;
764 jb_lt_max= jb < jb_max; 764 jb_lt_max= jb < jb_max;
765 } 765 }
766 else 766 else
767 { 767 {
768 double tmp = atan2 (x.data(ja), y.data(jb)); 768 double tmp = atan2 (x.data (ja), y.data (jb));
769 if (tmp != 0.) 769 if (tmp != 0.)
770 { 770 {
771 r.data(jx) = tmp; 771 r.data (jx) = tmp;
772 r.ridx(jx) = x.ridx(ja); 772 r.ridx (jx) = x.ridx (ja);
773 jx++; 773 jx++;
774 } 774 }
775 ja++; 775 ja++;
776 ja_lt_max= ja < ja_max; 776 ja_lt_max= ja < ja_max;
777 jb++; 777 jb++;
778 jb_lt_max= jb < jb_max; 778 jb_lt_max= jb < jb_max;
779 } 779 }
780 } 780 }
781 r.cidx(i+1) = jx; 781 r.cidx (i+1) = jx;
782 } 782 }
783 783
784 r.maybe_compress (); 784 r.maybe_compress ();
785 } 785 }
786 } 786 }
847 if (calccond) 847 if (calccond)
848 { 848 {
849 double dmax = 0., dmin = octave_Inf; 849 double dmax = 0., dmin = octave_Inf;
850 for (octave_idx_type i = 0; i < nr; i++) 850 for (octave_idx_type i = 0; i < nr; i++)
851 { 851 {
852 double tmp = fabs(v[i]); 852 double tmp = fabs (v[i]);
853 if (tmp > dmax) 853 if (tmp > dmax)
854 dmax = tmp; 854 dmax = tmp;
855 if (tmp < dmin) 855 if (tmp < dmin)
856 dmin = tmp; 856 dmin = tmp;
857 } 857 }
897 { 897 {
898 // Calculate the 1-norm of matrix for rcond calculation 898 // Calculate the 1-norm of matrix for rcond calculation
899 for (octave_idx_type j = 0; j < nr; j++) 899 for (octave_idx_type j = 0; j < nr; j++)
900 { 900 {
901 double atmp = 0.; 901 double atmp = 0.;
902 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 902 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
903 atmp += fabs(data(i)); 903 atmp += fabs (data (i));
904 if (atmp > anorm) 904 if (atmp > anorm)
905 anorm = atmp; 905 anorm = atmp;
906 } 906 }
907 } 907 }
908 908
923 { 923 {
924 nz2 *= 2; 924 nz2 *= 2;
925 retval.change_capacity (nz2); 925 retval.change_capacity (nz2);
926 } 926 }
927 927
928 retval.xcidx(i) = cx; 928 retval.xcidx (i) = cx;
929 retval.xridx(cx) = i; 929 retval.xridx (cx) = i;
930 retval.xdata(cx) = 1.0; 930 retval.xdata (cx) = 1.0;
931 cx++; 931 cx++;
932 932
933 // iterate accross columns of input matrix 933 // iterate accross columns of input matrix
934 for (octave_idx_type j = i+1; j < nr; j++) 934 for (octave_idx_type j = i+1; j < nr; j++)
935 { 935 {
936 double v = 0.; 936 double v = 0.;
937 // iterate to calculate sum 937 // iterate to calculate sum
938 octave_idx_type colXp = retval.xcidx(i); 938 octave_idx_type colXp = retval.xcidx (i);
939 octave_idx_type colUp = cidx(j); 939 octave_idx_type colUp = cidx (j);
940 octave_idx_type rpX, rpU; 940 octave_idx_type rpX, rpU;
941 941
942 if (cidx(j) == cidx(j+1)) 942 if (cidx (j) == cidx (j+1))
943 { 943 {
944 (*current_liboctave_error_handler) 944 (*current_liboctave_error_handler)
945 ("division by zero"); 945 ("division by zero");
946 goto inverse_singular; 946 goto inverse_singular;
947 } 947 }
948 948
949 do 949 do
950 { 950 {
951 octave_quit (); 951 octave_quit ();
952 rpX = retval.xridx(colXp); 952 rpX = retval.xridx (colXp);
953 rpU = ridx(colUp); 953 rpU = ridx (colUp);
954 954
955 if (rpX < rpU) 955 if (rpX < rpU)
956 colXp++; 956 colXp++;
957 else if (rpX > rpU) 957 else if (rpX > rpU)
958 colUp++; 958 colUp++;
959 else 959 else
960 { 960 {
961 v -= retval.xdata(colXp) * data(colUp); 961 v -= retval.xdata (colXp) * data (colUp);
962 colXp++; 962 colXp++;
963 colUp++; 963 colUp++;
964 } 964 }
965 } while ((rpX<j) && (rpU<j) && 965 } while ((rpX<j) && (rpU<j) &&
966 (colXp<cx) && (colUp<nz)); 966 (colXp<cx) && (colUp<nz));
967 967
968 // get A(m,m) 968 // get A(m,m)
969 if (typ == MatrixType::Upper) 969 if (typ == MatrixType::Upper)
970 colUp = cidx(j+1) - 1; 970 colUp = cidx (j+1) - 1;
971 else 971 else
972 colUp = cidx(j); 972 colUp = cidx (j);
973 double pivot = data(colUp); 973 double pivot = data (colUp);
974 if (pivot == 0. || ridx(colUp) != j) 974 if (pivot == 0. || ridx (colUp) != j)
975 { 975 {
976 (*current_liboctave_error_handler) 976 (*current_liboctave_error_handler)
977 ("division by zero"); 977 ("division by zero");
978 goto inverse_singular; 978 goto inverse_singular;
979 } 979 }
984 { 984 {
985 nz2 *= 2; 985 nz2 *= 2;
986 retval.change_capacity (nz2); 986 retval.change_capacity (nz2);
987 } 987 }
988 988
989 retval.xridx(cx) = j; 989 retval.xridx (cx) = j;
990 retval.xdata(cx) = v / pivot; 990 retval.xdata (cx) = v / pivot;
991 cx++; 991 cx++;
992 } 992 }
993 } 993 }
994 994
995 // get A(m,m) 995 // get A(m,m)
996 octave_idx_type colUp; 996 octave_idx_type colUp;
997 if (typ == MatrixType::Upper) 997 if (typ == MatrixType::Upper)
998 colUp = cidx(i+1) - 1; 998 colUp = cidx (i+1) - 1;
999 else 999 else
1000 colUp = cidx(i); 1000 colUp = cidx (i);
1001 double pivot = data(colUp); 1001 double pivot = data (colUp);
1002 if (pivot == 0. || ridx(colUp) != i) 1002 if (pivot == 0. || ridx (colUp) != i)
1003 { 1003 {
1004 (*current_liboctave_error_handler) ("division by zero"); 1004 (*current_liboctave_error_handler) ("division by zero");
1005 goto inverse_singular; 1005 goto inverse_singular;
1006 } 1006 }
1007 1007
1008 if (pivot != 1.0) 1008 if (pivot != 1.0)
1009 for (octave_idx_type j = cx_colstart; j < cx; j++) 1009 for (octave_idx_type j = cx_colstart; j < cx; j++)
1010 retval.xdata(j) /= pivot; 1010 retval.xdata (j) /= pivot;
1011 } 1011 }
1012 retval.xcidx(nr) = cx; 1012 retval.xcidx (nr) = cx;
1013 retval.maybe_compress (); 1013 retval.maybe_compress ();
1014 } 1014 }
1015 else 1015 else
1016 { 1016 {
1017 octave_idx_type nz = nnz (); 1017 octave_idx_type nz = nnz ();
1051 for (octave_idx_type j = iidx+1; j < nr; j++) 1051 for (octave_idx_type j = iidx+1; j < nr; j++)
1052 { 1052 {
1053 double v = 0.; 1053 double v = 0.;
1054 octave_idx_type jidx = perm[j]; 1054 octave_idx_type jidx = perm[j];
1055 // iterate to calculate sum 1055 // iterate to calculate sum
1056 for (octave_idx_type k = cidx(jidx); 1056 for (octave_idx_type k = cidx (jidx);
1057 k < cidx(jidx+1); k++) 1057 k < cidx (jidx+1); k++)
1058 { 1058 {
1059 octave_quit (); 1059 octave_quit ();
1060 v -= work[ridx(k)] * data(k); 1060 v -= work[ridx (k)] * data (k);
1061 } 1061 }
1062 1062
1063 // get A(m,m) 1063 // get A(m,m)
1064 double pivot; 1064 double pivot;
1065 if (typ == MatrixType::Permuted_Upper) 1065 if (typ == MatrixType::Permuted_Upper)
1066 pivot = data(cidx(jidx+1) - 1); 1066 pivot = data (cidx (jidx+1) - 1);
1067 else 1067 else
1068 pivot = data(cidx(jidx)); 1068 pivot = data (cidx (jidx));
1069 if (pivot == 0.) 1069 if (pivot == 0.)
1070 { 1070 {
1071 (*current_liboctave_error_handler) 1071 (*current_liboctave_error_handler)
1072 ("division by zero"); 1072 ("division by zero");
1073 goto inverse_singular; 1073 goto inverse_singular;
1077 } 1077 }
1078 1078
1079 // get A(m,m) 1079 // get A(m,m)
1080 octave_idx_type colUp; 1080 octave_idx_type colUp;
1081 if (typ == MatrixType::Permuted_Upper) 1081 if (typ == MatrixType::Permuted_Upper)
1082 colUp = cidx(perm[iidx]+1) - 1; 1082 colUp = cidx (perm[iidx]+1) - 1;
1083 else 1083 else
1084 colUp = cidx(perm[iidx]); 1084 colUp = cidx (perm[iidx]);
1085 1085
1086 double pivot = data(colUp); 1086 double pivot = data (colUp);
1087 if (pivot == 0.) 1087 if (pivot == 0.)
1088 { 1088 {
1089 (*current_liboctave_error_handler) 1089 (*current_liboctave_error_handler)
1090 ("division by zero"); 1090 ("division by zero");
1091 goto inverse_singular; 1091 goto inverse_singular;
1104 { 1104 {
1105 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2); 1105 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1106 retval.change_capacity (nz2); 1106 retval.change_capacity (nz2);
1107 } 1107 }
1108 1108
1109 retval.xcidx(i) = cx; 1109 retval.xcidx (i) = cx;
1110 for (octave_idx_type j = iidx; j < nr; j++) 1110 for (octave_idx_type j = iidx; j < nr; j++)
1111 if (work[j] != 0.) 1111 if (work[j] != 0.)
1112 { 1112 {
1113 retval.xridx(cx) = j; 1113 retval.xridx (cx) = j;
1114 retval.xdata(cx++) = work[j]; 1114 retval.xdata (cx++) = work[j];
1115 } 1115 }
1116 } 1116 }
1117 1117
1118 retval.xcidx(nr) = cx; 1118 retval.xcidx (nr) = cx;
1119 retval.maybe_compress (); 1119 retval.maybe_compress ();
1120 } 1120 }
1121 1121
1122 if (calccond) 1122 if (calccond)
1123 { 1123 {
1124 // Calculate the 1-norm of inverse matrix for rcond calculation 1124 // Calculate the 1-norm of inverse matrix for rcond calculation
1125 for (octave_idx_type j = 0; j < nr; j++) 1125 for (octave_idx_type j = 0; j < nr; j++)
1126 { 1126 {
1127 double atmp = 0.; 1127 double atmp = 0.;
1128 for (octave_idx_type i = retval.cidx(j); 1128 for (octave_idx_type i = retval.cidx (j);
1129 i < retval.cidx(j+1); i++) 1129 i < retval.cidx (j+1); i++)
1130 atmp += fabs(retval.data(i)); 1130 atmp += fabs (retval.data (i));
1131 if (atmp > ainvnorm) 1131 if (atmp > ainvnorm)
1132 ainvnorm = atmp; 1132 ainvnorm = atmp;
1133 } 1133 }
1134 1134
1135 rcond = 1. / ainvnorm / anorm; 1135 rcond = 1. / ainvnorm / anorm;
1173 rcond = fact.rcond (); 1173 rcond = fact.rcond ();
1174 if (info == 0) 1174 if (info == 0)
1175 { 1175 {
1176 double rcond2; 1176 double rcond2;
1177 SparseMatrix Q = fact.Q (); 1177 SparseMatrix Q = fact.Q ();
1178 SparseMatrix InvL = fact.L ().transpose ().tinverse(tmp_typ, 1178 SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
1179 info, rcond2, true, false); 1179 info, rcond2, true, false);
1180 ret = Q * InvL.transpose () * InvL * Q.transpose (); 1180 ret = Q * InvL.transpose () * InvL * Q.transpose ();
1181 } 1181 }
1182 else 1182 else
1183 { 1183 {
1196 1196
1197 MatrixType tmp_typ (MatrixType::Upper); 1197 MatrixType tmp_typ (MatrixType::Upper);
1198 SparseLU fact (*this, Qinit, Matrix (), false, false); 1198 SparseLU fact (*this, Qinit, Matrix (), false, false);
1199 rcond = fact.rcond (); 1199 rcond = fact.rcond ();
1200 double rcond2; 1200 double rcond2;
1201 SparseMatrix InvL = fact.L ().transpose ().tinverse(tmp_typ, 1201 SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
1202 info, rcond2, true, false); 1202 info, rcond2, true, false);
1203 SparseMatrix InvU = fact.U ().tinverse(tmp_typ, info, rcond2, 1203 SparseMatrix InvU = fact.U ().tinverse (tmp_typ, info, rcond2,
1204 true, false).transpose (); 1204 true, false).transpose ();
1205 ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr (); 1205 ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
1206 } 1206 }
1207 } 1207 }
1208 1208
1372 for (octave_idx_type i = 0; i < nm; i++) 1372 for (octave_idx_type i = 0; i < nm; i++)
1373 retval(i,j) = b(i,j) / data (i); 1373 retval(i,j) = b(i,j) / data (i);
1374 else 1374 else
1375 for (octave_idx_type j = 0; j < b.cols (); j++) 1375 for (octave_idx_type j = 0; j < b.cols (); j++)
1376 for (octave_idx_type k = 0; k < nc; k++) 1376 for (octave_idx_type k = 0; k < nc; k++)
1377 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 1377 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1378 retval(k,j) = b(ridx(i),j) / data (i); 1378 retval(k,j) = b(ridx (i),j) / data (i);
1379 1379
1380 if (calc_cond) 1380 if (calc_cond)
1381 { 1381 {
1382 double dmax = 0., dmin = octave_Inf; 1382 double dmax = 0., dmin = octave_Inf;
1383 for (octave_idx_type i = 0; i < nm; i++) 1383 for (octave_idx_type i = 0; i < nm; i++)
1384 { 1384 {
1385 double tmp = fabs(data(i)); 1385 double tmp = fabs (data (i));
1386 if (tmp > dmax) 1386 if (tmp > dmax)
1387 dmax = tmp; 1387 dmax = tmp;
1388 if (tmp < dmin) 1388 if (tmp < dmin)
1389 dmin = tmp; 1389 dmin = tmp;
1390 } 1390 }
1428 { 1428 {
1429 octave_idx_type b_nc = b.cols (); 1429 octave_idx_type b_nc = b.cols ();
1430 octave_idx_type b_nz = b.nnz (); 1430 octave_idx_type b_nz = b.nnz ();
1431 retval = SparseMatrix (nc, b_nc, b_nz); 1431 retval = SparseMatrix (nc, b_nc, b_nz);
1432 1432
1433 retval.xcidx(0) = 0; 1433 retval.xcidx (0) = 0;
1434 octave_idx_type ii = 0; 1434 octave_idx_type ii = 0;
1435 if (typ == MatrixType::Diagonal) 1435 if (typ == MatrixType::Diagonal)
1436 for (octave_idx_type j = 0; j < b_nc; j++) 1436 for (octave_idx_type j = 0; j < b_nc; j++)
1437 { 1437 {
1438 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 1438 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1439 { 1439 {
1440 if (b.ridx(i) >= nm) 1440 if (b.ridx (i) >= nm)
1441 break; 1441 break;
1442 retval.xridx (ii) = b.ridx(i); 1442 retval.xridx (ii) = b.ridx (i);
1443 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); 1443 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1444 } 1444 }
1445 retval.xcidx(j+1) = ii; 1445 retval.xcidx (j+1) = ii;
1446 } 1446 }
1447 else 1447 else
1448 for (octave_idx_type j = 0; j < b_nc; j++) 1448 for (octave_idx_type j = 0; j < b_nc; j++)
1449 { 1449 {
1450 for (octave_idx_type l = 0; l < nc; l++) 1450 for (octave_idx_type l = 0; l < nc; l++)
1451 for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) 1451 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1452 { 1452 {
1453 bool found = false; 1453 bool found = false;
1454 octave_idx_type k; 1454 octave_idx_type k;
1455 for (k = b.cidx(j); k < b.cidx(j+1); k++) 1455 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1456 if (ridx(i) == b.ridx(k)) 1456 if (ridx (i) == b.ridx (k))
1457 { 1457 {
1458 found = true; 1458 found = true;
1459 break; 1459 break;
1460 } 1460 }
1461 if (found) 1461 if (found)
1462 { 1462 {
1463 retval.xridx (ii) = l; 1463 retval.xridx (ii) = l;
1464 retval.xdata (ii++) = b.data(k) / data (i); 1464 retval.xdata (ii++) = b.data (k) / data (i);
1465 } 1465 }
1466 } 1466 }
1467 retval.xcidx(j+1) = ii; 1467 retval.xcidx (j+1) = ii;
1468 } 1468 }
1469 1469
1470 if (calc_cond) 1470 if (calc_cond)
1471 { 1471 {
1472 double dmax = 0., dmin = octave_Inf; 1472 double dmax = 0., dmin = octave_Inf;
1473 for (octave_idx_type i = 0; i < nm; i++) 1473 for (octave_idx_type i = 0; i < nm; i++)
1474 { 1474 {
1475 double tmp = fabs(data(i)); 1475 double tmp = fabs (data (i));
1476 if (tmp > dmax) 1476 if (tmp > dmax)
1477 dmax = tmp; 1477 dmax = tmp;
1478 if (tmp < dmin) 1478 if (tmp < dmin)
1479 dmin = tmp; 1479 dmin = tmp;
1480 } 1480 }
1522 for (octave_idx_type i = 0; i < nm; i++) 1522 for (octave_idx_type i = 0; i < nm; i++)
1523 retval(i,j) = b(i,j) / data (i); 1523 retval(i,j) = b(i,j) / data (i);
1524 else 1524 else
1525 for (octave_idx_type j = 0; j < b.cols (); j++) 1525 for (octave_idx_type j = 0; j < b.cols (); j++)
1526 for (octave_idx_type k = 0; k < nc; k++) 1526 for (octave_idx_type k = 0; k < nc; k++)
1527 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 1527 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1528 retval(k,j) = b(ridx(i),j) / data (i); 1528 retval(k,j) = b(ridx (i),j) / data (i);
1529 1529
1530 if (calc_cond) 1530 if (calc_cond)
1531 { 1531 {
1532 double dmax = 0., dmin = octave_Inf; 1532 double dmax = 0., dmin = octave_Inf;
1533 for (octave_idx_type i = 0; i < nm; i++) 1533 for (octave_idx_type i = 0; i < nm; i++)
1534 { 1534 {
1535 double tmp = fabs(data(i)); 1535 double tmp = fabs (data (i));
1536 if (tmp > dmax) 1536 if (tmp > dmax)
1537 dmax = tmp; 1537 dmax = tmp;
1538 if (tmp < dmin) 1538 if (tmp < dmin)
1539 dmin = tmp; 1539 dmin = tmp;
1540 } 1540 }
1578 { 1578 {
1579 octave_idx_type b_nc = b.cols (); 1579 octave_idx_type b_nc = b.cols ();
1580 octave_idx_type b_nz = b.nnz (); 1580 octave_idx_type b_nz = b.nnz ();
1581 retval = SparseComplexMatrix (nc, b_nc, b_nz); 1581 retval = SparseComplexMatrix (nc, b_nc, b_nz);
1582 1582
1583 retval.xcidx(0) = 0; 1583 retval.xcidx (0) = 0;
1584 octave_idx_type ii = 0; 1584 octave_idx_type ii = 0;
1585 if (typ == MatrixType::Diagonal) 1585 if (typ == MatrixType::Diagonal)
1586 for (octave_idx_type j = 0; j < b.cols (); j++) 1586 for (octave_idx_type j = 0; j < b.cols (); j++)
1587 { 1587 {
1588 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 1588 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1589 { 1589 {
1590 if (b.ridx(i) >= nm) 1590 if (b.ridx (i) >= nm)
1591 break; 1591 break;
1592 retval.xridx (ii) = b.ridx(i); 1592 retval.xridx (ii) = b.ridx (i);
1593 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); 1593 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1594 } 1594 }
1595 retval.xcidx(j+1) = ii; 1595 retval.xcidx (j+1) = ii;
1596 } 1596 }
1597 else 1597 else
1598 for (octave_idx_type j = 0; j < b.cols (); j++) 1598 for (octave_idx_type j = 0; j < b.cols (); j++)
1599 { 1599 {
1600 for (octave_idx_type l = 0; l < nc; l++) 1600 for (octave_idx_type l = 0; l < nc; l++)
1601 for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) 1601 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1602 { 1602 {
1603 bool found = false; 1603 bool found = false;
1604 octave_idx_type k; 1604 octave_idx_type k;
1605 for (k = b.cidx(j); k < b.cidx(j+1); k++) 1605 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1606 if (ridx(i) == b.ridx(k)) 1606 if (ridx (i) == b.ridx (k))
1607 { 1607 {
1608 found = true; 1608 found = true;
1609 break; 1609 break;
1610 } 1610 }
1611 if (found) 1611 if (found)
1612 { 1612 {
1613 retval.xridx (ii) = l; 1613 retval.xridx (ii) = l;
1614 retval.xdata (ii++) = b.data(k) / data (i); 1614 retval.xdata (ii++) = b.data (k) / data (i);
1615 } 1615 }
1616 } 1616 }
1617 retval.xcidx(j+1) = ii; 1617 retval.xcidx (j+1) = ii;
1618 } 1618 }
1619 1619
1620 if (calc_cond) 1620 if (calc_cond)
1621 { 1621 {
1622 double dmax = 0., dmin = octave_Inf; 1622 double dmax = 0., dmin = octave_Inf;
1623 for (octave_idx_type i = 0; i < nm; i++) 1623 for (octave_idx_type i = 0; i < nm; i++)
1624 { 1624 {
1625 double tmp = fabs(data(i)); 1625 double tmp = fabs (data (i));
1626 if (tmp > dmax) 1626 if (tmp > dmax)
1627 dmax = tmp; 1627 dmax = tmp;
1628 if (tmp < dmin) 1628 if (tmp < dmin)
1629 dmin = tmp; 1629 dmin = tmp;
1630 } 1630 }
1676 { 1676 {
1677 // Calculate the 1-norm of matrix for rcond calculation 1677 // Calculate the 1-norm of matrix for rcond calculation
1678 for (octave_idx_type j = 0; j < nc; j++) 1678 for (octave_idx_type j = 0; j < nc; j++)
1679 { 1679 {
1680 double atmp = 0.; 1680 double atmp = 0.;
1681 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 1681 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1682 atmp += fabs(data(i)); 1682 atmp += fabs (data (i));
1683 if (atmp > anorm) 1683 if (atmp > anorm)
1684 anorm = atmp; 1684 anorm = atmp;
1685 } 1685 }
1686 } 1686 }
1687 1687
1702 { 1702 {
1703 octave_idx_type kidx = perm[k]; 1703 octave_idx_type kidx = perm[k];
1704 1704
1705 if (work[k] != 0.) 1705 if (work[k] != 0.)
1706 { 1706 {
1707 if (ridx(cidx(kidx+1)-1) != k || 1707 if (ridx (cidx (kidx+1)-1) != k ||
1708 data(cidx(kidx+1)-1) == 0.) 1708 data (cidx (kidx+1)-1) == 0.)
1709 { 1709 {
1710 err = -2; 1710 err = -2;
1711 goto triangular_error; 1711 goto triangular_error;
1712 } 1712 }
1713 1713
1714 double tmp = work[k] / data(cidx(kidx+1)-1); 1714 double tmp = work[k] / data (cidx (kidx+1)-1);
1715 work[k] = tmp; 1715 work[k] = tmp;
1716 for (octave_idx_type i = cidx(kidx); 1716 for (octave_idx_type i = cidx (kidx);
1717 i < cidx(kidx+1)-1; i++) 1717 i < cidx (kidx+1)-1; i++)
1718 { 1718 {
1719 octave_idx_type iidx = ridx(i); 1719 octave_idx_type iidx = ridx (i);
1720 work[iidx] = work[iidx] - tmp * data(i); 1720 work[iidx] = work[iidx] - tmp * data (i);
1721 } 1721 }
1722 } 1722 }
1723 } 1723 }
1724 1724
1725 for (octave_idx_type i = 0; i < nc; i++) 1725 for (octave_idx_type i = 0; i < nc; i++)
1740 { 1740 {
1741 octave_idx_type iidx = perm[k]; 1741 octave_idx_type iidx = perm[k];
1742 1742
1743 if (work[k] != 0.) 1743 if (work[k] != 0.)
1744 { 1744 {
1745 double tmp = work[k] / data(cidx(iidx+1)-1); 1745 double tmp = work[k] / data (cidx (iidx+1)-1);
1746 work[k] = tmp; 1746 work[k] = tmp;
1747 for (octave_idx_type i = cidx(iidx); 1747 for (octave_idx_type i = cidx (iidx);
1748 i < cidx(iidx+1)-1; i++) 1748 i < cidx (iidx+1)-1; i++)
1749 { 1749 {
1750 octave_idx_type idx2 = ridx(i); 1750 octave_idx_type idx2 = ridx (i);
1751 work[idx2] = work[idx2] - tmp * data(i); 1751 work[idx2] = work[idx2] - tmp * data (i);
1752 } 1752 }
1753 } 1753 }
1754 } 1754 }
1755 double atmp = 0; 1755 double atmp = 0;
1756 for (octave_idx_type i = 0; i < j+1; i++) 1756 for (octave_idx_type i = 0; i < j+1; i++)
1757 { 1757 {
1758 atmp += fabs(work[i]); 1758 atmp += fabs (work[i]);
1759 work[i] = 0.; 1759 work[i] = 0.;
1760 } 1760 }
1761 if (atmp > ainvnorm) 1761 if (atmp > ainvnorm)
1762 ainvnorm = atmp; 1762 ainvnorm = atmp;
1763 } 1763 }
1778 1778
1779 for (octave_idx_type k = nc-1; k >= 0; k--) 1779 for (octave_idx_type k = nc-1; k >= 0; k--)
1780 { 1780 {
1781 if (work[k] != 0.) 1781 if (work[k] != 0.)
1782 { 1782 {
1783 if (ridx(cidx(k+1)-1) != k || 1783 if (ridx (cidx (k+1)-1) != k ||
1784 data(cidx(k+1)-1) == 0.) 1784 data (cidx (k+1)-1) == 0.)
1785 { 1785 {
1786 err = -2; 1786 err = -2;
1787 goto triangular_error; 1787 goto triangular_error;
1788 } 1788 }
1789 1789
1790 double tmp = work[k] / data(cidx(k+1)-1); 1790 double tmp = work[k] / data (cidx (k+1)-1);
1791 work[k] = tmp; 1791 work[k] = tmp;
1792 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) 1792 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1793 { 1793 {
1794 octave_idx_type iidx = ridx(i); 1794 octave_idx_type iidx = ridx (i);
1795 work[iidx] = work[iidx] - tmp * data(i); 1795 work[iidx] = work[iidx] - tmp * data (i);
1796 } 1796 }
1797 } 1797 }
1798 } 1798 }
1799 1799
1800 for (octave_idx_type i = 0; i < nc; i++) 1800 for (octave_idx_type i = 0; i < nc; i++)
1813 1813
1814 for (octave_idx_type k = j; k >= 0; k--) 1814 for (octave_idx_type k = j; k >= 0; k--)
1815 { 1815 {
1816 if (work[k] != 0.) 1816 if (work[k] != 0.)
1817 { 1817 {
1818 double tmp = work[k] / data(cidx(k+1)-1); 1818 double tmp = work[k] / data (cidx (k+1)-1);
1819 work[k] = tmp; 1819 work[k] = tmp;
1820 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) 1820 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1821 { 1821 {
1822 octave_idx_type iidx = ridx(i); 1822 octave_idx_type iidx = ridx (i);
1823 work[iidx] = work[iidx] - tmp * data(i); 1823 work[iidx] = work[iidx] - tmp * data (i);
1824 } 1824 }
1825 } 1825 }
1826 } 1826 }
1827 double atmp = 0; 1827 double atmp = 0;
1828 for (octave_idx_type i = 0; i < j+1; i++) 1828 for (octave_idx_type i = 0; i < j+1; i++)
1829 { 1829 {
1830 atmp += fabs(work[i]); 1830 atmp += fabs (work[i]);
1831 work[i] = 0.; 1831 work[i] = 0.;
1832 } 1832 }
1833 if (atmp > ainvnorm) 1833 if (atmp > ainvnorm)
1834 ainvnorm = atmp; 1834 ainvnorm = atmp;
1835 } 1835 }
1910 { 1910 {
1911 // Calculate the 1-norm of matrix for rcond calculation 1911 // Calculate the 1-norm of matrix for rcond calculation
1912 for (octave_idx_type j = 0; j < nc; j++) 1912 for (octave_idx_type j = 0; j < nc; j++)
1913 { 1913 {
1914 double atmp = 0.; 1914 double atmp = 0.;
1915 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 1915 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1916 atmp += fabs(data(i)); 1916 atmp += fabs (data (i));
1917 if (atmp > anorm) 1917 if (atmp > anorm)
1918 anorm = atmp; 1918 anorm = atmp;
1919 } 1919 }
1920 } 1920 }
1921 1921
1922 octave_idx_type b_nc = b.cols (); 1922 octave_idx_type b_nc = b.cols ();
1923 octave_idx_type b_nz = b.nnz (); 1923 octave_idx_type b_nz = b.nnz ();
1924 retval = SparseMatrix (nc, b_nc, b_nz); 1924 retval = SparseMatrix (nc, b_nc, b_nz);
1925 retval.xcidx(0) = 0; 1925 retval.xcidx (0) = 0;
1926 octave_idx_type ii = 0; 1926 octave_idx_type ii = 0;
1927 octave_idx_type x_nz = b_nz; 1927 octave_idx_type x_nz = b_nz;
1928 1928
1929 if (typ == MatrixType::Permuted_Upper) 1929 if (typ == MatrixType::Permuted_Upper)
1930 { 1930 {
1937 1937
1938 for (octave_idx_type j = 0; j < b_nc; j++) 1938 for (octave_idx_type j = 0; j < b_nc; j++)
1939 { 1939 {
1940 for (octave_idx_type i = 0; i < nm; i++) 1940 for (octave_idx_type i = 0; i < nm; i++)
1941 work[i] = 0.; 1941 work[i] = 0.;
1942 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 1942 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1943 work[b.ridx(i)] = b.data(i); 1943 work[b.ridx (i)] = b.data (i);
1944 1944
1945 for (octave_idx_type k = nc-1; k >= 0; k--) 1945 for (octave_idx_type k = nc-1; k >= 0; k--)
1946 { 1946 {
1947 octave_idx_type kidx = perm[k]; 1947 octave_idx_type kidx = perm[k];
1948 1948
1949 if (work[k] != 0.) 1949 if (work[k] != 0.)
1950 { 1950 {
1951 if (ridx(cidx(kidx+1)-1) != k || 1951 if (ridx (cidx (kidx+1)-1) != k ||
1952 data(cidx(kidx+1)-1) == 0.) 1952 data (cidx (kidx+1)-1) == 0.)
1953 { 1953 {
1954 err = -2; 1954 err = -2;
1955 goto triangular_error; 1955 goto triangular_error;
1956 } 1956 }
1957 1957
1958 double tmp = work[k] / data(cidx(kidx+1)-1); 1958 double tmp = work[k] / data (cidx (kidx+1)-1);
1959 work[k] = tmp; 1959 work[k] = tmp;
1960 for (octave_idx_type i = cidx(kidx); 1960 for (octave_idx_type i = cidx (kidx);
1961 i < cidx(kidx+1)-1; i++) 1961 i < cidx (kidx+1)-1; i++)
1962 { 1962 {
1963 octave_idx_type iidx = ridx(i); 1963 octave_idx_type iidx = ridx (i);
1964 work[iidx] = work[iidx] - tmp * data(i); 1964 work[iidx] = work[iidx] - tmp * data (i);
1965 } 1965 }
1966 } 1966 }
1967 } 1967 }
1968 1968
1969 // Count non-zeros in work vector and adjust space in 1969 // Count non-zeros in work vector and adjust space in
1982 } 1982 }
1983 1983
1984 for (octave_idx_type i = 0; i < nc; i++) 1984 for (octave_idx_type i = 0; i < nc; i++)
1985 if (work[rperm[i]] != 0.) 1985 if (work[rperm[i]] != 0.)
1986 { 1986 {
1987 retval.xridx(ii) = i; 1987 retval.xridx (ii) = i;
1988 retval.xdata(ii++) = work[rperm[i]]; 1988 retval.xdata (ii++) = work[rperm[i]];
1989 } 1989 }
1990 retval.xcidx(j+1) = ii; 1990 retval.xcidx (j+1) = ii;
1991 } 1991 }
1992 1992
1993 retval.maybe_compress (); 1993 retval.maybe_compress ();
1994 1994
1995 if (calc_cond) 1995 if (calc_cond)
2006 { 2006 {
2007 octave_idx_type iidx = perm[k]; 2007 octave_idx_type iidx = perm[k];
2008 2008
2009 if (work[k] != 0.) 2009 if (work[k] != 0.)
2010 { 2010 {
2011 double tmp = work[k] / data(cidx(iidx+1)-1); 2011 double tmp = work[k] / data (cidx (iidx+1)-1);
2012 work[k] = tmp; 2012 work[k] = tmp;
2013 for (octave_idx_type i = cidx(iidx); 2013 for (octave_idx_type i = cidx (iidx);
2014 i < cidx(iidx+1)-1; i++) 2014 i < cidx (iidx+1)-1; i++)
2015 { 2015 {
2016 octave_idx_type idx2 = ridx(i); 2016 octave_idx_type idx2 = ridx (i);
2017 work[idx2] = work[idx2] - tmp * data(i); 2017 work[idx2] = work[idx2] - tmp * data (i);
2018 } 2018 }
2019 } 2019 }
2020 } 2020 }
2021 double atmp = 0; 2021 double atmp = 0;
2022 for (octave_idx_type i = 0; i < j+1; i++) 2022 for (octave_idx_type i = 0; i < j+1; i++)
2023 { 2023 {
2024 atmp += fabs(work[i]); 2024 atmp += fabs (work[i]);
2025 work[i] = 0.; 2025 work[i] = 0.;
2026 } 2026 }
2027 if (atmp > ainvnorm) 2027 if (atmp > ainvnorm)
2028 ainvnorm = atmp; 2028 ainvnorm = atmp;
2029 } 2029 }
2036 2036
2037 for (octave_idx_type j = 0; j < b_nc; j++) 2037 for (octave_idx_type j = 0; j < b_nc; j++)
2038 { 2038 {
2039 for (octave_idx_type i = 0; i < nm; i++) 2039 for (octave_idx_type i = 0; i < nm; i++)
2040 work[i] = 0.; 2040 work[i] = 0.;
2041 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 2041 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2042 work[b.ridx(i)] = b.data(i); 2042 work[b.ridx (i)] = b.data (i);
2043 2043
2044 for (octave_idx_type k = nc-1; k >= 0; k--) 2044 for (octave_idx_type k = nc-1; k >= 0; k--)
2045 { 2045 {
2046 if (work[k] != 0.) 2046 if (work[k] != 0.)
2047 { 2047 {
2048 if (ridx(cidx(k+1)-1) != k || 2048 if (ridx (cidx (k+1)-1) != k ||
2049 data(cidx(k+1)-1) == 0.) 2049 data (cidx (k+1)-1) == 0.)
2050 { 2050 {
2051 err = -2; 2051 err = -2;
2052 goto triangular_error; 2052 goto triangular_error;
2053 } 2053 }
2054 2054
2055 double tmp = work[k] / data(cidx(k+1)-1); 2055 double tmp = work[k] / data (cidx (k+1)-1);
2056 work[k] = tmp; 2056 work[k] = tmp;
2057 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) 2057 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2058 { 2058 {
2059 octave_idx_type iidx = ridx(i); 2059 octave_idx_type iidx = ridx (i);
2060 work[iidx] = work[iidx] - tmp * data(i); 2060 work[iidx] = work[iidx] - tmp * data (i);
2061 } 2061 }
2062 } 2062 }
2063 } 2063 }
2064 2064
2065 // Count non-zeros in work vector and adjust space in 2065 // Count non-zeros in work vector and adjust space in
2078 } 2078 }
2079 2079
2080 for (octave_idx_type i = 0; i < nc; i++) 2080 for (octave_idx_type i = 0; i < nc; i++)
2081 if (work[i] != 0.) 2081 if (work[i] != 0.)
2082 { 2082 {
2083 retval.xridx(ii) = i; 2083 retval.xridx (ii) = i;
2084 retval.xdata(ii++) = work[i]; 2084 retval.xdata (ii++) = work[i];
2085 } 2085 }
2086 retval.xcidx(j+1) = ii; 2086 retval.xcidx (j+1) = ii;
2087 } 2087 }
2088 2088
2089 retval.maybe_compress (); 2089 retval.maybe_compress ();
2090 2090
2091 if (calc_cond) 2091 if (calc_cond)
2100 2100
2101 for (octave_idx_type k = j; k >= 0; k--) 2101 for (octave_idx_type k = j; k >= 0; k--)
2102 { 2102 {
2103 if (work[k] != 0.) 2103 if (work[k] != 0.)
2104 { 2104 {
2105 double tmp = work[k] / data(cidx(k+1)-1); 2105 double tmp = work[k] / data (cidx (k+1)-1);
2106 work[k] = tmp; 2106 work[k] = tmp;
2107 for (octave_idx_type i = cidx(k); 2107 for (octave_idx_type i = cidx (k);
2108 i < cidx(k+1)-1; i++) 2108 i < cidx (k+1)-1; i++)
2109 { 2109 {
2110 octave_idx_type iidx = ridx(i); 2110 octave_idx_type iidx = ridx (i);
2111 work[iidx] = work[iidx] - tmp * data(i); 2111 work[iidx] = work[iidx] - tmp * data (i);
2112 } 2112 }
2113 } 2113 }
2114 } 2114 }
2115 double atmp = 0; 2115 double atmp = 0;
2116 for (octave_idx_type i = 0; i < j+1; i++) 2116 for (octave_idx_type i = 0; i < j+1; i++)
2117 { 2117 {
2118 atmp += fabs(work[i]); 2118 atmp += fabs (work[i]);
2119 work[i] = 0.; 2119 work[i] = 0.;
2120 } 2120 }
2121 if (atmp > ainvnorm) 2121 if (atmp > ainvnorm)
2122 ainvnorm = atmp; 2122 ainvnorm = atmp;
2123 } 2123 }
2198 { 2198 {
2199 // Calculate the 1-norm of matrix for rcond calculation 2199 // Calculate the 1-norm of matrix for rcond calculation
2200 for (octave_idx_type j = 0; j < nc; j++) 2200 for (octave_idx_type j = 0; j < nc; j++)
2201 { 2201 {
2202 double atmp = 0.; 2202 double atmp = 0.;
2203 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 2203 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2204 atmp += fabs(data(i)); 2204 atmp += fabs (data (i));
2205 if (atmp > anorm) 2205 if (atmp > anorm)
2206 anorm = atmp; 2206 anorm = atmp;
2207 } 2207 }
2208 } 2208 }
2209 2209
2224 { 2224 {
2225 octave_idx_type kidx = perm[k]; 2225 octave_idx_type kidx = perm[k];
2226 2226
2227 if (cwork[k] != 0.) 2227 if (cwork[k] != 0.)
2228 { 2228 {
2229 if (ridx(cidx(kidx+1)-1) != k || 2229 if (ridx (cidx (kidx+1)-1) != k ||
2230 data(cidx(kidx+1)-1) == 0.) 2230 data (cidx (kidx+1)-1) == 0.)
2231 { 2231 {
2232 err = -2; 2232 err = -2;
2233 goto triangular_error; 2233 goto triangular_error;
2234 } 2234 }
2235 2235
2236 Complex tmp = cwork[k] / data(cidx(kidx+1)-1); 2236 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2237 cwork[k] = tmp; 2237 cwork[k] = tmp;
2238 for (octave_idx_type i = cidx(kidx); 2238 for (octave_idx_type i = cidx (kidx);
2239 i < cidx(kidx+1)-1; i++) 2239 i < cidx (kidx+1)-1; i++)
2240 { 2240 {
2241 octave_idx_type iidx = ridx(i); 2241 octave_idx_type iidx = ridx (i);
2242 cwork[iidx] = cwork[iidx] - tmp * data(i); 2242 cwork[iidx] = cwork[iidx] - tmp * data (i);
2243 } 2243 }
2244 } 2244 }
2245 } 2245 }
2246 2246
2247 for (octave_idx_type i = 0; i < nc; i++) 2247 for (octave_idx_type i = 0; i < nc; i++)
2263 { 2263 {
2264 octave_idx_type iidx = perm[k]; 2264 octave_idx_type iidx = perm[k];
2265 2265
2266 if (work[k] != 0.) 2266 if (work[k] != 0.)
2267 { 2267 {
2268 double tmp = work[k] / data(cidx(iidx+1)-1); 2268 double tmp = work[k] / data (cidx (iidx+1)-1);
2269 work[k] = tmp; 2269 work[k] = tmp;
2270 for (octave_idx_type i = cidx(iidx); 2270 for (octave_idx_type i = cidx (iidx);
2271 i < cidx(iidx+1)-1; i++) 2271 i < cidx (iidx+1)-1; i++)
2272 { 2272 {
2273 octave_idx_type idx2 = ridx(i); 2273 octave_idx_type idx2 = ridx (i);
2274 work[idx2] = work[idx2] - tmp * data(i); 2274 work[idx2] = work[idx2] - tmp * data (i);
2275 } 2275 }
2276 } 2276 }
2277 } 2277 }
2278 double atmp = 0; 2278 double atmp = 0;
2279 for (octave_idx_type i = 0; i < j+1; i++) 2279 for (octave_idx_type i = 0; i < j+1; i++)
2280 { 2280 {
2281 atmp += fabs(work[i]); 2281 atmp += fabs (work[i]);
2282 work[i] = 0.; 2282 work[i] = 0.;
2283 } 2283 }
2284 if (atmp > ainvnorm) 2284 if (atmp > ainvnorm)
2285 ainvnorm = atmp; 2285 ainvnorm = atmp;
2286 } 2286 }
2301 2301
2302 for (octave_idx_type k = nc-1; k >= 0; k--) 2302 for (octave_idx_type k = nc-1; k >= 0; k--)
2303 { 2303 {
2304 if (cwork[k] != 0.) 2304 if (cwork[k] != 0.)
2305 { 2305 {
2306 if (ridx(cidx(k+1)-1) != k || 2306 if (ridx (cidx (k+1)-1) != k ||
2307 data(cidx(k+1)-1) == 0.) 2307 data (cidx (k+1)-1) == 0.)
2308 { 2308 {
2309 err = -2; 2309 err = -2;
2310 goto triangular_error; 2310 goto triangular_error;
2311 } 2311 }
2312 2312
2313 Complex tmp = cwork[k] / data(cidx(k+1)-1); 2313 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2314 cwork[k] = tmp; 2314 cwork[k] = tmp;
2315 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) 2315 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2316 { 2316 {
2317 octave_idx_type iidx = ridx(i); 2317 octave_idx_type iidx = ridx (i);
2318 cwork[iidx] = cwork[iidx] - tmp * data(i); 2318 cwork[iidx] = cwork[iidx] - tmp * data (i);
2319 } 2319 }
2320 } 2320 }
2321 } 2321 }
2322 2322
2323 for (octave_idx_type i = 0; i < nc; i++) 2323 for (octave_idx_type i = 0; i < nc; i++)
2337 2337
2338 for (octave_idx_type k = j; k >= 0; k--) 2338 for (octave_idx_type k = j; k >= 0; k--)
2339 { 2339 {
2340 if (work[k] != 0.) 2340 if (work[k] != 0.)
2341 { 2341 {
2342 double tmp = work[k] / data(cidx(k+1)-1); 2342 double tmp = work[k] / data (cidx (k+1)-1);
2343 work[k] = tmp; 2343 work[k] = tmp;
2344 for (octave_idx_type i = cidx(k); 2344 for (octave_idx_type i = cidx (k);
2345 i < cidx(k+1)-1; i++) 2345 i < cidx (k+1)-1; i++)
2346 { 2346 {
2347 octave_idx_type iidx = ridx(i); 2347 octave_idx_type iidx = ridx (i);
2348 work[iidx] = work[iidx] - tmp * data(i); 2348 work[iidx] = work[iidx] - tmp * data (i);
2349 } 2349 }
2350 } 2350 }
2351 } 2351 }
2352 double atmp = 0; 2352 double atmp = 0;
2353 for (octave_idx_type i = 0; i < j+1; i++) 2353 for (octave_idx_type i = 0; i < j+1; i++)
2354 { 2354 {
2355 atmp += fabs(work[i]); 2355 atmp += fabs (work[i]);
2356 work[i] = 0.; 2356 work[i] = 0.;
2357 } 2357 }
2358 if (atmp > ainvnorm) 2358 if (atmp > ainvnorm)
2359 ainvnorm = atmp; 2359 ainvnorm = atmp;
2360 } 2360 }
2435 { 2435 {
2436 // Calculate the 1-norm of matrix for rcond calculation 2436 // Calculate the 1-norm of matrix for rcond calculation
2437 for (octave_idx_type j = 0; j < nc; j++) 2437 for (octave_idx_type j = 0; j < nc; j++)
2438 { 2438 {
2439 double atmp = 0.; 2439 double atmp = 0.;
2440 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 2440 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2441 atmp += fabs(data(i)); 2441 atmp += fabs (data (i));
2442 if (atmp > anorm) 2442 if (atmp > anorm)
2443 anorm = atmp; 2443 anorm = atmp;
2444 } 2444 }
2445 } 2445 }
2446 2446
2447 octave_idx_type b_nc = b.cols (); 2447 octave_idx_type b_nc = b.cols ();
2448 octave_idx_type b_nz = b.nnz (); 2448 octave_idx_type b_nz = b.nnz ();
2449 retval = SparseComplexMatrix (nc, b_nc, b_nz); 2449 retval = SparseComplexMatrix (nc, b_nc, b_nz);
2450 retval.xcidx(0) = 0; 2450 retval.xcidx (0) = 0;
2451 octave_idx_type ii = 0; 2451 octave_idx_type ii = 0;
2452 octave_idx_type x_nz = b_nz; 2452 octave_idx_type x_nz = b_nz;
2453 2453
2454 if (typ == MatrixType::Permuted_Upper) 2454 if (typ == MatrixType::Permuted_Upper)
2455 { 2455 {
2462 2462
2463 for (octave_idx_type j = 0; j < b_nc; j++) 2463 for (octave_idx_type j = 0; j < b_nc; j++)
2464 { 2464 {
2465 for (octave_idx_type i = 0; i < nm; i++) 2465 for (octave_idx_type i = 0; i < nm; i++)
2466 cwork[i] = 0.; 2466 cwork[i] = 0.;
2467 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 2467 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2468 cwork[b.ridx(i)] = b.data(i); 2468 cwork[b.ridx (i)] = b.data (i);
2469 2469
2470 for (octave_idx_type k = nc-1; k >= 0; k--) 2470 for (octave_idx_type k = nc-1; k >= 0; k--)
2471 { 2471 {
2472 octave_idx_type kidx = perm[k]; 2472 octave_idx_type kidx = perm[k];
2473 2473
2474 if (cwork[k] != 0.) 2474 if (cwork[k] != 0.)
2475 { 2475 {
2476 if (ridx(cidx(kidx+1)-1) != k || 2476 if (ridx (cidx (kidx+1)-1) != k ||
2477 data(cidx(kidx+1)-1) == 0.) 2477 data (cidx (kidx+1)-1) == 0.)
2478 { 2478 {
2479 err = -2; 2479 err = -2;
2480 goto triangular_error; 2480 goto triangular_error;
2481 } 2481 }
2482 2482
2483 Complex tmp = cwork[k] / data(cidx(kidx+1)-1); 2483 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2484 cwork[k] = tmp; 2484 cwork[k] = tmp;
2485 for (octave_idx_type i = cidx(kidx); 2485 for (octave_idx_type i = cidx (kidx);
2486 i < cidx(kidx+1)-1; i++) 2486 i < cidx (kidx+1)-1; i++)
2487 { 2487 {
2488 octave_idx_type iidx = ridx(i); 2488 octave_idx_type iidx = ridx (i);
2489 cwork[iidx] = cwork[iidx] - tmp * data(i); 2489 cwork[iidx] = cwork[iidx] - tmp * data (i);
2490 } 2490 }
2491 } 2491 }
2492 } 2492 }
2493 2493
2494 // Count non-zeros in work vector and adjust space in 2494 // Count non-zeros in work vector and adjust space in
2507 } 2507 }
2508 2508
2509 for (octave_idx_type i = 0; i < nc; i++) 2509 for (octave_idx_type i = 0; i < nc; i++)
2510 if (cwork[rperm[i]] != 0.) 2510 if (cwork[rperm[i]] != 0.)
2511 { 2511 {
2512 retval.xridx(ii) = i; 2512 retval.xridx (ii) = i;
2513 retval.xdata(ii++) = cwork[rperm[i]]; 2513 retval.xdata (ii++) = cwork[rperm[i]];
2514 } 2514 }
2515 retval.xcidx(j+1) = ii; 2515 retval.xcidx (j+1) = ii;
2516 } 2516 }
2517 2517
2518 retval.maybe_compress (); 2518 retval.maybe_compress ();
2519 2519
2520 if (calc_cond) 2520 if (calc_cond)
2532 { 2532 {
2533 octave_idx_type iidx = perm[k]; 2533 octave_idx_type iidx = perm[k];
2534 2534
2535 if (work[k] != 0.) 2535 if (work[k] != 0.)
2536 { 2536 {
2537 double tmp = work[k] / data(cidx(iidx+1)-1); 2537 double tmp = work[k] / data (cidx (iidx+1)-1);
2538 work[k] = tmp; 2538 work[k] = tmp;
2539 for (octave_idx_type i = cidx(iidx); 2539 for (octave_idx_type i = cidx (iidx);
2540 i < cidx(iidx+1)-1; i++) 2540 i < cidx (iidx+1)-1; i++)
2541 { 2541 {
2542 octave_idx_type idx2 = ridx(i); 2542 octave_idx_type idx2 = ridx (i);
2543 work[idx2] = work[idx2] - tmp * data(i); 2543 work[idx2] = work[idx2] - tmp * data (i);
2544 } 2544 }
2545 } 2545 }
2546 } 2546 }
2547 double atmp = 0; 2547 double atmp = 0;
2548 for (octave_idx_type i = 0; i < j+1; i++) 2548 for (octave_idx_type i = 0; i < j+1; i++)
2549 { 2549 {
2550 atmp += fabs(work[i]); 2550 atmp += fabs (work[i]);
2551 work[i] = 0.; 2551 work[i] = 0.;
2552 } 2552 }
2553 if (atmp > ainvnorm) 2553 if (atmp > ainvnorm)
2554 ainvnorm = atmp; 2554 ainvnorm = atmp;
2555 } 2555 }
2562 2562
2563 for (octave_idx_type j = 0; j < b_nc; j++) 2563 for (octave_idx_type j = 0; j < b_nc; j++)
2564 { 2564 {
2565 for (octave_idx_type i = 0; i < nm; i++) 2565 for (octave_idx_type i = 0; i < nm; i++)
2566 cwork[i] = 0.; 2566 cwork[i] = 0.;
2567 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 2567 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2568 cwork[b.ridx(i)] = b.data(i); 2568 cwork[b.ridx (i)] = b.data (i);
2569 2569
2570 for (octave_idx_type k = nc-1; k >= 0; k--) 2570 for (octave_idx_type k = nc-1; k >= 0; k--)
2571 { 2571 {
2572 if (cwork[k] != 0.) 2572 if (cwork[k] != 0.)
2573 { 2573 {
2574 if (ridx(cidx(k+1)-1) != k || 2574 if (ridx (cidx (k+1)-1) != k ||
2575 data(cidx(k+1)-1) == 0.) 2575 data (cidx (k+1)-1) == 0.)
2576 { 2576 {
2577 err = -2; 2577 err = -2;
2578 goto triangular_error; 2578 goto triangular_error;
2579 } 2579 }
2580 2580
2581 Complex tmp = cwork[k] / data(cidx(k+1)-1); 2581 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2582 cwork[k] = tmp; 2582 cwork[k] = tmp;
2583 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) 2583 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2584 { 2584 {
2585 octave_idx_type iidx = ridx(i); 2585 octave_idx_type iidx = ridx (i);
2586 cwork[iidx] = cwork[iidx] - tmp * data(i); 2586 cwork[iidx] = cwork[iidx] - tmp * data (i);
2587 } 2587 }
2588 } 2588 }
2589 } 2589 }
2590 2590
2591 // Count non-zeros in work vector and adjust space in 2591 // Count non-zeros in work vector and adjust space in
2604 } 2604 }
2605 2605
2606 for (octave_idx_type i = 0; i < nc; i++) 2606 for (octave_idx_type i = 0; i < nc; i++)
2607 if (cwork[i] != 0.) 2607 if (cwork[i] != 0.)
2608 { 2608 {
2609 retval.xridx(ii) = i; 2609 retval.xridx (ii) = i;
2610 retval.xdata(ii++) = cwork[i]; 2610 retval.xdata (ii++) = cwork[i];
2611 } 2611 }
2612 retval.xcidx(j+1) = ii; 2612 retval.xcidx (j+1) = ii;
2613 } 2613 }
2614 2614
2615 retval.maybe_compress (); 2615 retval.maybe_compress ();
2616 2616
2617 if (calc_cond) 2617 if (calc_cond)
2627 2627
2628 for (octave_idx_type k = j; k >= 0; k--) 2628 for (octave_idx_type k = j; k >= 0; k--)
2629 { 2629 {
2630 if (work[k] != 0.) 2630 if (work[k] != 0.)
2631 { 2631 {
2632 double tmp = work[k] / data(cidx(k+1)-1); 2632 double tmp = work[k] / data (cidx (k+1)-1);
2633 work[k] = tmp; 2633 work[k] = tmp;
2634 for (octave_idx_type i = cidx(k); 2634 for (octave_idx_type i = cidx (k);
2635 i < cidx(k+1)-1; i++) 2635 i < cidx (k+1)-1; i++)
2636 { 2636 {
2637 octave_idx_type iidx = ridx(i); 2637 octave_idx_type iidx = ridx (i);
2638 work[iidx] = work[iidx] - tmp * data(i); 2638 work[iidx] = work[iidx] - tmp * data (i);
2639 } 2639 }
2640 } 2640 }
2641 } 2641 }
2642 double atmp = 0; 2642 double atmp = 0;
2643 for (octave_idx_type i = 0; i < j+1; i++) 2643 for (octave_idx_type i = 0; i < j+1; i++)
2644 { 2644 {
2645 atmp += fabs(work[i]); 2645 atmp += fabs (work[i]);
2646 work[i] = 0.; 2646 work[i] = 0.;
2647 } 2647 }
2648 if (atmp > ainvnorm) 2648 if (atmp > ainvnorm)
2649 ainvnorm = atmp; 2649 ainvnorm = atmp;
2650 } 2650 }
2726 { 2726 {
2727 // Calculate the 1-norm of matrix for rcond calculation 2727 // Calculate the 1-norm of matrix for rcond calculation
2728 for (octave_idx_type j = 0; j < nc; j++) 2728 for (octave_idx_type j = 0; j < nc; j++)
2729 { 2729 {
2730 double atmp = 0.; 2730 double atmp = 0.;
2731 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 2731 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2732 atmp += fabs(data(i)); 2732 atmp += fabs (data (i));
2733 if (atmp > anorm) 2733 if (atmp > anorm)
2734 anorm = atmp; 2734 anorm = atmp;
2735 } 2735 }
2736 } 2736 }
2737 2737
2754 if (work[k] != 0.) 2754 if (work[k] != 0.)
2755 { 2755 {
2756 octave_idx_type minr = nr; 2756 octave_idx_type minr = nr;
2757 octave_idx_type mini = 0; 2757 octave_idx_type mini = 0;
2758 2758
2759 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 2759 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2760 if (perm[ridx(i)] < minr) 2760 if (perm[ridx (i)] < minr)
2761 { 2761 {
2762 minr = perm[ridx(i)]; 2762 minr = perm[ridx (i)];
2763 mini = i; 2763 mini = i;
2764 } 2764 }
2765 2765
2766 if (minr != k || data(mini) == 0) 2766 if (minr != k || data (mini) == 0)
2767 { 2767 {
2768 err = -2; 2768 err = -2;
2769 goto triangular_error; 2769 goto triangular_error;
2770 } 2770 }
2771 2771
2772 double tmp = work[k] / data(mini); 2772 double tmp = work[k] / data (mini);
2773 work[k] = tmp; 2773 work[k] = tmp;
2774 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 2774 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2775 { 2775 {
2776 if (i == mini) 2776 if (i == mini)
2777 continue; 2777 continue;
2778 2778
2779 octave_idx_type iidx = perm[ridx(i)]; 2779 octave_idx_type iidx = perm[ridx (i)];
2780 work[iidx] = work[iidx] - tmp * data(i); 2780 work[iidx] = work[iidx] - tmp * data (i);
2781 } 2781 }
2782 } 2782 }
2783 } 2783 }
2784 2784
2785 for (octave_idx_type i = 0; i < nc; i++) 2785 for (octave_idx_type i = 0; i < nc; i++)
2801 if (work[k] != 0.) 2801 if (work[k] != 0.)
2802 { 2802 {
2803 octave_idx_type minr = nr; 2803 octave_idx_type minr = nr;
2804 octave_idx_type mini = 0; 2804 octave_idx_type mini = 0;
2805 2805
2806 for (octave_idx_type i = cidx(k); 2806 for (octave_idx_type i = cidx (k);
2807 i < cidx(k+1); i++) 2807 i < cidx (k+1); i++)
2808 if (perm[ridx(i)] < minr) 2808 if (perm[ridx (i)] < minr)
2809 { 2809 {
2810 minr = perm[ridx(i)]; 2810 minr = perm[ridx (i)];
2811 mini = i; 2811 mini = i;
2812 } 2812 }
2813 2813
2814 double tmp = work[k] / data(mini); 2814 double tmp = work[k] / data (mini);
2815 work[k] = tmp; 2815 work[k] = tmp;
2816 for (octave_idx_type i = cidx(k); 2816 for (octave_idx_type i = cidx (k);
2817 i < cidx(k+1); i++) 2817 i < cidx (k+1); i++)
2818 { 2818 {
2819 if (i == mini) 2819 if (i == mini)
2820 continue; 2820 continue;
2821 2821
2822 octave_idx_type iidx = perm[ridx(i)]; 2822 octave_idx_type iidx = perm[ridx (i)];
2823 work[iidx] = work[iidx] - tmp * data(i); 2823 work[iidx] = work[iidx] - tmp * data (i);
2824 } 2824 }
2825 } 2825 }
2826 } 2826 }
2827 2827
2828 double atmp = 0; 2828 double atmp = 0;
2829 for (octave_idx_type i = j; i < nc; i++) 2829 for (octave_idx_type i = j; i < nc; i++)
2830 { 2830 {
2831 atmp += fabs(work[i]); 2831 atmp += fabs (work[i]);
2832 work[i] = 0.; 2832 work[i] = 0.;
2833 } 2833 }
2834 if (atmp > ainvnorm) 2834 if (atmp > ainvnorm)
2835 ainvnorm = atmp; 2835 ainvnorm = atmp;
2836 } 2836 }
2850 work[i] = 0.; 2850 work[i] = 0.;
2851 for (octave_idx_type k = 0; k < nc; k++) 2851 for (octave_idx_type k = 0; k < nc; k++)
2852 { 2852 {
2853 if (work[k] != 0.) 2853 if (work[k] != 0.)
2854 { 2854 {
2855 if (ridx(cidx(k)) != k || 2855 if (ridx (cidx (k)) != k ||
2856 data(cidx(k)) == 0.) 2856 data (cidx (k)) == 0.)
2857 { 2857 {
2858 err = -2; 2858 err = -2;
2859 goto triangular_error; 2859 goto triangular_error;
2860 } 2860 }
2861 2861
2862 double tmp = work[k] / data(cidx(k)); 2862 double tmp = work[k] / data (cidx (k));
2863 work[k] = tmp; 2863 work[k] = tmp;
2864 for (octave_idx_type i = cidx(k)+1; 2864 for (octave_idx_type i = cidx (k)+1;
2865 i < cidx(k+1); i++) 2865 i < cidx (k+1); i++)
2866 { 2866 {
2867 octave_idx_type iidx = ridx(i); 2867 octave_idx_type iidx = ridx (i);
2868 work[iidx] = work[iidx] - tmp * data(i); 2868 work[iidx] = work[iidx] - tmp * data (i);
2869 } 2869 }
2870 } 2870 }
2871 } 2871 }
2872 2872
2873 for (octave_idx_type i = 0; i < nc; i++) 2873 for (octave_idx_type i = 0; i < nc; i++)
2887 for (octave_idx_type k = j; k < nc; k++) 2887 for (octave_idx_type k = j; k < nc; k++)
2888 { 2888 {
2889 2889
2890 if (work[k] != 0.) 2890 if (work[k] != 0.)
2891 { 2891 {
2892 double tmp = work[k] / data(cidx(k)); 2892 double tmp = work[k] / data (cidx (k));
2893 work[k] = tmp; 2893 work[k] = tmp;
2894 for (octave_idx_type i = cidx(k)+1; 2894 for (octave_idx_type i = cidx (k)+1;
2895 i < cidx(k+1); i++) 2895 i < cidx (k+1); i++)
2896 { 2896 {
2897 octave_idx_type iidx = ridx(i); 2897 octave_idx_type iidx = ridx (i);
2898 work[iidx] = work[iidx] - tmp * data(i); 2898 work[iidx] = work[iidx] - tmp * data (i);
2899 } 2899 }
2900 } 2900 }
2901 } 2901 }
2902 double atmp = 0; 2902 double atmp = 0;
2903 for (octave_idx_type i = j; i < nc; i++) 2903 for (octave_idx_type i = j; i < nc; i++)
2904 { 2904 {
2905 atmp += fabs(work[i]); 2905 atmp += fabs (work[i]);
2906 work[i] = 0.; 2906 work[i] = 0.;
2907 } 2907 }
2908 if (atmp > ainvnorm) 2908 if (atmp > ainvnorm)
2909 ainvnorm = atmp; 2909 ainvnorm = atmp;
2910 } 2910 }
2985 { 2985 {
2986 // Calculate the 1-norm of matrix for rcond calculation 2986 // Calculate the 1-norm of matrix for rcond calculation
2987 for (octave_idx_type j = 0; j < nc; j++) 2987 for (octave_idx_type j = 0; j < nc; j++)
2988 { 2988 {
2989 double atmp = 0.; 2989 double atmp = 0.;
2990 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 2990 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2991 atmp += fabs(data(i)); 2991 atmp += fabs (data (i));
2992 if (atmp > anorm) 2992 if (atmp > anorm)
2993 anorm = atmp; 2993 anorm = atmp;
2994 } 2994 }
2995 } 2995 }
2996 2996
2997 octave_idx_type b_nc = b.cols (); 2997 octave_idx_type b_nc = b.cols ();
2998 octave_idx_type b_nz = b.nnz (); 2998 octave_idx_type b_nz = b.nnz ();
2999 retval = SparseMatrix (nc, b_nc, b_nz); 2999 retval = SparseMatrix (nc, b_nc, b_nz);
3000 retval.xcidx(0) = 0; 3000 retval.xcidx (0) = 0;
3001 octave_idx_type ii = 0; 3001 octave_idx_type ii = 0;
3002 octave_idx_type x_nz = b_nz; 3002 octave_idx_type x_nz = b_nz;
3003 3003
3004 if (typ == MatrixType::Permuted_Lower) 3004 if (typ == MatrixType::Permuted_Lower)
3005 { 3005 {
3008 3008
3009 for (octave_idx_type j = 0; j < b_nc; j++) 3009 for (octave_idx_type j = 0; j < b_nc; j++)
3010 { 3010 {
3011 for (octave_idx_type i = 0; i < nm; i++) 3011 for (octave_idx_type i = 0; i < nm; i++)
3012 work[i] = 0.; 3012 work[i] = 0.;
3013 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 3013 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3014 work[perm[b.ridx(i)]] = b.data(i); 3014 work[perm[b.ridx (i)]] = b.data (i);
3015 3015
3016 for (octave_idx_type k = 0; k < nc; k++) 3016 for (octave_idx_type k = 0; k < nc; k++)
3017 { 3017 {
3018 if (work[k] != 0.) 3018 if (work[k] != 0.)
3019 { 3019 {
3020 octave_idx_type minr = nr; 3020 octave_idx_type minr = nr;
3021 octave_idx_type mini = 0; 3021 octave_idx_type mini = 0;
3022 3022
3023 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 3023 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3024 if (perm[ridx(i)] < minr) 3024 if (perm[ridx (i)] < minr)
3025 { 3025 {
3026 minr = perm[ridx(i)]; 3026 minr = perm[ridx (i)];
3027 mini = i; 3027 mini = i;
3028 } 3028 }
3029 3029
3030 if (minr != k || data(mini) == 0) 3030 if (minr != k || data (mini) == 0)
3031 { 3031 {
3032 err = -2; 3032 err = -2;
3033 goto triangular_error; 3033 goto triangular_error;
3034 } 3034 }
3035 3035
3036 double tmp = work[k] / data(mini); 3036 double tmp = work[k] / data (mini);
3037 work[k] = tmp; 3037 work[k] = tmp;
3038 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 3038 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3039 { 3039 {
3040 if (i == mini) 3040 if (i == mini)
3041 continue; 3041 continue;
3042 3042
3043 octave_idx_type iidx = perm[ridx(i)]; 3043 octave_idx_type iidx = perm[ridx (i)];
3044 work[iidx] = work[iidx] - tmp * data(i); 3044 work[iidx] = work[iidx] - tmp * data (i);
3045 } 3045 }
3046 } 3046 }
3047 } 3047 }
3048 3048
3049 // Count non-zeros in work vector and adjust space in 3049 // Count non-zeros in work vector and adjust space in
3062 } 3062 }
3063 3063
3064 for (octave_idx_type i = 0; i < nc; i++) 3064 for (octave_idx_type i = 0; i < nc; i++)
3065 if (work[i] != 0.) 3065 if (work[i] != 0.)
3066 { 3066 {
3067 retval.xridx(ii) = i; 3067 retval.xridx (ii) = i;
3068 retval.xdata(ii++) = work[i]; 3068 retval.xdata (ii++) = work[i];
3069 } 3069 }
3070 retval.xcidx(j+1) = ii; 3070 retval.xcidx (j+1) = ii;
3071 } 3071 }
3072 3072
3073 retval.maybe_compress (); 3073 retval.maybe_compress ();
3074 3074
3075 if (calc_cond) 3075 if (calc_cond)
3087 if (work[k] != 0.) 3087 if (work[k] != 0.)
3088 { 3088 {
3089 octave_idx_type minr = nr; 3089 octave_idx_type minr = nr;
3090 octave_idx_type mini = 0; 3090 octave_idx_type mini = 0;
3091 3091
3092 for (octave_idx_type i = cidx(k); 3092 for (octave_idx_type i = cidx (k);
3093 i < cidx(k+1); i++) 3093 i < cidx (k+1); i++)
3094 if (perm[ridx(i)] < minr) 3094 if (perm[ridx (i)] < minr)
3095 { 3095 {
3096 minr = perm[ridx(i)]; 3096 minr = perm[ridx (i)];
3097 mini = i; 3097 mini = i;
3098 } 3098 }
3099 3099
3100 double tmp = work[k] / data(mini); 3100 double tmp = work[k] / data (mini);
3101 work[k] = tmp; 3101 work[k] = tmp;
3102 for (octave_idx_type i = cidx(k); 3102 for (octave_idx_type i = cidx (k);
3103 i < cidx(k+1); i++) 3103 i < cidx (k+1); i++)
3104 { 3104 {
3105 if (i == mini) 3105 if (i == mini)
3106 continue; 3106 continue;
3107 3107
3108 octave_idx_type iidx = perm[ridx(i)]; 3108 octave_idx_type iidx = perm[ridx (i)];
3109 work[iidx] = work[iidx] - tmp * data(i); 3109 work[iidx] = work[iidx] - tmp * data (i);
3110 } 3110 }
3111 } 3111 }
3112 } 3112 }
3113 3113
3114 double atmp = 0; 3114 double atmp = 0;
3115 for (octave_idx_type i = j; i < nr; i++) 3115 for (octave_idx_type i = j; i < nr; i++)
3116 { 3116 {
3117 atmp += fabs(work[i]); 3117 atmp += fabs (work[i]);
3118 work[i] = 0.; 3118 work[i] = 0.;
3119 } 3119 }
3120 if (atmp > ainvnorm) 3120 if (atmp > ainvnorm)
3121 ainvnorm = atmp; 3121 ainvnorm = atmp;
3122 } 3122 }
3129 3129
3130 for (octave_idx_type j = 0; j < b_nc; j++) 3130 for (octave_idx_type j = 0; j < b_nc; j++)
3131 { 3131 {
3132 for (octave_idx_type i = 0; i < nm; i++) 3132 for (octave_idx_type i = 0; i < nm; i++)
3133 work[i] = 0.; 3133 work[i] = 0.;
3134 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 3134 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3135 work[b.ridx(i)] = b.data(i); 3135 work[b.ridx (i)] = b.data (i);
3136 3136
3137 for (octave_idx_type k = 0; k < nc; k++) 3137 for (octave_idx_type k = 0; k < nc; k++)
3138 { 3138 {
3139 if (work[k] != 0.) 3139 if (work[k] != 0.)
3140 { 3140 {
3141 if (ridx(cidx(k)) != k || 3141 if (ridx (cidx (k)) != k ||
3142 data(cidx(k)) == 0.) 3142 data (cidx (k)) == 0.)
3143 { 3143 {
3144 err = -2; 3144 err = -2;
3145 goto triangular_error; 3145 goto triangular_error;
3146 } 3146 }
3147 3147
3148 double tmp = work[k] / data(cidx(k)); 3148 double tmp = work[k] / data (cidx (k));
3149 work[k] = tmp; 3149 work[k] = tmp;
3150 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) 3150 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3151 { 3151 {
3152 octave_idx_type iidx = ridx(i); 3152 octave_idx_type iidx = ridx (i);
3153 work[iidx] = work[iidx] - tmp * data(i); 3153 work[iidx] = work[iidx] - tmp * data (i);
3154 } 3154 }
3155 } 3155 }
3156 } 3156 }
3157 3157
3158 // Count non-zeros in work vector and adjust space in 3158 // Count non-zeros in work vector and adjust space in
3171 } 3171 }
3172 3172
3173 for (octave_idx_type i = 0; i < nc; i++) 3173 for (octave_idx_type i = 0; i < nc; i++)
3174 if (work[i] != 0.) 3174 if (work[i] != 0.)
3175 { 3175 {
3176 retval.xridx(ii) = i; 3176 retval.xridx (ii) = i;
3177 retval.xdata(ii++) = work[i]; 3177 retval.xdata (ii++) = work[i];
3178 } 3178 }
3179 retval.xcidx(j+1) = ii; 3179 retval.xcidx (j+1) = ii;
3180 } 3180 }
3181 3181
3182 retval.maybe_compress (); 3182 retval.maybe_compress ();
3183 3183
3184 if (calc_cond) 3184 if (calc_cond)
3194 for (octave_idx_type k = j; k < nc; k++) 3194 for (octave_idx_type k = j; k < nc; k++)
3195 { 3195 {
3196 3196
3197 if (work[k] != 0.) 3197 if (work[k] != 0.)
3198 { 3198 {
3199 double tmp = work[k] / data(cidx(k)); 3199 double tmp = work[k] / data (cidx (k));
3200 work[k] = tmp; 3200 work[k] = tmp;
3201 for (octave_idx_type i = cidx(k)+1; 3201 for (octave_idx_type i = cidx (k)+1;
3202 i < cidx(k+1); i++) 3202 i < cidx (k+1); i++)
3203 { 3203 {
3204 octave_idx_type iidx = ridx(i); 3204 octave_idx_type iidx = ridx (i);
3205 work[iidx] = work[iidx] - tmp * data(i); 3205 work[iidx] = work[iidx] - tmp * data (i);
3206 } 3206 }
3207 } 3207 }
3208 } 3208 }
3209 double atmp = 0; 3209 double atmp = 0;
3210 for (octave_idx_type i = j; i < nc; i++) 3210 for (octave_idx_type i = j; i < nc; i++)
3211 { 3211 {
3212 atmp += fabs(work[i]); 3212 atmp += fabs (work[i]);
3213 work[i] = 0.; 3213 work[i] = 0.;
3214 } 3214 }
3215 if (atmp > ainvnorm) 3215 if (atmp > ainvnorm)
3216 ainvnorm = atmp; 3216 ainvnorm = atmp;
3217 } 3217 }
3293 { 3293 {
3294 // Calculate the 1-norm of matrix for rcond calculation 3294 // Calculate the 1-norm of matrix for rcond calculation
3295 for (octave_idx_type j = 0; j < nc; j++) 3295 for (octave_idx_type j = 0; j < nc; j++)
3296 { 3296 {
3297 double atmp = 0.; 3297 double atmp = 0.;
3298 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 3298 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3299 atmp += fabs(data(i)); 3299 atmp += fabs (data (i));
3300 if (atmp > anorm) 3300 if (atmp > anorm)
3301 anorm = atmp; 3301 anorm = atmp;
3302 } 3302 }
3303 } 3303 }
3304 3304
3320 if (cwork[k] != 0.) 3320 if (cwork[k] != 0.)
3321 { 3321 {
3322 octave_idx_type minr = nr; 3322 octave_idx_type minr = nr;
3323 octave_idx_type mini = 0; 3323 octave_idx_type mini = 0;
3324 3324
3325 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 3325 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3326 if (perm[ridx(i)] < minr) 3326 if (perm[ridx (i)] < minr)
3327 { 3327 {
3328 minr = perm[ridx(i)]; 3328 minr = perm[ridx (i)];
3329 mini = i; 3329 mini = i;
3330 } 3330 }
3331 3331
3332 if (minr != k || data(mini) == 0) 3332 if (minr != k || data (mini) == 0)
3333 { 3333 {
3334 err = -2; 3334 err = -2;
3335 goto triangular_error; 3335 goto triangular_error;
3336 } 3336 }
3337 3337
3338 Complex tmp = cwork[k] / data(mini); 3338 Complex tmp = cwork[k] / data (mini);
3339 cwork[k] = tmp; 3339 cwork[k] = tmp;
3340 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 3340 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3341 { 3341 {
3342 if (i == mini) 3342 if (i == mini)
3343 continue; 3343 continue;
3344 3344
3345 octave_idx_type iidx = perm[ridx(i)]; 3345 octave_idx_type iidx = perm[ridx (i)];
3346 cwork[iidx] = cwork[iidx] - tmp * data(i); 3346 cwork[iidx] = cwork[iidx] - tmp * data (i);
3347 } 3347 }
3348 } 3348 }
3349 } 3349 }
3350 3350
3351 for (octave_idx_type i = 0; i < nc; i++) 3351 for (octave_idx_type i = 0; i < nc; i++)
3368 if (work[k] != 0.) 3368 if (work[k] != 0.)
3369 { 3369 {
3370 octave_idx_type minr = nr; 3370 octave_idx_type minr = nr;
3371 octave_idx_type mini = 0; 3371 octave_idx_type mini = 0;
3372 3372
3373 for (octave_idx_type i = cidx(k); 3373 for (octave_idx_type i = cidx (k);
3374 i < cidx(k+1); i++) 3374 i < cidx (k+1); i++)
3375 if (perm[ridx(i)] < minr) 3375 if (perm[ridx (i)] < minr)
3376 { 3376 {
3377 minr = perm[ridx(i)]; 3377 minr = perm[ridx (i)];
3378 mini = i; 3378 mini = i;
3379 } 3379 }
3380 3380
3381 double tmp = work[k] / data(mini); 3381 double tmp = work[k] / data (mini);
3382 work[k] = tmp; 3382 work[k] = tmp;
3383 for (octave_idx_type i = cidx(k); 3383 for (octave_idx_type i = cidx (k);
3384 i < cidx(k+1); i++) 3384 i < cidx (k+1); i++)
3385 { 3385 {
3386 if (i == mini) 3386 if (i == mini)
3387 continue; 3387 continue;
3388 3388
3389 octave_idx_type iidx = perm[ridx(i)]; 3389 octave_idx_type iidx = perm[ridx (i)];
3390 work[iidx] = work[iidx] - tmp * data(i); 3390 work[iidx] = work[iidx] - tmp * data (i);
3391 } 3391 }
3392 } 3392 }
3393 } 3393 }
3394 3394
3395 double atmp = 0; 3395 double atmp = 0;
3396 for (octave_idx_type i = j; i < nc; i++) 3396 for (octave_idx_type i = j; i < nc; i++)
3397 { 3397 {
3398 atmp += fabs(work[i]); 3398 atmp += fabs (work[i]);
3399 work[i] = 0.; 3399 work[i] = 0.;
3400 } 3400 }
3401 if (atmp > ainvnorm) 3401 if (atmp > ainvnorm)
3402 ainvnorm = atmp; 3402 ainvnorm = atmp;
3403 } 3403 }
3418 3418
3419 for (octave_idx_type k = 0; k < nc; k++) 3419 for (octave_idx_type k = 0; k < nc; k++)
3420 { 3420 {
3421 if (cwork[k] != 0.) 3421 if (cwork[k] != 0.)
3422 { 3422 {
3423 if (ridx(cidx(k)) != k || 3423 if (ridx (cidx (k)) != k ||
3424 data(cidx(k)) == 0.) 3424 data (cidx (k)) == 0.)
3425 { 3425 {
3426 err = -2; 3426 err = -2;
3427 goto triangular_error; 3427 goto triangular_error;
3428 } 3428 }
3429 3429
3430 Complex tmp = cwork[k] / data(cidx(k)); 3430 Complex tmp = cwork[k] / data (cidx (k));
3431 cwork[k] = tmp; 3431 cwork[k] = tmp;
3432 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) 3432 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3433 { 3433 {
3434 octave_idx_type iidx = ridx(i); 3434 octave_idx_type iidx = ridx (i);
3435 cwork[iidx] = cwork[iidx] - tmp * data(i); 3435 cwork[iidx] = cwork[iidx] - tmp * data (i);
3436 } 3436 }
3437 } 3437 }
3438 } 3438 }
3439 3439
3440 for (octave_idx_type i = 0; i < nc; i++) 3440 for (octave_idx_type i = 0; i < nc; i++)
3455 for (octave_idx_type k = j; k < nc; k++) 3455 for (octave_idx_type k = j; k < nc; k++)
3456 { 3456 {
3457 3457
3458 if (work[k] != 0.) 3458 if (work[k] != 0.)
3459 { 3459 {
3460 double tmp = work[k] / data(cidx(k)); 3460 double tmp = work[k] / data (cidx (k));
3461 work[k] = tmp; 3461 work[k] = tmp;
3462 for (octave_idx_type i = cidx(k)+1; 3462 for (octave_idx_type i = cidx (k)+1;
3463 i < cidx(k+1); i++) 3463 i < cidx (k+1); i++)
3464 { 3464 {
3465 octave_idx_type iidx = ridx(i); 3465 octave_idx_type iidx = ridx (i);
3466 work[iidx] = work[iidx] - tmp * data(i); 3466 work[iidx] = work[iidx] - tmp * data (i);
3467 } 3467 }
3468 } 3468 }
3469 } 3469 }
3470 double atmp = 0; 3470 double atmp = 0;
3471 for (octave_idx_type i = j; i < nc; i++) 3471 for (octave_idx_type i = j; i < nc; i++)
3472 { 3472 {
3473 atmp += fabs(work[i]); 3473 atmp += fabs (work[i]);
3474 work[i] = 0.; 3474 work[i] = 0.;
3475 } 3475 }
3476 if (atmp > ainvnorm) 3476 if (atmp > ainvnorm)
3477 ainvnorm = atmp; 3477 ainvnorm = atmp;
3478 } 3478 }
3553 { 3553 {
3554 // Calculate the 1-norm of matrix for rcond calculation 3554 // Calculate the 1-norm of matrix for rcond calculation
3555 for (octave_idx_type j = 0; j < nc; j++) 3555 for (octave_idx_type j = 0; j < nc; j++)
3556 { 3556 {
3557 double atmp = 0.; 3557 double atmp = 0.;
3558 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 3558 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3559 atmp += fabs(data(i)); 3559 atmp += fabs (data (i));
3560 if (atmp > anorm) 3560 if (atmp > anorm)
3561 anorm = atmp; 3561 anorm = atmp;
3562 } 3562 }
3563 } 3563 }
3564 3564
3565 octave_idx_type b_nc = b.cols (); 3565 octave_idx_type b_nc = b.cols ();
3566 octave_idx_type b_nz = b.nnz (); 3566 octave_idx_type b_nz = b.nnz ();
3567 retval = SparseComplexMatrix (nc, b_nc, b_nz); 3567 retval = SparseComplexMatrix (nc, b_nc, b_nz);
3568 retval.xcidx(0) = 0; 3568 retval.xcidx (0) = 0;
3569 octave_idx_type ii = 0; 3569 octave_idx_type ii = 0;
3570 octave_idx_type x_nz = b_nz; 3570 octave_idx_type x_nz = b_nz;
3571 3571
3572 if (typ == MatrixType::Permuted_Lower) 3572 if (typ == MatrixType::Permuted_Lower)
3573 { 3573 {
3576 3576
3577 for (octave_idx_type j = 0; j < b_nc; j++) 3577 for (octave_idx_type j = 0; j < b_nc; j++)
3578 { 3578 {
3579 for (octave_idx_type i = 0; i < nm; i++) 3579 for (octave_idx_type i = 0; i < nm; i++)
3580 cwork[i] = 0.; 3580 cwork[i] = 0.;
3581 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 3581 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3582 cwork[perm[b.ridx(i)]] = b.data(i); 3582 cwork[perm[b.ridx (i)]] = b.data (i);
3583 3583
3584 for (octave_idx_type k = 0; k < nc; k++) 3584 for (octave_idx_type k = 0; k < nc; k++)
3585 { 3585 {
3586 if (cwork[k] != 0.) 3586 if (cwork[k] != 0.)
3587 { 3587 {
3588 octave_idx_type minr = nr; 3588 octave_idx_type minr = nr;
3589 octave_idx_type mini = 0; 3589 octave_idx_type mini = 0;
3590 3590
3591 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 3591 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3592 if (perm[ridx(i)] < minr) 3592 if (perm[ridx (i)] < minr)
3593 { 3593 {
3594 minr = perm[ridx(i)]; 3594 minr = perm[ridx (i)];
3595 mini = i; 3595 mini = i;
3596 } 3596 }
3597 3597
3598 if (minr != k || data(mini) == 0) 3598 if (minr != k || data (mini) == 0)
3599 { 3599 {
3600 err = -2; 3600 err = -2;
3601 goto triangular_error; 3601 goto triangular_error;
3602 } 3602 }
3603 3603
3604 Complex tmp = cwork[k] / data(mini); 3604 Complex tmp = cwork[k] / data (mini);
3605 cwork[k] = tmp; 3605 cwork[k] = tmp;
3606 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) 3606 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3607 { 3607 {
3608 if (i == mini) 3608 if (i == mini)
3609 continue; 3609 continue;
3610 3610
3611 octave_idx_type iidx = perm[ridx(i)]; 3611 octave_idx_type iidx = perm[ridx (i)];
3612 cwork[iidx] = cwork[iidx] - tmp * data(i); 3612 cwork[iidx] = cwork[iidx] - tmp * data (i);
3613 } 3613 }
3614 } 3614 }
3615 } 3615 }
3616 3616
3617 // Count non-zeros in work vector and adjust space in 3617 // Count non-zeros in work vector and adjust space in
3630 } 3630 }
3631 3631
3632 for (octave_idx_type i = 0; i < nc; i++) 3632 for (octave_idx_type i = 0; i < nc; i++)
3633 if (cwork[i] != 0.) 3633 if (cwork[i] != 0.)
3634 { 3634 {
3635 retval.xridx(ii) = i; 3635 retval.xridx (ii) = i;
3636 retval.xdata(ii++) = cwork[i]; 3636 retval.xdata (ii++) = cwork[i];
3637 } 3637 }
3638 retval.xcidx(j+1) = ii; 3638 retval.xcidx (j+1) = ii;
3639 } 3639 }
3640 3640
3641 retval.maybe_compress (); 3641 retval.maybe_compress ();
3642 3642
3643 if (calc_cond) 3643 if (calc_cond)
3656 if (work[k] != 0.) 3656 if (work[k] != 0.)
3657 { 3657 {
3658 octave_idx_type minr = nr; 3658 octave_idx_type minr = nr;
3659 octave_idx_type mini = 0; 3659 octave_idx_type mini = 0;
3660 3660
3661 for (octave_idx_type i = cidx(k); 3661 for (octave_idx_type i = cidx (k);
3662 i < cidx(k+1); i++) 3662 i < cidx (k+1); i++)
3663 if (perm[ridx(i)] < minr) 3663 if (perm[ridx (i)] < minr)
3664 { 3664 {
3665 minr = perm[ridx(i)]; 3665 minr = perm[ridx (i)];
3666 mini = i; 3666 mini = i;
3667 } 3667 }
3668 3668
3669 double tmp = work[k] / data(mini); 3669 double tmp = work[k] / data (mini);
3670 work[k] = tmp; 3670 work[k] = tmp;
3671 for (octave_idx_type i = cidx(k); 3671 for (octave_idx_type i = cidx (k);
3672 i < cidx(k+1); i++) 3672 i < cidx (k+1); i++)
3673 { 3673 {
3674 if (i == mini) 3674 if (i == mini)
3675 continue; 3675 continue;
3676 3676
3677 octave_idx_type iidx = perm[ridx(i)]; 3677 octave_idx_type iidx = perm[ridx (i)];
3678 work[iidx] = work[iidx] - tmp * data(i); 3678 work[iidx] = work[iidx] - tmp * data (i);
3679 } 3679 }
3680 } 3680 }
3681 } 3681 }
3682 3682
3683 double atmp = 0; 3683 double atmp = 0;
3684 for (octave_idx_type i = j; i < nc; i++) 3684 for (octave_idx_type i = j; i < nc; i++)
3685 { 3685 {
3686 atmp += fabs(work[i]); 3686 atmp += fabs (work[i]);
3687 work[i] = 0.; 3687 work[i] = 0.;
3688 } 3688 }
3689 if (atmp > ainvnorm) 3689 if (atmp > ainvnorm)
3690 ainvnorm = atmp; 3690 ainvnorm = atmp;
3691 } 3691 }
3698 3698
3699 for (octave_idx_type j = 0; j < b_nc; j++) 3699 for (octave_idx_type j = 0; j < b_nc; j++)
3700 { 3700 {
3701 for (octave_idx_type i = 0; i < nm; i++) 3701 for (octave_idx_type i = 0; i < nm; i++)
3702 cwork[i] = 0.; 3702 cwork[i] = 0.;
3703 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 3703 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3704 cwork[b.ridx(i)] = b.data(i); 3704 cwork[b.ridx (i)] = b.data (i);
3705 3705
3706 for (octave_idx_type k = 0; k < nc; k++) 3706 for (octave_idx_type k = 0; k < nc; k++)
3707 { 3707 {
3708 if (cwork[k] != 0.) 3708 if (cwork[k] != 0.)
3709 { 3709 {
3710 if (ridx(cidx(k)) != k || 3710 if (ridx (cidx (k)) != k ||
3711 data(cidx(k)) == 0.) 3711 data (cidx (k)) == 0.)
3712 { 3712 {
3713 err = -2; 3713 err = -2;
3714 goto triangular_error; 3714 goto triangular_error;
3715 } 3715 }
3716 3716
3717 Complex tmp = cwork[k] / data(cidx(k)); 3717 Complex tmp = cwork[k] / data (cidx (k));
3718 cwork[k] = tmp; 3718 cwork[k] = tmp;
3719 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) 3719 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3720 { 3720 {
3721 octave_idx_type iidx = ridx(i); 3721 octave_idx_type iidx = ridx (i);
3722 cwork[iidx] = cwork[iidx] - tmp * data(i); 3722 cwork[iidx] = cwork[iidx] - tmp * data (i);
3723 } 3723 }
3724 } 3724 }
3725 } 3725 }
3726 3726
3727 // Count non-zeros in work vector and adjust space in 3727 // Count non-zeros in work vector and adjust space in
3740 } 3740 }
3741 3741
3742 for (octave_idx_type i = 0; i < nc; i++) 3742 for (octave_idx_type i = 0; i < nc; i++)
3743 if (cwork[i] != 0.) 3743 if (cwork[i] != 0.)
3744 { 3744 {
3745 retval.xridx(ii) = i; 3745 retval.xridx (ii) = i;
3746 retval.xdata(ii++) = cwork[i]; 3746 retval.xdata (ii++) = cwork[i];
3747 } 3747 }
3748 retval.xcidx(j+1) = ii; 3748 retval.xcidx (j+1) = ii;
3749 } 3749 }
3750 3750
3751 retval.maybe_compress (); 3751 retval.maybe_compress ();
3752 3752
3753 if (calc_cond) 3753 if (calc_cond)
3764 for (octave_idx_type k = j; k < nc; k++) 3764 for (octave_idx_type k = j; k < nc; k++)
3765 { 3765 {
3766 3766
3767 if (work[k] != 0.) 3767 if (work[k] != 0.)
3768 { 3768 {
3769 double tmp = work[k] / data(cidx(k)); 3769 double tmp = work[k] / data (cidx (k));
3770 work[k] = tmp; 3770 work[k] = tmp;
3771 for (octave_idx_type i = cidx(k)+1; 3771 for (octave_idx_type i = cidx (k)+1;
3772 i < cidx(k+1); i++) 3772 i < cidx (k+1); i++)
3773 { 3773 {
3774 octave_idx_type iidx = ridx(i); 3774 octave_idx_type iidx = ridx (i);
3775 work[iidx] = work[iidx] - tmp * data(i); 3775 work[iidx] = work[iidx] - tmp * data (i);
3776 } 3776 }
3777 } 3777 }
3778 } 3778 }
3779 double atmp = 0; 3779 double atmp = 0;
3780 for (octave_idx_type i = j; i < nc; i++) 3780 for (octave_idx_type i = j; i < nc; i++)
3781 { 3781 {
3782 atmp += fabs(work[i]); 3782 atmp += fabs (work[i]);
3783 work[i] = 0.; 3783 work[i] = 0.;
3784 } 3784 }
3785 if (atmp > ainvnorm) 3785 if (atmp > ainvnorm)
3786 ainvnorm = atmp; 3786 ainvnorm = atmp;
3787 } 3787 }
3862 { 3862 {
3863 octave_idx_type ii = 0; 3863 octave_idx_type ii = 0;
3864 3864
3865 for (octave_idx_type j = 0; j < nc-1; j++) 3865 for (octave_idx_type j = 0; j < nc-1; j++)
3866 { 3866 {
3867 D[j] = data(ii++); 3867 D[j] = data (ii++);
3868 DL[j] = data(ii); 3868 DL[j] = data (ii);
3869 ii += 2; 3869 ii += 2;
3870 } 3870 }
3871 D[nc-1] = data(ii); 3871 D[nc-1] = data (ii);
3872 } 3872 }
3873 else 3873 else
3874 { 3874 {
3875 D[0] = 0.; 3875 D[0] = 0.;
3876 for (octave_idx_type i = 0; i < nr - 1; i++) 3876 for (octave_idx_type i = 0; i < nr - 1; i++)
3878 D[i+1] = 0.; 3878 D[i+1] = 0.;
3879 DL[i] = 0.; 3879 DL[i] = 0.;
3880 } 3880 }
3881 3881
3882 for (octave_idx_type j = 0; j < nc; j++) 3882 for (octave_idx_type j = 0; j < nc; j++)
3883 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 3883 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3884 { 3884 {
3885 if (ridx(i) == j) 3885 if (ridx (i) == j)
3886 D[j] = data(i); 3886 D[j] = data (i);
3887 else if (ridx(i) == j + 1) 3887 else if (ridx (i) == j + 1)
3888 DL[j] = data(i); 3888 DL[j] = data (i);
3889 } 3889 }
3890 } 3890 }
3891 3891
3892 octave_idx_type b_nc = b.cols (); 3892 octave_idx_type b_nc = b.cols ();
3893 retval = b; 3893 retval = b;
3916 { 3916 {
3917 octave_idx_type ii = 0; 3917 octave_idx_type ii = 0;
3918 3918
3919 for (octave_idx_type j = 0; j < nc-1; j++) 3919 for (octave_idx_type j = 0; j < nc-1; j++)
3920 { 3920 {
3921 D[j] = data(ii++); 3921 D[j] = data (ii++);
3922 DL[j] = data(ii++); 3922 DL[j] = data (ii++);
3923 DU[j] = data(ii++); 3923 DU[j] = data (ii++);
3924 } 3924 }
3925 D[nc-1] = data(ii); 3925 D[nc-1] = data (ii);
3926 } 3926 }
3927 else 3927 else
3928 { 3928 {
3929 D[0] = 0.; 3929 D[0] = 0.;
3930 for (octave_idx_type i = 0; i < nr - 1; i++) 3930 for (octave_idx_type i = 0; i < nr - 1; i++)
3933 DL[i] = 0.; 3933 DL[i] = 0.;
3934 DU[i] = 0.; 3934 DU[i] = 0.;
3935 } 3935 }
3936 3936
3937 for (octave_idx_type j = 0; j < nc; j++) 3937 for (octave_idx_type j = 0; j < nc; j++)
3938 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 3938 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3939 { 3939 {
3940 if (ridx(i) == j) 3940 if (ridx (i) == j)
3941 D[j] = data(i); 3941 D[j] = data (i);
3942 else if (ridx(i) == j + 1) 3942 else if (ridx (i) == j + 1)
3943 DL[j] = data(i); 3943 DL[j] = data (i);
3944 else if (ridx(i) == j - 1) 3944 else if (ridx (i) == j - 1)
3945 DU[j-1] = data(i); 3945 DU[j-1] = data (i);
3946 } 3946 }
3947 } 3947 }
3948 3948
3949 octave_idx_type b_nc = b.cols (); 3949 octave_idx_type b_nc = b.cols ();
3950 retval = b; 3950 retval = b;
4019 { 4019 {
4020 octave_idx_type ii = 0; 4020 octave_idx_type ii = 0;
4021 4021
4022 for (octave_idx_type j = 0; j < nc-1; j++) 4022 for (octave_idx_type j = 0; j < nc-1; j++)
4023 { 4023 {
4024 D[j] = data(ii++); 4024 D[j] = data (ii++);
4025 DL[j] = data(ii++); 4025 DL[j] = data (ii++);
4026 DU[j] = data(ii++); 4026 DU[j] = data (ii++);
4027 } 4027 }
4028 D[nc-1] = data(ii); 4028 D[nc-1] = data (ii);
4029 } 4029 }
4030 else 4030 else
4031 { 4031 {
4032 D[0] = 0.; 4032 D[0] = 0.;
4033 for (octave_idx_type i = 0; i < nr - 1; i++) 4033 for (octave_idx_type i = 0; i < nr - 1; i++)
4036 DL[i] = 0.; 4036 DL[i] = 0.;
4037 DU[i] = 0.; 4037 DU[i] = 0.;
4038 } 4038 }
4039 4039
4040 for (octave_idx_type j = 0; j < nc; j++) 4040 for (octave_idx_type j = 0; j < nc; j++)
4041 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4041 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4042 { 4042 {
4043 if (ridx(i) == j) 4043 if (ridx (i) == j)
4044 D[j] = data(i); 4044 D[j] = data (i);
4045 else if (ridx(i) == j + 1) 4045 else if (ridx (i) == j + 1)
4046 DL[j] = data(i); 4046 DL[j] = data (i);
4047 else if (ridx(i) == j - 1) 4047 else if (ridx (i) == j - 1)
4048 DU[j-1] = data(i); 4048 DU[j-1] = data (i);
4049 } 4049 }
4050 } 4050 }
4051 4051
4052 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); 4052 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4053 4053
4071 rcond = 1.0; 4071 rcond = 1.0;
4072 char job = 'N'; 4072 char job = 'N';
4073 volatile octave_idx_type x_nz = b.nnz (); 4073 volatile octave_idx_type x_nz = b.nnz ();
4074 octave_idx_type b_nc = b.cols (); 4074 octave_idx_type b_nc = b.cols ();
4075 retval = SparseMatrix (nr, b_nc, x_nz); 4075 retval = SparseMatrix (nr, b_nc, x_nz);
4076 retval.xcidx(0) = 0; 4076 retval.xcidx (0) = 0;
4077 volatile octave_idx_type ii = 0; 4077 volatile octave_idx_type ii = 0;
4078 4078
4079 OCTAVE_LOCAL_BUFFER (double, work, nr); 4079 OCTAVE_LOCAL_BUFFER (double, work, nr);
4080 4080
4081 for (volatile octave_idx_type j = 0; j < b_nc; j++) 4081 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4082 { 4082 {
4083 for (octave_idx_type i = 0; i < nr; i++) 4083 for (octave_idx_type i = 0; i < nr; i++)
4084 work[i] = 0.; 4084 work[i] = 0.;
4085 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 4085 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
4086 work[b.ridx(i)] = b.data(i); 4086 work[b.ridx (i)] = b.data (i);
4087 4087
4088 F77_XFCN (dgttrs, DGTTRS, 4088 F77_XFCN (dgttrs, DGTTRS,
4089 (F77_CONST_CHAR_ARG2 (&job, 1), 4089 (F77_CONST_CHAR_ARG2 (&job, 1),
4090 nr, 1, DL, D, DU, DU2, pipvt, 4090 nr, 1, DL, D, DU, DU2, pipvt,
4091 work, b.rows (), err 4091 work, b.rows (), err
4107 } 4107 }
4108 4108
4109 for (octave_idx_type i = 0; i < nr; i++) 4109 for (octave_idx_type i = 0; i < nr; i++)
4110 if (work[i] != 0.) 4110 if (work[i] != 0.)
4111 { 4111 {
4112 retval.xridx(ii) = i; 4112 retval.xridx (ii) = i;
4113 retval.xdata(ii++) = work[i]; 4113 retval.xdata (ii++) = work[i];
4114 } 4114 }
4115 retval.xcidx(j+1) = ii; 4115 retval.xcidx (j+1) = ii;
4116 } 4116 }
4117 4117
4118 retval.maybe_compress (); 4118 retval.maybe_compress ();
4119 } 4119 }
4120 } 4120 }
4160 { 4160 {
4161 octave_idx_type ii = 0; 4161 octave_idx_type ii = 0;
4162 4162
4163 for (octave_idx_type j = 0; j < nc-1; j++) 4163 for (octave_idx_type j = 0; j < nc-1; j++)
4164 { 4164 {
4165 D[j] = data(ii++); 4165 D[j] = data (ii++);
4166 DL[j] = data(ii); 4166 DL[j] = data (ii);
4167 ii += 2; 4167 ii += 2;
4168 } 4168 }
4169 D[nc-1] = data(ii); 4169 D[nc-1] = data (ii);
4170 } 4170 }
4171 else 4171 else
4172 { 4172 {
4173 D[0] = 0.; 4173 D[0] = 0.;
4174 for (octave_idx_type i = 0; i < nr - 1; i++) 4174 for (octave_idx_type i = 0; i < nr - 1; i++)
4176 D[i+1] = 0.; 4176 D[i+1] = 0.;
4177 DL[i] = 0.; 4177 DL[i] = 0.;
4178 } 4178 }
4179 4179
4180 for (octave_idx_type j = 0; j < nc; j++) 4180 for (octave_idx_type j = 0; j < nc; j++)
4181 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4181 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4182 { 4182 {
4183 if (ridx(i) == j) 4183 if (ridx (i) == j)
4184 D[j] = data(i); 4184 D[j] = data (i);
4185 else if (ridx(i) == j + 1) 4185 else if (ridx (i) == j + 1)
4186 DL[j] = data(i); 4186 DL[j] = data (i);
4187 } 4187 }
4188 } 4188 }
4189 4189
4190 octave_idx_type b_nr = b.rows (); 4190 octave_idx_type b_nr = b.rows ();
4191 octave_idx_type b_nc = b.cols (); 4191 octave_idx_type b_nc = b.cols ();
4215 { 4215 {
4216 octave_idx_type ii = 0; 4216 octave_idx_type ii = 0;
4217 4217
4218 for (octave_idx_type j = 0; j < nc-1; j++) 4218 for (octave_idx_type j = 0; j < nc-1; j++)
4219 { 4219 {
4220 D[j] = data(ii++); 4220 D[j] = data (ii++);
4221 DL[j] = data(ii++); 4221 DL[j] = data (ii++);
4222 DU[j] = data(ii++); 4222 DU[j] = data (ii++);
4223 } 4223 }
4224 D[nc-1] = data(ii); 4224 D[nc-1] = data (ii);
4225 } 4225 }
4226 else 4226 else
4227 { 4227 {
4228 D[0] = 0.; 4228 D[0] = 0.;
4229 for (octave_idx_type i = 0; i < nr - 1; i++) 4229 for (octave_idx_type i = 0; i < nr - 1; i++)
4232 DL[i] = 0.; 4232 DL[i] = 0.;
4233 DU[i] = 0.; 4233 DU[i] = 0.;
4234 } 4234 }
4235 4235
4236 for (octave_idx_type j = 0; j < nc; j++) 4236 for (octave_idx_type j = 0; j < nc; j++)
4237 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4237 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4238 { 4238 {
4239 if (ridx(i) == j) 4239 if (ridx (i) == j)
4240 D[j] = data(i); 4240 D[j] = data (i);
4241 else if (ridx(i) == j + 1) 4241 else if (ridx (i) == j + 1)
4242 DL[j] = data(i); 4242 DL[j] = data (i);
4243 else if (ridx(i) == j - 1) 4243 else if (ridx (i) == j - 1)
4244 DU[j-1] = data(i); 4244 DU[j-1] = data (i);
4245 } 4245 }
4246 } 4246 }
4247 4247
4248 octave_idx_type b_nr = b.rows (); 4248 octave_idx_type b_nr = b.rows ();
4249 octave_idx_type b_nc = b.cols (); 4249 octave_idx_type b_nc = b.cols ();
4318 { 4318 {
4319 octave_idx_type ii = 0; 4319 octave_idx_type ii = 0;
4320 4320
4321 for (octave_idx_type j = 0; j < nc-1; j++) 4321 for (octave_idx_type j = 0; j < nc-1; j++)
4322 { 4322 {
4323 D[j] = data(ii++); 4323 D[j] = data (ii++);
4324 DL[j] = data(ii++); 4324 DL[j] = data (ii++);
4325 DU[j] = data(ii++); 4325 DU[j] = data (ii++);
4326 } 4326 }
4327 D[nc-1] = data(ii); 4327 D[nc-1] = data (ii);
4328 } 4328 }
4329 else 4329 else
4330 { 4330 {
4331 D[0] = 0.; 4331 D[0] = 0.;
4332 for (octave_idx_type i = 0; i < nr - 1; i++) 4332 for (octave_idx_type i = 0; i < nr - 1; i++)
4335 DL[i] = 0.; 4335 DL[i] = 0.;
4336 DU[i] = 0.; 4336 DU[i] = 0.;
4337 } 4337 }
4338 4338
4339 for (octave_idx_type j = 0; j < nc; j++) 4339 for (octave_idx_type j = 0; j < nc; j++)
4340 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4340 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4341 { 4341 {
4342 if (ridx(i) == j) 4342 if (ridx (i) == j)
4343 D[j] = data(i); 4343 D[j] = data (i);
4344 else if (ridx(i) == j + 1) 4344 else if (ridx (i) == j + 1)
4345 DL[j] = data(i); 4345 DL[j] = data (i);
4346 else if (ridx(i) == j - 1) 4346 else if (ridx (i) == j - 1)
4347 DU[j-1] = data(i); 4347 DU[j-1] = data (i);
4348 } 4348 }
4349 } 4349 }
4350 4350
4351 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); 4351 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4352 4352
4377 // will be as many as in b 4377 // will be as many as in b
4378 volatile octave_idx_type x_nz = b.nnz (); 4378 volatile octave_idx_type x_nz = b.nnz ();
4379 volatile octave_idx_type ii = 0; 4379 volatile octave_idx_type ii = 0;
4380 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); 4380 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4381 4381
4382 retval.xcidx(0) = 0; 4382 retval.xcidx (0) = 0;
4383 for (volatile octave_idx_type j = 0; j < b_nc; j++) 4383 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4384 { 4384 {
4385 4385
4386 for (octave_idx_type i = 0; i < b_nr; i++) 4386 for (octave_idx_type i = 0; i < b_nr; i++)
4387 { 4387 {
4436 } 4436 }
4437 4437
4438 for (octave_idx_type i = 0; i < nr; i++) 4438 for (octave_idx_type i = 0; i < nr; i++)
4439 if (Bx[i] != 0. || Bz[i] != 0.) 4439 if (Bx[i] != 0. || Bz[i] != 0.)
4440 { 4440 {
4441 retval.xridx(ii) = i; 4441 retval.xridx (ii) = i;
4442 retval.xdata(ii++) = 4442 retval.xdata (ii++) =
4443 Complex (Bx[i], Bz[i]); 4443 Complex (Bx[i], Bz[i]);
4444 } 4444 }
4445 4445
4446 retval.xcidx(j+1) = ii; 4446 retval.xcidx (j+1) = ii;
4447 } 4447 }
4448 4448
4449 retval.maybe_compress (); 4449 retval.maybe_compress ();
4450 } 4450 }
4451 } 4451 }
4494 for (octave_idx_type i = 0; i < nc; i++) 4494 for (octave_idx_type i = 0; i < nc; i++)
4495 tmp_data[ii++] = 0.; 4495 tmp_data[ii++] = 0.;
4496 } 4496 }
4497 4497
4498 for (octave_idx_type j = 0; j < nc; j++) 4498 for (octave_idx_type j = 0; j < nc; j++)
4499 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4499 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4500 { 4500 {
4501 octave_idx_type ri = ridx (i); 4501 octave_idx_type ri = ridx (i);
4502 if (ri >= j) 4502 if (ri >= j)
4503 m_band(ri - j, j) = data(i); 4503 m_band(ri - j, j) = data (i);
4504 } 4504 }
4505 4505
4506 // Calculate the norm of the matrix, for later use. 4506 // Calculate the norm of the matrix, for later use.
4507 double anorm; 4507 double anorm;
4508 if (calc_cond) 4508 if (calc_cond)
4509 anorm = m_band.abs ().sum ().row(0).max (); 4509 anorm = m_band.abs ().sum ().row (0).max ();
4510 4510
4511 char job = 'L'; 4511 char job = 'L';
4512 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 4512 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4513 nr, n_lower, tmp_data, ldm, err 4513 nr, n_lower, tmp_data, ldm, err
4514 F77_CHAR_ARG_LEN (1))); 4514 F77_CHAR_ARG_LEN (1)));
4601 for (octave_idx_type i = 0; i < nc; i++) 4601 for (octave_idx_type i = 0; i < nc; i++)
4602 tmp_data[ii++] = 0.; 4602 tmp_data[ii++] = 0.;
4603 } 4603 }
4604 4604
4605 for (octave_idx_type j = 0; j < nc; j++) 4605 for (octave_idx_type j = 0; j < nc; j++)
4606 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4606 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4607 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); 4607 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4608 4608
4609 // Calculate the norm of the matrix, for later use. 4609 // Calculate the norm of the matrix, for later use.
4610 double anorm; 4610 double anorm;
4611 if (calc_cond) 4611 if (calc_cond)
4612 { 4612 {
4613 for (octave_idx_type j = 0; j < nr; j++) 4613 for (octave_idx_type j = 0; j < nr; j++)
4614 { 4614 {
4615 double atmp = 0.; 4615 double atmp = 0.;
4616 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4616 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4617 atmp += fabs(data(i)); 4617 atmp += fabs (data (i));
4618 if (atmp > anorm) 4618 if (atmp > anorm)
4619 anorm = atmp; 4619 anorm = atmp;
4620 } 4620 }
4621 } 4621 }
4622 4622
4744 for (octave_idx_type i = 0; i < nc; i++) 4744 for (octave_idx_type i = 0; i < nc; i++)
4745 tmp_data[ii++] = 0.; 4745 tmp_data[ii++] = 0.;
4746 } 4746 }
4747 4747
4748 for (octave_idx_type j = 0; j < nc; j++) 4748 for (octave_idx_type j = 0; j < nc; j++)
4749 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4749 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4750 { 4750 {
4751 octave_idx_type ri = ridx (i); 4751 octave_idx_type ri = ridx (i);
4752 if (ri >= j) 4752 if (ri >= j)
4753 m_band(ri - j, j) = data(i); 4753 m_band(ri - j, j) = data (i);
4754 } 4754 }
4755 4755
4756 // Calculate the norm of the matrix, for later use. 4756 // Calculate the norm of the matrix, for later use.
4757 double anorm; 4757 double anorm;
4758 if (calc_cond) 4758 if (calc_cond)
4759 anorm = m_band.abs ().sum ().row(0).max (); 4759 anorm = m_band.abs ().sum ().row (0).max ();
4760 4760
4761 char job = 'L'; 4761 char job = 'L';
4762 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 4762 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4763 nr, n_lower, tmp_data, ldm, err 4763 nr, n_lower, tmp_data, ldm, err
4764 F77_CHAR_ARG_LEN (1))); 4764 F77_CHAR_ARG_LEN (1)));
4818 // will be as many as in b 4818 // will be as many as in b
4819 volatile octave_idx_type x_nz = b.nnz (); 4819 volatile octave_idx_type x_nz = b.nnz ();
4820 volatile octave_idx_type ii = 0; 4820 volatile octave_idx_type ii = 0;
4821 retval = SparseMatrix (b_nr, b_nc, x_nz); 4821 retval = SparseMatrix (b_nr, b_nc, x_nz);
4822 4822
4823 retval.xcidx(0) = 0; 4823 retval.xcidx (0) = 0;
4824 for (volatile octave_idx_type j = 0; j < b_nc; j++) 4824 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4825 { 4825 {
4826 for (octave_idx_type i = 0; i < b_nr; i++) 4826 for (octave_idx_type i = 0; i < b_nr; i++)
4827 Bx[i] = b.elem (i, j); 4827 Bx[i] = b.elem (i, j);
4828 4828
4852 (b_nc - j) / b_nc; 4852 (b_nc - j) / b_nc;
4853 sz = (sz > 10 ? sz : 10) + x_nz; 4853 sz = (sz > 10 ? sz : 10) + x_nz;
4854 retval.change_capacity (sz); 4854 retval.change_capacity (sz);
4855 x_nz = sz; 4855 x_nz = sz;
4856 } 4856 }
4857 retval.xdata(ii) = tmp; 4857 retval.xdata (ii) = tmp;
4858 retval.xridx(ii++) = i; 4858 retval.xridx (ii++) = i;
4859 } 4859 }
4860 } 4860 }
4861 retval.xcidx(j+1) = ii; 4861 retval.xcidx (j+1) = ii;
4862 } 4862 }
4863 4863
4864 retval.maybe_compress (); 4864 retval.maybe_compress ();
4865 } 4865 }
4866 } 4866 }
4884 for (octave_idx_type i = 0; i < nc; i++) 4884 for (octave_idx_type i = 0; i < nc; i++)
4885 tmp_data[ii++] = 0.; 4885 tmp_data[ii++] = 0.;
4886 } 4886 }
4887 4887
4888 for (octave_idx_type j = 0; j < nc; j++) 4888 for (octave_idx_type j = 0; j < nc; j++)
4889 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4889 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4890 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); 4890 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4891 4891
4892 // Calculate the norm of the matrix, for later use. 4892 // Calculate the norm of the matrix, for later use.
4893 double anorm; 4893 double anorm;
4894 if (calc_cond) 4894 if (calc_cond)
4895 { 4895 {
4896 for (octave_idx_type j = 0; j < nr; j++) 4896 for (octave_idx_type j = 0; j < nr; j++)
4897 { 4897 {
4898 double atmp = 0.; 4898 double atmp = 0.;
4899 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 4899 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4900 atmp += fabs(data(i)); 4900 atmp += fabs (data (i));
4901 if (atmp > anorm) 4901 if (atmp > anorm)
4902 anorm = atmp; 4902 anorm = atmp;
4903 } 4903 }
4904 } 4904 }
4905 4905
4967 { 4967 {
4968 char job = 'N'; 4968 char job = 'N';
4969 volatile octave_idx_type x_nz = b.nnz (); 4969 volatile octave_idx_type x_nz = b.nnz ();
4970 octave_idx_type b_nc = b.cols (); 4970 octave_idx_type b_nc = b.cols ();
4971 retval = SparseMatrix (nr, b_nc, x_nz); 4971 retval = SparseMatrix (nr, b_nc, x_nz);
4972 retval.xcidx(0) = 0; 4972 retval.xcidx (0) = 0;
4973 volatile octave_idx_type ii = 0; 4973 volatile octave_idx_type ii = 0;
4974 4974
4975 OCTAVE_LOCAL_BUFFER (double, work, nr); 4975 OCTAVE_LOCAL_BUFFER (double, work, nr);
4976 4976
4977 for (volatile octave_idx_type j = 0; j < b_nc; j++) 4977 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4978 { 4978 {
4979 for (octave_idx_type i = 0; i < nr; i++) 4979 for (octave_idx_type i = 0; i < nr; i++)
4980 work[i] = 0.; 4980 work[i] = 0.;
4981 for (octave_idx_type i = b.cidx(j); 4981 for (octave_idx_type i = b.cidx (j);
4982 i < b.cidx(j+1); i++) 4982 i < b.cidx (j+1); i++)
4983 work[b.ridx(i)] = b.data(i); 4983 work[b.ridx (i)] = b.data (i);
4984 4984
4985 F77_XFCN (dgbtrs, DGBTRS, 4985 F77_XFCN (dgbtrs, DGBTRS,
4986 (F77_CONST_CHAR_ARG2 (&job, 1), 4986 (F77_CONST_CHAR_ARG2 (&job, 1),
4987 nr, n_lower, n_upper, 1, tmp_data, 4987 nr, n_lower, n_upper, 1, tmp_data,
4988 ldm, pipvt, work, b.rows (), err 4988 ldm, pipvt, work, b.rows (), err
5004 } 5004 }
5005 5005
5006 for (octave_idx_type i = 0; i < nr; i++) 5006 for (octave_idx_type i = 0; i < nr; i++)
5007 if (work[i] != 0.) 5007 if (work[i] != 0.)
5008 { 5008 {
5009 retval.xridx(ii) = i; 5009 retval.xridx (ii) = i;
5010 retval.xdata(ii++) = work[i]; 5010 retval.xdata (ii++) = work[i];
5011 } 5011 }
5012 retval.xcidx(j+1) = ii; 5012 retval.xcidx (j+1) = ii;
5013 } 5013 }
5014 5014
5015 retval.maybe_compress (); 5015 retval.maybe_compress ();
5016 } 5016 }
5017 } 5017 }
5062 for (octave_idx_type i = 0; i < nc; i++) 5062 for (octave_idx_type i = 0; i < nc; i++)
5063 tmp_data[ii++] = 0.; 5063 tmp_data[ii++] = 0.;
5064 } 5064 }
5065 5065
5066 for (octave_idx_type j = 0; j < nc; j++) 5066 for (octave_idx_type j = 0; j < nc; j++)
5067 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 5067 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5068 { 5068 {
5069 octave_idx_type ri = ridx (i); 5069 octave_idx_type ri = ridx (i);
5070 if (ri >= j) 5070 if (ri >= j)
5071 m_band(ri - j, j) = data(i); 5071 m_band(ri - j, j) = data (i);
5072 } 5072 }
5073 5073
5074 // Calculate the norm of the matrix, for later use. 5074 // Calculate the norm of the matrix, for later use.
5075 double anorm; 5075 double anorm;
5076 if (calc_cond) 5076 if (calc_cond)
5077 anorm = m_band.abs ().sum ().row(0).max (); 5077 anorm = m_band.abs ().sum ().row (0).max ();
5078 5078
5079 char job = 'L'; 5079 char job = 'L';
5080 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 5080 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5081 nr, n_lower, tmp_data, ldm, err 5081 nr, n_lower, tmp_data, ldm, err
5082 F77_CHAR_ARG_LEN (1))); 5082 F77_CHAR_ARG_LEN (1)));
5200 for (octave_idx_type i = 0; i < nc; i++) 5200 for (octave_idx_type i = 0; i < nc; i++)
5201 tmp_data[ii++] = 0.; 5201 tmp_data[ii++] = 0.;
5202 } 5202 }
5203 5203
5204 for (octave_idx_type j = 0; j < nc; j++) 5204 for (octave_idx_type j = 0; j < nc; j++)
5205 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 5205 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5206 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); 5206 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5207 5207
5208 // Calculate the norm of the matrix, for later use. 5208 // Calculate the norm of the matrix, for later use.
5209 double anorm; 5209 double anorm;
5210 if (calc_cond) 5210 if (calc_cond)
5211 { 5211 {
5212 for (octave_idx_type j = 0; j < nr; j++) 5212 for (octave_idx_type j = 0; j < nr; j++)
5213 { 5213 {
5214 double atmp = 0.; 5214 double atmp = 0.;
5215 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 5215 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5216 atmp += fabs(data(i)); 5216 atmp += fabs (data (i));
5217 if (atmp > anorm) 5217 if (atmp > anorm)
5218 anorm = atmp; 5218 anorm = atmp;
5219 } 5219 }
5220 } 5220 }
5221 5221
5361 for (octave_idx_type i = 0; i < nc; i++) 5361 for (octave_idx_type i = 0; i < nc; i++)
5362 tmp_data[ii++] = 0.; 5362 tmp_data[ii++] = 0.;
5363 } 5363 }
5364 5364
5365 for (octave_idx_type j = 0; j < nc; j++) 5365 for (octave_idx_type j = 0; j < nc; j++)
5366 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 5366 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5367 { 5367 {
5368 octave_idx_type ri = ridx (i); 5368 octave_idx_type ri = ridx (i);
5369 if (ri >= j) 5369 if (ri >= j)
5370 m_band(ri - j, j) = data(i); 5370 m_band(ri - j, j) = data (i);
5371 } 5371 }
5372 5372
5373 // Calculate the norm of the matrix, for later use. 5373 // Calculate the norm of the matrix, for later use.
5374 double anorm; 5374 double anorm;
5375 if (calc_cond) 5375 if (calc_cond)
5376 anorm = m_band.abs ().sum ().row(0).max (); 5376 anorm = m_band.abs ().sum ().row (0).max ();
5377 5377
5378 char job = 'L'; 5378 char job = 'L';
5379 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 5379 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5380 nr, n_lower, tmp_data, ldm, err 5380 nr, n_lower, tmp_data, ldm, err
5381 F77_CHAR_ARG_LEN (1))); 5381 F77_CHAR_ARG_LEN (1)));
5439 // will be as many as in b 5439 // will be as many as in b
5440 volatile octave_idx_type x_nz = b.nnz (); 5440 volatile octave_idx_type x_nz = b.nnz ();
5441 volatile octave_idx_type ii = 0; 5441 volatile octave_idx_type ii = 0;
5442 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); 5442 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
5443 5443
5444 retval.xcidx(0) = 0; 5444 retval.xcidx (0) = 0;
5445 for (volatile octave_idx_type j = 0; j < b_nc; j++) 5445 for (volatile octave_idx_type j = 0; j < b_nc; j++)
5446 { 5446 {
5447 5447
5448 for (octave_idx_type i = 0; i < b_nr; i++) 5448 for (octave_idx_type i = 0; i < b_nr; i++)
5449 { 5449 {
5497 } 5497 }
5498 5498
5499 for (octave_idx_type i = 0; i < nr; i++) 5499 for (octave_idx_type i = 0; i < nr; i++)
5500 if (Bx[i] != 0. || Bz[i] != 0.) 5500 if (Bx[i] != 0. || Bz[i] != 0.)
5501 { 5501 {
5502 retval.xridx(ii) = i; 5502 retval.xridx (ii) = i;
5503 retval.xdata(ii++) = 5503 retval.xdata (ii++) =
5504 Complex (Bx[i], Bz[i]); 5504 Complex (Bx[i], Bz[i]);
5505 } 5505 }
5506 5506
5507 retval.xcidx(j+1) = ii; 5507 retval.xcidx (j+1) = ii;
5508 } 5508 }
5509 5509
5510 retval.maybe_compress (); 5510 retval.maybe_compress ();
5511 } 5511 }
5512 } 5512 }
5530 for (octave_idx_type i = 0; i < nc; i++) 5530 for (octave_idx_type i = 0; i < nc; i++)
5531 tmp_data[ii++] = 0.; 5531 tmp_data[ii++] = 0.;
5532 } 5532 }
5533 5533
5534 for (octave_idx_type j = 0; j < nc; j++) 5534 for (octave_idx_type j = 0; j < nc; j++)
5535 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 5535 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5536 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); 5536 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5537 5537
5538 // Calculate the norm of the matrix, for later use. 5538 // Calculate the norm of the matrix, for later use.
5539 double anorm; 5539 double anorm;
5540 if (calc_cond) 5540 if (calc_cond)
5541 { 5541 {
5542 for (octave_idx_type j = 0; j < nr; j++) 5542 for (octave_idx_type j = 0; j < nr; j++)
5543 { 5543 {
5544 double atmp = 0.; 5544 double atmp = 0.;
5545 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 5545 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5546 atmp += fabs(data(i)); 5546 atmp += fabs (data (i));
5547 if (atmp > anorm) 5547 if (atmp > anorm)
5548 anorm = atmp; 5548 anorm = atmp;
5549 } 5549 }
5550 } 5550 }
5551 5551
5613 { 5613 {
5614 char job = 'N'; 5614 char job = 'N';
5615 volatile octave_idx_type x_nz = b.nnz (); 5615 volatile octave_idx_type x_nz = b.nnz ();
5616 octave_idx_type b_nc = b.cols (); 5616 octave_idx_type b_nc = b.cols ();
5617 retval = SparseComplexMatrix (nr, b_nc, x_nz); 5617 retval = SparseComplexMatrix (nr, b_nc, x_nz);
5618 retval.xcidx(0) = 0; 5618 retval.xcidx (0) = 0;
5619 volatile octave_idx_type ii = 0; 5619 volatile octave_idx_type ii = 0;
5620 5620
5621 OCTAVE_LOCAL_BUFFER (double, Bx, nr); 5621 OCTAVE_LOCAL_BUFFER (double, Bx, nr);
5622 OCTAVE_LOCAL_BUFFER (double, Bz, nr); 5622 OCTAVE_LOCAL_BUFFER (double, Bz, nr);
5623 5623
5626 for (octave_idx_type i = 0; i < nr; i++) 5626 for (octave_idx_type i = 0; i < nr; i++)
5627 { 5627 {
5628 Bx[i] = 0.; 5628 Bx[i] = 0.;
5629 Bz[i] = 0.; 5629 Bz[i] = 0.;
5630 } 5630 }
5631 for (octave_idx_type i = b.cidx(j); 5631 for (octave_idx_type i = b.cidx (j);
5632 i < b.cidx(j+1); i++) 5632 i < b.cidx (j+1); i++)
5633 { 5633 {
5634 Complex c = b.data(i); 5634 Complex c = b.data (i);
5635 Bx[b.ridx(i)] = std::real (c); 5635 Bx[b.ridx (i)] = std::real (c);
5636 Bz[b.ridx(i)] = std::imag (c); 5636 Bz[b.ridx (i)] = std::imag (c);
5637 } 5637 }
5638 5638
5639 F77_XFCN (dgbtrs, DGBTRS, 5639 F77_XFCN (dgbtrs, DGBTRS,
5640 (F77_CONST_CHAR_ARG2 (&job, 1), 5640 (F77_CONST_CHAR_ARG2 (&job, 1),
5641 nr, n_lower, n_upper, 1, tmp_data, 5641 nr, n_lower, n_upper, 1, tmp_data,
5664 } 5664 }
5665 5665
5666 for (octave_idx_type i = 0; i < nr; i++) 5666 for (octave_idx_type i = 0; i < nr; i++)
5667 if (Bx[i] != 0. || Bz[i] != 0.) 5667 if (Bx[i] != 0. || Bz[i] != 0.)
5668 { 5668 {
5669 retval.xridx(ii) = i; 5669 retval.xridx (ii) = i;
5670 retval.xdata(ii++) = 5670 retval.xdata (ii++) =
5671 Complex (Bx[i], Bz[i]); 5671 Complex (Bx[i], Bz[i]);
5672 } 5672 }
5673 retval.xcidx(j+1) = ii; 5673 retval.xcidx (j+1) = ii;
5674 } 5674 }
5675 5675
5676 retval.maybe_compress (); 5676 retval.maybe_compress ();
5677 } 5677 }
5678 } 5678 }
5933 retval.resize (b.rows (), b.cols ()); 5933 retval.resize (b.rows (), b.cols ());
5934 for (octave_idx_type j = 0; j < b.cols (); j++) 5934 for (octave_idx_type j = 0; j < b.cols (); j++)
5935 { 5935 {
5936 octave_idx_type jr = j * b.rows (); 5936 octave_idx_type jr = j * b.rows ();
5937 for (octave_idx_type i = 0; i < b.rows (); i++) 5937 for (octave_idx_type i = 0; i < b.rows (); i++)
5938 retval.xelem(i,j) = static_cast<double *>(X->x)[jr + i]; 5938 retval.xelem (i,j) = static_cast<double *>(X->x)[jr + i];
5939 } 5939 }
5940 5940
5941 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 5941 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5942 CHOLMOD_NAME(free_dense) (&X, cm); 5942 CHOLMOD_NAME(free_dense) (&X, cm);
5943 CHOLMOD_NAME(free_factor) (&L, cm); 5943 CHOLMOD_NAME(free_factor) (&L, cm);
6158 retval = SparseMatrix (static_cast<octave_idx_type>(X->nrow), 6158 retval = SparseMatrix (static_cast<octave_idx_type>(X->nrow),
6159 static_cast<octave_idx_type>(X->ncol), 6159 static_cast<octave_idx_type>(X->ncol),
6160 static_cast<octave_idx_type>(X->nzmax)); 6160 static_cast<octave_idx_type>(X->nzmax));
6161 for (octave_idx_type j = 0; 6161 for (octave_idx_type j = 0;
6162 j <= static_cast<octave_idx_type>(X->ncol); j++) 6162 j <= static_cast<octave_idx_type>(X->ncol); j++)
6163 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j]; 6163 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6164 for (octave_idx_type j = 0; 6164 for (octave_idx_type j = 0;
6165 j < static_cast<octave_idx_type>(X->nzmax); j++) 6165 j < static_cast<octave_idx_type>(X->nzmax); j++)
6166 { 6166 {
6167 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j]; 6167 retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
6168 retval.xdata(j) = static_cast<double *>(X->x)[j]; 6168 retval.xdata (j) = static_cast<double *>(X->x)[j];
6169 } 6169 }
6170 6170
6171 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6171 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6172 CHOLMOD_NAME(free_sparse) (&X, cm); 6172 CHOLMOD_NAME(free_sparse) (&X, cm);
6173 CHOLMOD_NAME(free_factor) (&L, cm); 6173 CHOLMOD_NAME(free_factor) (&L, cm);
6210 // will be as many as in b 6210 // will be as many as in b
6211 octave_idx_type x_nz = b.nnz (); 6211 octave_idx_type x_nz = b.nnz ();
6212 octave_idx_type ii = 0; 6212 octave_idx_type ii = 0;
6213 retval = SparseMatrix (b_nr, b_nc, x_nz); 6213 retval = SparseMatrix (b_nr, b_nc, x_nz);
6214 6214
6215 retval.xcidx(0) = 0; 6215 retval.xcidx (0) = 0;
6216 for (octave_idx_type j = 0; j < b_nc; j++) 6216 for (octave_idx_type j = 0; j < b_nc; j++)
6217 { 6217 {
6218 6218
6219 for (octave_idx_type i = 0; i < b_nr; i++) 6219 for (octave_idx_type i = 0; i < b_nr; i++)
6220 Bx[i] = b.elem (i, j); 6220 Bx[i] = b.elem (i, j);
6245 octave_idx_type sz = x_nz * (b_nc - j) / b_nc; 6245 octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
6246 sz = (sz > 10 ? sz : 10) + x_nz; 6246 sz = (sz > 10 ? sz : 10) + x_nz;
6247 retval.change_capacity (sz); 6247 retval.change_capacity (sz);
6248 x_nz = sz; 6248 x_nz = sz;
6249 } 6249 }
6250 retval.xdata(ii) = tmp; 6250 retval.xdata (ii) = tmp;
6251 retval.xridx(ii++) = i; 6251 retval.xridx (ii++) = i;
6252 } 6252 }
6253 } 6253 }
6254 retval.xcidx(j+1) = ii; 6254 retval.xcidx (j+1) = ii;
6255 } 6255 }
6256 6256
6257 retval.maybe_compress (); 6257 retval.maybe_compress ();
6258 6258
6259 UMFPACK_DNAME (report_info) (control, info); 6259 UMFPACK_DNAME (report_info) (control, info);
6410 retval.resize (b.rows (), b.cols ()); 6410 retval.resize (b.rows (), b.cols ());
6411 for (octave_idx_type j = 0; j < b.cols (); j++) 6411 for (octave_idx_type j = 0; j < b.cols (); j++)
6412 { 6412 {
6413 octave_idx_type jr = j * b.rows (); 6413 octave_idx_type jr = j * b.rows ();
6414 for (octave_idx_type i = 0; i < b.rows (); i++) 6414 for (octave_idx_type i = 0; i < b.rows (); i++)
6415 retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i]; 6415 retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i];
6416 } 6416 }
6417 6417
6418 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6418 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6419 CHOLMOD_NAME(free_dense) (&X, cm); 6419 CHOLMOD_NAME(free_dense) (&X, cm);
6420 CHOLMOD_NAME(free_factor) (&L, cm); 6420 CHOLMOD_NAME(free_factor) (&L, cm);
6655 (static_cast<octave_idx_type>(X->nrow), 6655 (static_cast<octave_idx_type>(X->nrow),
6656 static_cast<octave_idx_type>(X->ncol), 6656 static_cast<octave_idx_type>(X->ncol),
6657 static_cast<octave_idx_type>(X->nzmax)); 6657 static_cast<octave_idx_type>(X->nzmax));
6658 for (octave_idx_type j = 0; 6658 for (octave_idx_type j = 0;
6659 j <= static_cast<octave_idx_type>(X->ncol); j++) 6659 j <= static_cast<octave_idx_type>(X->ncol); j++)
6660 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j]; 6660 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6661 for (octave_idx_type j = 0; 6661 for (octave_idx_type j = 0;
6662 j < static_cast<octave_idx_type>(X->nzmax); j++) 6662 j < static_cast<octave_idx_type>(X->nzmax); j++)
6663 { 6663 {
6664 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j]; 6664 retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
6665 retval.xdata(j) = static_cast<Complex *>(X->x)[j]; 6665 retval.xdata (j) = static_cast<Complex *>(X->x)[j];
6666 } 6666 }
6667 6667
6668 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6668 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6669 CHOLMOD_NAME(free_sparse) (&X, cm); 6669 CHOLMOD_NAME(free_sparse) (&X, cm);
6670 CHOLMOD_NAME(free_factor) (&L, cm); 6670 CHOLMOD_NAME(free_factor) (&L, cm);
6710 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); 6710 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
6711 6711
6712 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr); 6712 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6713 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr); 6713 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6714 6714
6715 retval.xcidx(0) = 0; 6715 retval.xcidx (0) = 0;
6716 for (octave_idx_type j = 0; j < b_nc; j++) 6716 for (octave_idx_type j = 0; j < b_nc; j++)
6717 { 6717 {
6718 for (octave_idx_type i = 0; i < b_nr; i++) 6718 for (octave_idx_type i = 0; i < b_nr; i++)
6719 { 6719 {
6720 Complex c = b (i,j); 6720 Complex c = b (i,j);
6752 octave_idx_type sz = x_nz * (b_nc - j) / b_nc; 6752 octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
6753 sz = (sz > 10 ? sz : 10) + x_nz; 6753 sz = (sz > 10 ? sz : 10) + x_nz;
6754 retval.change_capacity (sz); 6754 retval.change_capacity (sz);
6755 x_nz = sz; 6755 x_nz = sz;
6756 } 6756 }
6757 retval.xdata(ii) = tmp; 6757 retval.xdata (ii) = tmp;
6758 retval.xridx(ii++) = i; 6758 retval.xridx (ii++) = i;
6759 } 6759 }
6760 } 6760 }
6761 retval.xcidx(j+1) = ii; 6761 retval.xcidx (j+1) = ii;
6762 } 6762 }
6763 6763
6764 retval.maybe_compress (); 6764 retval.maybe_compress ();
6765 6765
6766 UMFPACK_DNAME (report_info) (control, info); 6766 UMFPACK_DNAME (report_info) (control, info);
6968 { 6968 {
6969 (*current_liboctave_error_handler) ("unknown matrix type"); 6969 (*current_liboctave_error_handler) ("unknown matrix type");
6970 return ComplexMatrix (); 6970 return ComplexMatrix ();
6971 } 6971 }
6972 6972
6973 if (singular_fallback && mattype.type(false) == MatrixType::Rectangular) 6973 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6974 { 6974 {
6975 rcond = 1.; 6975 rcond = 1.;
6976 #ifdef USE_QRSOLVE 6976 #ifdef USE_QRSOLVE
6977 retval = qrsolve (*this, b, err); 6977 retval = qrsolve (*this, b, err);
6978 #else 6978 #else
7036 { 7036 {
7037 (*current_liboctave_error_handler) ("unknown matrix type"); 7037 (*current_liboctave_error_handler) ("unknown matrix type");
7038 return SparseComplexMatrix (); 7038 return SparseComplexMatrix ();
7039 } 7039 }
7040 7040
7041 if (singular_fallback && mattype.type(false) == MatrixType::Rectangular) 7041 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
7042 { 7042 {
7043 rcond = 1.; 7043 rcond = 1.;
7044 #ifdef USE_QRSOLVE 7044 #ifdef USE_QRSOLVE
7045 retval = qrsolve (*this, b, err); 7045 retval = qrsolve (*this, b, err);
7046 #else 7046 #else
7449 r.cidx (0) = 0; 7449 r.cidx (0) = 0;
7450 for (octave_idx_type i = 0; i < nc; i++) 7450 for (octave_idx_type i = 0; i < nc; i++)
7451 { 7451 {
7452 for (octave_idx_type j = 0; j < nr; j++) 7452 for (octave_idx_type j = 0; j < nr; j++)
7453 { 7453 {
7454 if (jj < cidx(i+1) && ridx(jj) == j) 7454 if (jj < cidx (i+1) && ridx (jj) == j)
7455 jj++; 7455 jj++;
7456 else 7456 else
7457 { 7457 {
7458 r.data(ii) = true; 7458 r.data (ii) = true;
7459 r.ridx(ii++) = j; 7459 r.ridx (ii++) = j;
7460 } 7460 }
7461 } 7461 }
7462 r.cidx (i+1) = ii; 7462 r.cidx (i+1) = ii;
7463 } 7463 }
7464 7464
7498 if ((rows () == 1 && dim == -1) || dim == 1) 7498 if ((rows () == 1 && dim == -1) || dim == 1)
7499 return transpose (). prod (0). transpose (); 7499 return transpose (). prod (0). transpose ();
7500 else 7500 else
7501 { 7501 {
7502 SPARSE_REDUCTION_OP (SparseMatrix, double, *=, 7502 SPARSE_REDUCTION_OP (SparseMatrix, double, *=,
7503 (cidx(j+1) - cidx(j) < nr ? 0.0 : 1.0), 1.0); 7503 (cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0);
7504 } 7504 }
7505 } 7505 }
7506 7506
7507 SparseMatrix 7507 SparseMatrix
7508 SparseMatrix::sum (int dim) const 7508 SparseMatrix::sum (int dim) const
7513 SparseMatrix 7513 SparseMatrix
7514 SparseMatrix::sumsq (int dim) const 7514 SparseMatrix::sumsq (int dim) const
7515 { 7515 {
7516 #define ROW_EXPR \ 7516 #define ROW_EXPR \
7517 double d = data (i); \ 7517 double d = data (i); \
7518 tmp[ridx(i)] += d * d 7518 tmp[ridx (i)] += d * d
7519 7519
7520 #define COL_EXPR \ 7520 #define COL_EXPR \
7521 double d = data (i); \ 7521 double d = data (i); \
7522 tmp[j] += d * d 7522 tmp[j] += d * d
7523 7523
7534 octave_idx_type nz = nnz (); 7534 octave_idx_type nz = nnz ();
7535 7535
7536 SparseMatrix retval (*this); 7536 SparseMatrix retval (*this);
7537 7537
7538 for (octave_idx_type i = 0; i < nz; i++) 7538 for (octave_idx_type i = 0; i < nz; i++)
7539 retval.data(i) = fabs(retval.data(i)); 7539 retval.data (i) = fabs (retval.data (i));
7540 7540
7541 return retval; 7541 return retval;
7542 } 7542 }
7543 7543
7544 SparseMatrix 7544 SparseMatrix
7561 // add one to the printed indices to go from 7561 // add one to the printed indices to go from
7562 // zero-based to one-based arrays 7562 // zero-based to one-based arrays
7563 for (octave_idx_type j = 0; j < nc; j++) 7563 for (octave_idx_type j = 0; j < nc; j++)
7564 { 7564 {
7565 octave_quit (); 7565 octave_quit ();
7566 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) 7566 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
7567 { 7567 {
7568 os << a.ridx(i) + 1 << " " << j + 1 << " "; 7568 os << a.ridx (i) + 1 << " " << j + 1 << " ";
7569 octave_write_double (os, a.data(i)); 7569 octave_write_double (os, a.data (i));
7570 os << "\n"; 7570 os << "\n";
7571 } 7571 }
7572 } 7572 }
7573 7573
7574 return os; 7574 return os;
7710 // Count the number of non-zero elements 7710 // Count the number of non-zero elements
7711 if (d < 0.) 7711 if (d < 0.)
7712 { 7712 {
7713 result = SparseMatrix (nr, nc, d); 7713 result = SparseMatrix (nr, nc, d);
7714 for (octave_idx_type j = 0; j < nc; j++) 7714 for (octave_idx_type j = 0; j < nc; j++)
7715 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) 7715 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7716 { 7716 {
7717 double tmp = xmin (d, m.data (i)); 7717 double tmp = xmin (d, m.data (i));
7718 if (tmp != 0.) 7718 if (tmp != 0.)
7719 { 7719 {
7720 octave_idx_type idx = m.ridx(i) + j * nr; 7720 octave_idx_type idx = m.ridx (i) + j * nr;
7721 result.xdata(idx) = tmp; 7721 result.xdata (idx) = tmp;
7722 result.xridx(idx) = m.ridx(i); 7722 result.xridx (idx) = m.ridx (i);
7723 } 7723 }
7724 } 7724 }
7725 } 7725 }
7726 else 7726 else
7727 { 7727 {
7728 octave_idx_type nel = 0; 7728 octave_idx_type nel = 0;
7729 for (octave_idx_type j = 0; j < nc; j++) 7729 for (octave_idx_type j = 0; j < nc; j++)
7730 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) 7730 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7731 if (xmin (d, m.data (i)) != 0.) 7731 if (xmin (d, m.data (i)) != 0.)
7732 nel++; 7732 nel++;
7733 7733
7734 result = SparseMatrix (nr, nc, nel); 7734 result = SparseMatrix (nr, nc, nel);
7735 7735
7736 octave_idx_type ii = 0; 7736 octave_idx_type ii = 0;
7737 result.xcidx(0) = 0; 7737 result.xcidx (0) = 0;
7738 for (octave_idx_type j = 0; j < nc; j++) 7738 for (octave_idx_type j = 0; j < nc; j++)
7739 { 7739 {
7740 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) 7740 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7741 { 7741 {
7742 double tmp = xmin (d, m.data (i)); 7742 double tmp = xmin (d, m.data (i));
7743 7743
7744 if (tmp != 0.) 7744 if (tmp != 0.)
7745 { 7745 {
7746 result.xdata(ii) = tmp; 7746 result.xdata (ii) = tmp;
7747 result.xridx(ii++) = m.ridx(i); 7747 result.xridx (ii++) = m.ridx (i);
7748 } 7748 }
7749 } 7749 }
7750 result.xcidx(j+1) = ii; 7750 result.xcidx (j+1) = ii;
7751 } 7751 }
7752 } 7752 }
7753 7753
7754 return result; 7754 return result;
7755 } 7755 }
7781 7781
7782 octave_idx_type jx = 0; 7782 octave_idx_type jx = 0;
7783 r.cidx (0) = 0; 7783 r.cidx (0) = 0;
7784 for (octave_idx_type i = 0 ; i < a_nc ; i++) 7784 for (octave_idx_type i = 0 ; i < a_nc ; i++)
7785 { 7785 {
7786 octave_idx_type ja = a.cidx(i); 7786 octave_idx_type ja = a.cidx (i);
7787 octave_idx_type ja_max = a.cidx(i+1); 7787 octave_idx_type ja_max = a.cidx (i+1);
7788 bool ja_lt_max= ja < ja_max; 7788 bool ja_lt_max= ja < ja_max;
7789 7789
7790 octave_idx_type jb = b.cidx(i); 7790 octave_idx_type jb = b.cidx (i);
7791 octave_idx_type jb_max = b.cidx(i+1); 7791 octave_idx_type jb_max = b.cidx (i+1);
7792 bool jb_lt_max = jb < jb_max; 7792 bool jb_lt_max = jb < jb_max;
7793 7793
7794 while (ja_lt_max || jb_lt_max ) 7794 while (ja_lt_max || jb_lt_max )
7795 { 7795 {
7796 octave_quit (); 7796 octave_quit ();
7797 if ((! jb_lt_max) || 7797 if ((! jb_lt_max) ||
7798 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) 7798 (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7799 { 7799 {
7800 double tmp = xmin (a.data(ja), 0.); 7800 double tmp = xmin (a.data (ja), 0.);
7801 if (tmp != 0.) 7801 if (tmp != 0.)
7802 { 7802 {
7803 r.ridx(jx) = a.ridx(ja); 7803 r.ridx (jx) = a.ridx (ja);
7804 r.data(jx) = tmp; 7804 r.data (jx) = tmp;
7805 jx++; 7805 jx++;
7806 } 7806 }
7807 ja++; 7807 ja++;
7808 ja_lt_max= ja < ja_max; 7808 ja_lt_max= ja < ja_max;
7809 } 7809 }
7810 else if (( !ja_lt_max ) || 7810 else if (( !ja_lt_max ) ||
7811 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) 7811 (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
7812 { 7812 {
7813 double tmp = xmin (0., b.data(jb)); 7813 double tmp = xmin (0., b.data (jb));
7814 if (tmp != 0.) 7814 if (tmp != 0.)
7815 { 7815 {
7816 r.ridx(jx) = b.ridx(jb); 7816 r.ridx (jx) = b.ridx (jb);
7817 r.data(jx) = tmp; 7817 r.data (jx) = tmp;
7818 jx++; 7818 jx++;
7819 } 7819 }
7820 jb++; 7820 jb++;
7821 jb_lt_max= jb < jb_max; 7821 jb_lt_max= jb < jb_max;
7822 } 7822 }
7823 else 7823 else
7824 { 7824 {
7825 double tmp = xmin (a.data(ja), b.data(jb)); 7825 double tmp = xmin (a.data (ja), b.data (jb));
7826 if (tmp != 0.) 7826 if (tmp != 0.)
7827 { 7827 {
7828 r.data(jx) = tmp; 7828 r.data (jx) = tmp;
7829 r.ridx(jx) = a.ridx(ja); 7829 r.ridx (jx) = a.ridx (ja);
7830 jx++; 7830 jx++;
7831 } 7831 }
7832 ja++; 7832 ja++;
7833 ja_lt_max= ja < ja_max; 7833 ja_lt_max= ja < ja_max;
7834 jb++; 7834 jb++;
7835 jb_lt_max= jb < jb_max; 7835 jb_lt_max= jb < jb_max;
7836 } 7836 }
7837 } 7837 }
7838 r.cidx(i+1) = jx; 7838 r.cidx (i+1) = jx;
7839 } 7839 }
7840 7840
7841 r.maybe_compress (); 7841 r.maybe_compress ();
7842 } 7842 }
7843 } 7843 }
7860 // Count the number of non-zero elements 7860 // Count the number of non-zero elements
7861 if (d > 0.) 7861 if (d > 0.)
7862 { 7862 {
7863 result = SparseMatrix (nr, nc, d); 7863 result = SparseMatrix (nr, nc, d);
7864 for (octave_idx_type j = 0; j < nc; j++) 7864 for (octave_idx_type j = 0; j < nc; j++)
7865 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) 7865 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7866 { 7866 {
7867 double tmp = xmax (d, m.data (i)); 7867 double tmp = xmax (d, m.data (i));
7868 7868
7869 if (tmp != 0.) 7869 if (tmp != 0.)
7870 { 7870 {
7871 octave_idx_type idx = m.ridx(i) + j * nr; 7871 octave_idx_type idx = m.ridx (i) + j * nr;
7872 result.xdata(idx) = tmp; 7872 result.xdata (idx) = tmp;
7873 result.xridx(idx) = m.ridx(i); 7873 result.xridx (idx) = m.ridx (i);
7874 } 7874 }
7875 } 7875 }
7876 } 7876 }
7877 else 7877 else
7878 { 7878 {
7879 octave_idx_type nel = 0; 7879 octave_idx_type nel = 0;
7880 for (octave_idx_type j = 0; j < nc; j++) 7880 for (octave_idx_type j = 0; j < nc; j++)
7881 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) 7881 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7882 if (xmax (d, m.data (i)) != 0.) 7882 if (xmax (d, m.data (i)) != 0.)
7883 nel++; 7883 nel++;
7884 7884
7885 result = SparseMatrix (nr, nc, nel); 7885 result = SparseMatrix (nr, nc, nel);
7886 7886
7887 octave_idx_type ii = 0; 7887 octave_idx_type ii = 0;
7888 result.xcidx(0) = 0; 7888 result.xcidx (0) = 0;
7889 for (octave_idx_type j = 0; j < nc; j++) 7889 for (octave_idx_type j = 0; j < nc; j++)
7890 { 7890 {
7891 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) 7891 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7892 { 7892 {
7893 double tmp = xmax (d, m.data (i)); 7893 double tmp = xmax (d, m.data (i));
7894 if (tmp != 0.) 7894 if (tmp != 0.)
7895 { 7895 {
7896 result.xdata(ii) = tmp; 7896 result.xdata (ii) = tmp;
7897 result.xridx(ii++) = m.ridx(i); 7897 result.xridx (ii++) = m.ridx (i);
7898 } 7898 }
7899 } 7899 }
7900 result.xcidx(j+1) = ii; 7900 result.xcidx (j+1) = ii;
7901 } 7901 }
7902 } 7902 }
7903 7903
7904 return result; 7904 return result;
7905 } 7905 }
7931 7931
7932 octave_idx_type jx = 0; 7932 octave_idx_type jx = 0;
7933 r.cidx (0) = 0; 7933 r.cidx (0) = 0;
7934 for (octave_idx_type i = 0 ; i < a_nc ; i++) 7934 for (octave_idx_type i = 0 ; i < a_nc ; i++)
7935 { 7935 {
7936 octave_idx_type ja = a.cidx(i); 7936 octave_idx_type ja = a.cidx (i);
7937 octave_idx_type ja_max = a.cidx(i+1); 7937 octave_idx_type ja_max = a.cidx (i+1);
7938 bool ja_lt_max= ja < ja_max; 7938 bool ja_lt_max= ja < ja_max;
7939 7939
7940 octave_idx_type jb = b.cidx(i); 7940 octave_idx_type jb = b.cidx (i);
7941 octave_idx_type jb_max = b.cidx(i+1); 7941 octave_idx_type jb_max = b.cidx (i+1);
7942 bool jb_lt_max = jb < jb_max; 7942 bool jb_lt_max = jb < jb_max;
7943 7943
7944 while (ja_lt_max || jb_lt_max ) 7944 while (ja_lt_max || jb_lt_max )
7945 { 7945 {
7946 octave_quit (); 7946 octave_quit ();
7947 if ((! jb_lt_max) || 7947 if ((! jb_lt_max) ||
7948 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) 7948 (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7949 { 7949 {
7950 double tmp = xmax (a.data(ja), 0.); 7950 double tmp = xmax (a.data (ja), 0.);
7951 if (tmp != 0.) 7951 if (tmp != 0.)
7952 { 7952 {
7953 r.ridx(jx) = a.ridx(ja); 7953 r.ridx (jx) = a.ridx (ja);
7954 r.data(jx) = tmp; 7954 r.data (jx) = tmp;
7955 jx++; 7955 jx++;
7956 } 7956 }
7957 ja++; 7957 ja++;
7958 ja_lt_max= ja < ja_max; 7958 ja_lt_max= ja < ja_max;
7959 } 7959 }
7960 else if (( !ja_lt_max ) || 7960 else if (( !ja_lt_max ) ||
7961 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) 7961 (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
7962 { 7962 {
7963 double tmp = xmax (0., b.data(jb)); 7963 double tmp = xmax (0., b.data (jb));
7964 if (tmp != 0.) 7964 if (tmp != 0.)
7965 { 7965 {
7966 r.ridx(jx) = b.ridx(jb); 7966 r.ridx (jx) = b.ridx (jb);
7967 r.data(jx) = tmp; 7967 r.data (jx) = tmp;
7968 jx++; 7968 jx++;
7969 } 7969 }
7970 jb++; 7970 jb++;
7971 jb_lt_max= jb < jb_max; 7971 jb_lt_max= jb < jb_max;
7972 } 7972 }
7973 else 7973 else
7974 { 7974 {
7975 double tmp = xmax (a.data(ja), b.data(jb)); 7975 double tmp = xmax (a.data (ja), b.data (jb));
7976 if (tmp != 0.) 7976 if (tmp != 0.)
7977 { 7977 {
7978 r.data(jx) = tmp; 7978 r.data (jx) = tmp;
7979 r.ridx(jx) = a.ridx(ja); 7979 r.ridx (jx) = a.ridx (ja);
7980 jx++; 7980 jx++;
7981 } 7981 }
7982 ja++; 7982 ja++;
7983 ja_lt_max= ja < ja_max; 7983 ja_lt_max= ja < ja_max;
7984 jb++; 7984 jb++;
7985 jb_lt_max= jb < jb_max; 7985 jb_lt_max= jb < jb_max;
7986 } 7986 }
7987 } 7987 }
7988 r.cidx(i+1) = jx; 7988 r.cidx (i+1) = jx;
7989 } 7989 }
7990 7990
7991 r.maybe_compress (); 7991 r.maybe_compress ();
7992 } 7992 }
7993 } 7993 }