octave-nkf: liboctave/sparse-dmsolve.cc annotate

annotate liboctave/sparse-dmsolve.cc @ 16088:b29b10fbb744 stable release-3-6-4

Version 3.6.4 released. * configure.ac (AC_INIT): Version is now 3.6.4. (OCTAVE_RELEASE_DATE): Now 2013-02-21.

author	John W. Eaton <jwe@octave.org>
date	Thu, 21 Feb 2013 15:17:54 -0500
parents	1316bfc6e260
children	4328e28414aa

rev	line source
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	1 /*
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	2
14138 72c96de7a403 maint: update copyright notices for 2012 John W. Eaton <jwe@octave.org> parents: 11586 diff changeset	3 Copyright (C) 2006-2012 David Bateman
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	4
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6924 diff changeset	5 This file is part of Octave.
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6924 diff changeset	6
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	8 under the terms of the GNU General Public License as published by the
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6924 diff changeset	9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6924 diff changeset	10 option) any later version.
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	11
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	15 for more details.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	16
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	17 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6924 diff changeset	18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6924 diff changeset	19 <http://www.gnu.org/licenses/>.
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	20
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	21 */
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	22
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	23 #ifdef HAVE_CONFIG_H
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	24 #include <config.h>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	25 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	26
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5794 diff changeset	27 #include <vector>
5794 1138ced03f14 [project @ 2006-05-08 20:23:04 by jwe] jwe parents: 5792 diff changeset	28
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	29 #include "MArray.h"
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	30 #include "MSparse.h"
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	31 #include "SparseQR.h"
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	32 #include "SparseCmplxQR.h"
5785 6b9cec830d72 [project @ 2006-05-03 19:32:46 by dbateman] dbateman parents: 5775 diff changeset	33 #include "MatrixType.h"
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	34 #include "oct-sort.h"
8377 25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER Jaroslav Hajek <highegg@gmail.com> parents: 7520 diff changeset	35 #include "oct-locbuf.h"
14888 4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	36 #include "oct-inttypes.h"
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	37
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	38 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	39 static MSparse<T>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	40 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	41 const octave_idx_type *Q, octave_idx_type rst,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	42 octave_idx_type rend, octave_idx_type cst,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	43 octave_idx_type cend, octave_idx_type maxnz = -1,
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	44 bool lazy = false)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	45 {
14888 4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	46 octave_idx_type nr = rend - rst, nc = cend - cst;
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	47 maxnz = (maxnz < 0 ? A.nnz () : maxnz);
14888 4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	48 octave_idx_type nz;
4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	49
4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	50 // Cast to uint64 to handle overflow in this multiplication
4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	51 if (octave_uint64 (nr)*octave_uint64 (nc) < octave_uint64 (maxnz))
14999 1316bfc6e260 Fix think-o in 4315a39da4c9 Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14888 diff changeset	52 nz = nr*nc;
14888 4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	53 else
4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	54 nz = maxnz;
4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	55
4315a39da4c9 Do computations with octave_uint64 to avoid overflow Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14138 diff changeset	56 MSparse<T> B (nr, nc, (nz < maxnz ? nz : maxnz));
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	57 // Some sparse functions can support lazy indexing (where elements
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	58 // in the row are in no particular order), even though octave in
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	59 // general can't. For those functions that can using it is a big
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	60 // win here in terms of speed.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	61 if (lazy)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	62 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	63 nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	64 for (octave_idx_type j = cst ; j < cend ; j++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	65 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	66 octave_idx_type qq = (Q ? Q [j] : j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	67 B.xcidx (j - cst) = nz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	68 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	69 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	70 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	71 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	72 if (r >= rst && r < rend)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	73 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	74 B.xdata (nz) = A.data (p);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	75 B.xridx (nz++) = r - rst ;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	76 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	77 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	78 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	79 B.xcidx (cend - cst) = nz ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	80 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	81 else
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	82 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	83 OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	84 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	85 octave_idx_type *ri = B.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	86 nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	87 for (octave_idx_type j = cst ; j < cend ; j++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	88 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	89 octave_idx_type qq = (Q ? Q [j] : j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	90 B.xcidx (j - cst) = nz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	91 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	92 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	93 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	94 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	95 if (r >= rst && r < rend)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	96 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	97 X [r-rst] = A.data (p);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	98 B.xridx (nz++) = r - rst ;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	99 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	100 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	101 sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	102 for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	103 B.xdata (p) = X [B.xridx (p)];
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	104 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	105 B.xcidx (cend - cst) = nz ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	106 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	107
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	108 return B;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	109 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	110
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	111 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	112 static MSparse<double>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	113 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	114 const octave_idx_type *Q, octave_idx_type rst,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	115 octave_idx_type rend, octave_idx_type cst,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	116 octave_idx_type cend, octave_idx_type maxnz,
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	117 bool lazy);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	118
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	119 static MSparse<Complex>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	120 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	121 const octave_idx_type *Q, octave_idx_type rst,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	122 octave_idx_type rend, octave_idx_type cst,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	123 octave_idx_type cend, octave_idx_type maxnz,
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	124 bool lazy);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	125 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	126
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	127 template <class T>
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	128 static MArray<T>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	129 dmsolve_extract (const MArray<T> &m, const octave_idx_type *,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	130 const octave_idx_type *, octave_idx_type r1,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	131 octave_idx_type r2, octave_idx_type c1,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	132 octave_idx_type c2)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	133 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	134 r2 -= 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	135 c2 -= 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	136 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	137 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	138
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	139 octave_idx_type new_r = r2 - r1 + 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	140 octave_idx_type new_c = c2 - c1 + 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	141
11570 57632dea2446 attempt better backward compatibility for Array constructors John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	142 MArray<T> result (dim_vector (new_r, new_c));
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	143
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	144 for (octave_idx_type j = 0; j < new_c; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	145 for (octave_idx_type i = 0; i < new_r; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	146 result.xelem (i, j) = m.elem (r1+i, c1+j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	147
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	148 return result;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	149 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	150
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	151 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	152 static MArray<double>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	153 dmsolve_extract (const MArray<double> &m, const octave_idx_type *,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	154 const octave_idx_type *, octave_idx_type r1,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	155 octave_idx_type r2, octave_idx_type c1,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	156 octave_idx_type c2)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	157
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	158 static MArray<Complex>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	159 dmsolve_extract (const MArray<Complex> &m, const octave_idx_type *,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	160 const octave_idx_type *, octave_idx_type r1,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	161 octave_idx_type r2, octave_idx_type c1,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	162 octave_idx_type c2)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	163 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	164
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	165 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	166 static void
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	167 dmsolve_insert (MArray<T> &a, const MArray<T> &b, const octave_idx_type *Q,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	168 octave_idx_type r, octave_idx_type c)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	169 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	170 T *ax = a.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	171 const T *bx = b.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	172 octave_idx_type anr = a.rows();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	173 octave_idx_type nr = b.rows();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	174 octave_idx_type nc = b.cols();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	175 for (octave_idx_type j = 0; j < nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	176 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	177 octave_idx_type aoff = (c + j) * anr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	178 octave_idx_type boff = j * nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	179 for (octave_idx_type i = 0; i < nr; i++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	180 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	181 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	182 ax [Q [r + i] + aoff] = bx [i + boff];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	183 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	184 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	185 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	186
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	187 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	188 static void
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	189 dmsolve_insert (MArray<double> &a, const MArray<double> &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	190 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	191
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	192 static void
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	193 dmsolve_insert (MArray<Complex> &a, const MArray<Complex> &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	194 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	195 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	196
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	197 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	198 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	199 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	200 octave_idx_type r, octave_idx_type c)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	201 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	202 octave_idx_type b_rows = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	203 octave_idx_type b_cols = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	204 octave_idx_type nr = a.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	205 octave_idx_type nc = a.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	206
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	207 OCTAVE_LOCAL_BUFFER (octave_idx_type, Qinv, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	208 for (octave_idx_type i = 0; i < nr; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	209 Qinv [Q [i]] = i;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	210
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	211 // First count the number of elements in the final array
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	212 octave_idx_type nel = a.xcidx(c) + b.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	213
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	214 if (c + b_cols < nc)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	215 nel += a.xcidx(nc) - a.xcidx(c + b_cols);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	216
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	217 for (octave_idx_type i = c; i < c + b_cols; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	218 for (octave_idx_type j = a.xcidx(i); j < a.xcidx(i+1); j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	219 if (Qinv [a.xridx(j)] < r \|\| Qinv [a.xridx(j)] >= r + b_rows)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	220 nel++;
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	221
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	222 OCTAVE_LOCAL_BUFFER (T, X, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	223 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	224 MSparse<T> tmp (a);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	225 a = MSparse<T> (nr, nc, nel);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	226 octave_idx_type *ri = a.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	227
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	228 for (octave_idx_type i = 0; i < tmp.cidx(c); i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	229 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	230 a.xdata(i) = tmp.xdata(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	231 a.xridx(i) = tmp.xridx(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	232 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	233 for (octave_idx_type i = 0; i < c + 1; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	234 a.xcidx(i) = tmp.xcidx(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	235
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	236 octave_idx_type ii = a.xcidx(c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	237
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	238 for (octave_idx_type i = c; i < c + b_cols; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	239 {
10142 829e69ec3110 make OCTAVE_QUIT a function Jaroslav Hajek <highegg@gmail.com> parents: 8920 diff changeset	240 octave_quit ();
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	241
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	242 for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	243 if (Qinv [tmp.xridx(j)] < r \|\| Qinv [tmp.xridx(j)] >= r + b_rows)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	244 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	245 X [tmp.xridx(j)] = tmp.xdata(j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	246 a.xridx(ii++) = tmp.xridx(j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	247 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	248
10142 829e69ec3110 make OCTAVE_QUIT a function Jaroslav Hajek <highegg@gmail.com> parents: 8920 diff changeset	249 octave_quit ();
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	250
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	251 for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	252 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	253 X [Q [r + b.ridx(j)]] = b.data(j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	254 a.xridx(ii++) = Q [r + b.ridx(j)];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	255 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	256
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	257 sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	258 for (octave_idx_type p = a.xcidx (i); p < ii; p++)
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	259 a.xdata (p) = X [a.xridx (p)];
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	260 a.xcidx(i+1) = ii;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	261 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	262
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	263 for (octave_idx_type i = c + b_cols; i < nc; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	264 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	265 for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	266 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	267 a.xdata(ii) = tmp.xdata(j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	268 a.xridx(ii++) = tmp.xridx(j);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	269 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	270 a.xcidx(i+1) = ii;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	271 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	272 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	273
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	274 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	275 static void
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	276 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	277 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	278
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	279 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	280 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	281 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	282 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	283
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	284 template <class T, class RT>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	285 static void
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	286 dmsolve_permute (MArray<RT> &a, const MArray<T>& b, const octave_idx_type *p)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	287 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	288 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	289 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	290 const T *Bx = b.fortran_vec();
11574 a83bad07f7e3 attempt better backward compatibility for Array resize functions John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	291 a.resize (dim_vector (b_nr, b_nc));
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	292 RT *Btx = a.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	293 for (octave_idx_type j = 0; j < b_nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	294 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	295 octave_idx_type off = j * b_nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	296 for (octave_idx_type i = 0; i < b_nr; i++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	297 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	298 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	299 Btx [p [i] + off] = Bx [ i + off];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	300 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	301 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	302 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	303
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	304 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	305 static void
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	306 dmsolve_permute (MArray<double> &a, const MArray<double>& b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	307 const octave_idx_type *p);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	308
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	309 static void
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	310 dmsolve_permute (MArray<Complex> &a, const MArray<double>& b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	311 const octave_idx_type *p);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	312
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	313 static void
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	314 dmsolve_permute (MArray<Complex> &a, const MArray<Complex>& b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	315 const octave_idx_type *p);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	316 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	317
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	318 template <class T, class RT>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	319 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	320 dmsolve_permute (MSparse<RT> &a, const MSparse<T>& b, const octave_idx_type *p)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	321 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	322 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	323 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	324 octave_idx_type b_nz = b.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	325 octave_idx_type nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	326 a = MSparse<RT> (b_nr, b_nc, b_nz);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	327 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	328 octave_idx_type *ri = a.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	329 OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	330 a.xcidx(0) = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	331 for (octave_idx_type j = 0; j < b_nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	332 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	333 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	334 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	335 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	336 octave_idx_type r = p [b.ridx (i)];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	337 X [r] = b.data (i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	338 a.xridx(nz++) = p [b.ridx (i)];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	339 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	340 sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	341 for (octave_idx_type i = a.cidx (j); i < nz; i++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	342 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	343 octave_quit ();
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	344 a.xdata (i) = X [a.xridx (i)];
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	345 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	346 a.xcidx(j+1) = nz;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	347 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	348 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	349
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	350 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	351 static void
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	352 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	353 const octave_idx_type *p);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	354
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	355 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	356 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	357 const octave_idx_type *p);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	358
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	359 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	360 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	361 const octave_idx_type *p);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	362 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	363
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	364 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	365 solve_singularity_warning (double)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	366 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	367 // Dummy singularity handler so that LU solver doesn't flag
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	368 // an error for numerically rank defficient matrices
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	369 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	370
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	371 template <class RT, class ST, class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	372 RT
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	373 dmsolve (const ST &a, const T &b, octave_idx_type &info)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	374 {
5684 f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman] dbateman parents: 5683 diff changeset	375 #ifdef HAVE_CXSPARSE
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	376 octave_idx_type nr = a.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	377 octave_idx_type nc = a.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	378 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	379 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	380 RT retval;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	381
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	382 if (nr < 0 \|\| nc < 0 \|\| nr != b_nr)
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	383 (*current_liboctave_error_handler)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	384 ("matrix dimension mismatch in solution of minimum norm problem");
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	385 else if (nr == 0 \|\| nc == 0 \|\| b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	386 retval = RT (nc, b_nc, 0.0);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	387 else
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	388 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	389 octave_idx_type nnz_remaining = a.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	390 CXSPARSE_DNAME () csm;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	391 csm.m = nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	392 csm.n = nc;
7520 b166043585a8 use 0 instead of NULL John W. Eaton <jwe@octave.org> parents: 7017 diff changeset	393 csm.x = 0;
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	394 csm.nz = -1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	395 csm.nzmax = a.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	396 // Cast away const on A, with full knowledge that CSparse won't touch it.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	397 // Prevents the methods below making a copy of the data.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	398 csm.p = const_cast<octave_idx_type *>(a.cidx ());
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	399 csm.i = const_cast<octave_idx_type *>(a.ridx ());
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	400
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5785 diff changeset	401 #if defined(CS_VER) && (CS_VER >= 2)
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5785 diff changeset	402 CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm, 0);
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5785 diff changeset	403 octave_idx_type *p = dm->p;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5785 diff changeset	404 octave_idx_type *q = dm->q;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5785 diff changeset	405 #else
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	406 CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	407 octave_idx_type *p = dm->P;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	408 octave_idx_type *q = dm->Q;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5785 diff changeset	409 #endif
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	410 OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	411 for (octave_idx_type i = 0; i < nr; i++)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	412 pinv [p [i]] = i;
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	413 RT btmp;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	414 dmsolve_permute (btmp, b, pinv);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	415 info = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	416 retval.resize (nc, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	417
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	418 // Leading over-determined block
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	419 if (dm->rr [2] < nr && dm->cc [3] < nc)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	420 {
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	421 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	422 nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	423 nnz_remaining -= m.nnz();
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	424 RT mtmp =
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	425 qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	426 b_nc), info);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	427 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	428 if (dm->rr [2] > 0 && !info)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	429 {
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	430 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2],
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	431 dm->cc [3], nc, nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	432 nnz_remaining -= m.nnz();
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	433 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	434 dm->rr[2], 0, b_nc);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	435 btmp.insert (ctmp - m * mtmp, 0, 0);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	436 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	437 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	438
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	439 // Structurally non-singular blocks
5775 ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe] jwe parents: 5697 diff changeset	440 // FIXME Should use fine Dulmange-Mendelsohn decomposition here.
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5794 diff changeset	441 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	442 {
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	443 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2],
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	444 dm->cc [2], dm->cc [3], nnz_remaining, false);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	445 nnz_remaining -= m.nnz();
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	446 RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2],
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	447 0, b_nc);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	448 double rcond = 0.0;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	449 MatrixType mtyp (MatrixType::Full);
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	450 RT mtmp = m.solve (mtyp, btmp2, info, rcond,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	451 solve_singularity_warning, false);
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	452 if (info != 0)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	453 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	454 info = 0;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	455 mtmp = qrsolve (m, btmp2, info);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	456 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	457
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	458 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	459 if (dm->rr [1] > 0 && !info)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	460 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	461 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	462 dm->cc [3], nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	463 nnz_remaining -= m.nnz();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	464 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	465 dm->rr[1], 0, b_nc);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	466 btmp.insert (ctmp - m * mtmp, 0, 0);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	467 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	468 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	469
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	470 // Trailing under-determined block
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5794 diff changeset	471 if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	472 {
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	473 ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	474 dm->cc [2], nnz_remaining, true);
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	475 RT mtmp =
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	476 qrsolve (m, dmsolve_extract(btmp, 0, 0, 0, dm->rr [1] , 0,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	477 b_nc), info);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	478 dmsolve_insert (retval, mtmp, q, 0, 0);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	479 }
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	480
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	481 CXSPARSE_DNAME (_dfree) (dm);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	482 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	483 return retval;
5684 f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman] dbateman parents: 5683 diff changeset	484 #else
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman] dbateman parents: 5683 diff changeset	485 return RT ();
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman] dbateman parents: 5683 diff changeset	486 #endif
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	487 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	488
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	489 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	490 extern Matrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	491 dmsolve (const SparseMatrix &a, const Matrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	492 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	493
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	494 extern ComplexMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	495 dmsolve (const SparseMatrix &a, const ComplexMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	496 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	497
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	498 extern ComplexMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	499 dmsolve (const SparseComplexMatrix &a, const Matrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	500 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	501
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	502 extern ComplexMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	503 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	504 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	505
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	506 extern SparseMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	507 dmsolve (const SparseMatrix &a, const SparseMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	508 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	509
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	510 extern SparseComplexMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	511 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	512 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	513
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	514 extern SparseComplexMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	515 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	516 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	517
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	518 extern SparseComplexMatrix
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11574 diff changeset	519 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	520 octave_idx_type &info);
5683 a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman] dbateman parents: diff changeset	521 #endif

Mercurial > octave-nkf

annotate liboctave/sparse-dmsolve.cc @ 16088:b29b10fbb744 stable release-3-6-4