octave: liboctave/SparseQR.cc annotate

annotate liboctave/SparseQR.cc @ 11523:fd0a3ac60b0e

update copyright notices

author	John W. Eaton <jwe@octave.org>
date	Fri, 14 Jan 2011 05:47:45 -0500
parents	f3d0b69c4628
children	feb90c7cfc16

rev	line source
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	1 /*
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	2
11523 fd0a3ac60b0e update copyright notices John W. Eaton <jwe@octave.org> parents: 11510 diff changeset	3 Copyright (C) 2005-2011 David Bateman
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 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
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
9761b7d24e9e [project @ 2006-02-09 09:12:02 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.
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	11
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	15 for more details.
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	16
9761b7d24e9e [project @ 2006-02-09 09:12:02 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/>.
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	20
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	21 */
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	22
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	23 #ifdef HAVE_CONFIG_H
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	24 #include <config.h>
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	25 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	26 #include <vector>
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	27
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	28 #include "lo-error.h"
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	29 #include "SparseQR.h"
8377 25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER Jaroslav Hajek <highegg@gmail.com> parents: 7789 diff changeset	30 #include "oct-locbuf.h"
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	31
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	32 SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order)
11510 f3d0b69c4628 SparseQR.cc, SparseCmplxQR.cc: fix member initialization order John W. Eaton <jwe@octave.org> parents: 11498 diff changeset	33 : count (1), nrows (0), S (0), N (0)
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	34 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	35 #ifdef HAVE_CXSPARSE
5648 69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	36 CXSPARSE_DNAME () A;
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	37 A.nzmax = a.nnz ();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	38 A.m = a.rows ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	39 A.n = a.cols ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	40 nrows = A.m;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	41 // Cast away const on A, with full knowledge that CSparse won't touch it
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	42 // Prevents the methods below making a copy of the data.
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	43 A.p = const_cast<octave_idx_type *>(a.cidx ());
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	44 A.i = const_cast<octave_idx_type *>(a.ridx ());
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	45 A.x = const_cast<double *>(a.data ());
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	46 A.nz = -1;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	47 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	48 #if defined(CS_VER) && (CS_VER >= 2)
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	49 S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	50 #else
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	51 S = CXSPARSE_DNAME (_sqr) (&A, order - 1, 1);
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	52 #endif
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	53
5648 69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	54 N = CXSPARSE_DNAME (_qr) (&A, S);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	55 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	56 if (!N)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	57 (*current_liboctave_error_handler)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	58 ("SparseQR: sparse matrix QR factorization filled");
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	59 count = 1;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	60 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	61 (*current_liboctave_error_handler)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	62 ("SparseQR: sparse matrix QR factorization not implemented");
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	63 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	64 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	65
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	66 SparseQR::SparseQR_rep::~SparseQR_rep (void)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	67 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	68 #ifdef HAVE_CXSPARSE
5648 69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	69 CXSPARSE_DNAME (_sfree) (S);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	70 CXSPARSE_DNAME (_nfree) (N);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	71 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	72 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	73
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	74 SparseMatrix
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	75 SparseQR::SparseQR_rep::V (void) const
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	76 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	77 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	78 // Drop zeros from V and sort
5775 ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe] jwe parents: 5681 diff changeset	79 // FIXME Is the double transpose to sort necessary?
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	80 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5648 69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	81 CXSPARSE_DNAME (_dropzeros) (N->L);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	82 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	83 CXSPARSE_DNAME (_spfree) (N->L);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	84 N->L = CXSPARSE_DNAME (_transpose) (D, 1);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	85 CXSPARSE_DNAME (_spfree) (D);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	86 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	87
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	88 octave_idx_type nc = N->L->n;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	89 octave_idx_type nz = N->L->nzmax;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	90 SparseMatrix ret (N->L->m, nc, nz);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	91 for (octave_idx_type j = 0; j < nc+1; j++)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	92 ret.xcidx (j) = N->L->p[j];
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	93 for (octave_idx_type j = 0; j < nz; j++)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	94 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	95 ret.xridx (j) = N->L->i[j];
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	96 ret.xdata (j) = N->L->x[j];
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	97 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	98 return ret;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	99 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	100 return SparseMatrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	101 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	102 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	103
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	104 ColumnVector
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	105 SparseQR::SparseQR_rep::Pinv (void) const
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	106 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	107 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	108 ColumnVector ret(N->L->m);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	109 for (octave_idx_type i = 0; i < N->L->m; i++)
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	110 #if defined(CS_VER) && (CS_VER >= 2)
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	111 ret.xelem(i) = S->pinv[i];
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	112 #else
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	113 ret.xelem(i) = S->Pinv[i];
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	114 #endif
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	115 return ret;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	116 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	117 return ColumnVector ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	118 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	119 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	120
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	121 ColumnVector
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	122 SparseQR::SparseQR_rep::P (void) const
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	123 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	124 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	125 ColumnVector ret(N->L->m);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	126 for (octave_idx_type i = 0; i < N->L->m; i++)
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	127 #if defined(CS_VER) && (CS_VER >= 2)
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	128 ret.xelem(S->pinv[i]) = i;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	129 #else
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	130 ret.xelem(S->Pinv[i]) = i;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	131 #endif
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	132 return ret;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	133 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	134 return ColumnVector ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	135 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	136 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	137
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	138 SparseMatrix
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	139 SparseQR::SparseQR_rep::R (const bool econ) const
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	140 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	141 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	142 // Drop zeros from R and sort
5775 ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe] jwe parents: 5681 diff changeset	143 // FIXME Is the double transpose to sort necessary?
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	144 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5648 69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	145 CXSPARSE_DNAME (_dropzeros) (N->U);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	146 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	147 CXSPARSE_DNAME (_spfree) (N->U);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	148 N->U = CXSPARSE_DNAME (_transpose) (D, 1);
69a4f320d95a [project @ 2006-03-08 20:17:37 by dbateman] dbateman parents: 5610 diff changeset	149 CXSPARSE_DNAME (_spfree) (D);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	150 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	151
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	152 octave_idx_type nc = N->U->n;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	153 octave_idx_type nz = N->U->nzmax;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	154
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	155 SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	156
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	157 for (octave_idx_type j = 0; j < nc+1; j++)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	158 ret.xcidx (j) = N->U->p[j];
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	159 for (octave_idx_type j = 0; j < nz; j++)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	160 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	161 ret.xridx (j) = N->U->i[j];
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	162 ret.xdata (j) = N->U->x[j];
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	163 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	164 return ret;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	165 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	166 return SparseMatrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	167 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	168 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	169
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	170 Matrix
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	171 SparseQR::SparseQR_rep::C (const Matrix &b) const
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	172 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	173 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	174 octave_idx_type b_nr = b.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	175 octave_idx_type b_nc = b.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	176 octave_idx_type nc = N->L->n;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	177 octave_idx_type nr = nrows;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	178 const double *bvec = b.fortran_vec();
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	179 Matrix ret (b_nr, b_nc);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	180 double *vec = ret.fortran_vec();
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	181 if (nr < 0 \|\| nc < 0 \|\| nr != b_nr)
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	182 (*current_liboctave_error_handler) ("matrix dimension mismatch");
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	183 else if (nr == 0 \|\| nc == 0 \|\| b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	184 ret = Matrix (nc, b_nc, 0.0);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	185 else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	186 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	187 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	188 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	189 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	190 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	191 for (octave_idx_type i = nr; i < S->m2; i++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	192 buf[i] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	193 volatile octave_idx_type nm = (nr < nc ? nr : nc);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	194 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	195 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	196 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	197 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	198 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	199 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	200 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	201
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	202 for (volatile octave_idx_type i = 0; i < nm; i++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	203 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	204 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	205 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	206 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	207 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	208 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	209 for (octave_idx_type i = 0; i < b_nr; i++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	210 vec[i+idx] = buf[i];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	211 }
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	212 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	213 return ret;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	214 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	215 return Matrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	216 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	217 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	218
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	219 Matrix
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	220 SparseQR::SparseQR_rep::Q (void) const
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	221 {
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	222 #ifdef HAVE_CXSPARSE
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	223 octave_idx_type nc = N->L->n;
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	224 octave_idx_type nr = nrows;
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	225 Matrix ret (nr, nr);
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	226 double *vec = ret.fortran_vec();
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	227 if (nr < 0 \|\| nc < 0)
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	228 (*current_liboctave_error_handler) ("matrix dimension mismatch");
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	229 else if (nr == 0 \|\| nc == 0)
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	230 ret = Matrix (nc, nr, 0.0);
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	231 else
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	232 {
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	233 OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	234 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	235 bvec[i] = 0.;
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	236 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	237 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	238 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	239 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	240 bvec[j] = 1.0;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	241 for (octave_idx_type i = nr; i < S->m2; i++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	242 buf[i] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	243 volatile octave_idx_type nm = (nr < nc ? nr : nc);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	244 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	245 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	246 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	247 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	248 CXSPARSE_DNAME (_ipvec) (nr, S->Pinv, bvec, buf);
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	249 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	250 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	251
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	252 for (volatile octave_idx_type i = 0; i < nm; i++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	253 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	254 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	255 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	256 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	257 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	258 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	259 for (octave_idx_type i = 0; i < nr; i++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	260 vec[i+idx] = buf[i];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	261 bvec[j] = 0.0;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	262 }
7505 f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	263 }
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	264 return ret.transpose ();
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	265 #else
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	266 return Matrix ();
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	267 #endif
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	268 }
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	269
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	270 Matrix
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	271 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	272 {
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	273 info = -1;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	274 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	275 octave_idx_type nr = a.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	276 octave_idx_type nc = a.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	277 octave_idx_type b_nc = b.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	278 octave_idx_type b_nr = b.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	279 const double *bvec = b.fortran_vec();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	280 Matrix x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	281
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	282 if (nr < 0 \|\| nc < 0 \|\| nr != b_nr)
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	283 (*current_liboctave_error_handler)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	284 ("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	285 else if (nr == 0 \|\| nc == 0 \|\| b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	286 x = Matrix (nc, b_nc, 0.0);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	287 else if (nr >= nc)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	288 {
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	289 SparseQR q (a, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	290 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	291 return Matrix();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	292 x.resize(nc, b_nc);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	293 double *vec = x.fortran_vec();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	294 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	295 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	296 i++, idx+=nc, bidx+=b_nr)
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 for (octave_idx_type j = nr; j < q.S()->m2; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	300 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	301 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	302 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	303 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	304 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	305 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	306 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	307 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	308 for (volatile octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	309 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	310 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	311 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	312 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	313 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	314 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	315 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	316 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	317 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	318 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, vec + idx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	319 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	320 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	321 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	322 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	323 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	324 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	325 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	326 else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	327 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	328 SparseMatrix at = a.hermitian();
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	329 SparseQR q (at, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	330 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	331 return Matrix();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	332 x.resize(nc, b_nc);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	333 double *vec = x.fortran_vec();
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	334 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	335 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	336 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	337 i++, idx+=nc, bidx+=b_nr)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	338 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	339 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	340 for (octave_idx_type j = nr; j < nbuf; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	341 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	342 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	343 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	344 CXSPARSE_DNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	345 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	346 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	347 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	348 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	349 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	350 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	351 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	352 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	353 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	354 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	355 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	356 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	357 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	358 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	359 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	360 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	361 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	362 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	363 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	364 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	365 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	366 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	367
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	368 return x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	369 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	370 return Matrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	371 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	372 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	373
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	374 SparseMatrix
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	375 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	376 {
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	377 info = -1;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	378 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	379 octave_idx_type nr = a.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	380 octave_idx_type nc = a.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	381 octave_idx_type b_nr = b.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	382 octave_idx_type b_nc = b.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	383 SparseMatrix x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	384 volatile octave_idx_type ii, x_nz;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	385
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	386 if (nr < 0 \|\| nc < 0 \|\| nr != b_nr)
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	387 (*current_liboctave_error_handler)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	388 ("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	389 else if (nr == 0 \|\| nc == 0 \|\| b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	390 x = SparseMatrix (nc, b_nc);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	391 else if (nr >= nc)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	392 {
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	393 SparseQR q (a, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	394 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	395 return SparseMatrix();
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	396 x = SparseMatrix (nc, b_nc, b.nnz());
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	397 x.xcidx(0) = 0;
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	398 x_nz = b.nnz();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	399 ii = 0;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	400 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	401 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	402 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	403 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	404 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	405 for (octave_idx_type j = 0; j < b_nr; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	406 Xx[j] = b.xelem(j,i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	407 for (octave_idx_type j = nr; j < q.S()->m2; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	408 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	409 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	410 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	411 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	412 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	413 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	414 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	415 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	416 for (volatile octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	417 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	418 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	419 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	420 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	421 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	422 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	423 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	424 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	425 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	426 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	427 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	428 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	429 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	430 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	431
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	432 for (octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	433 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	434 double tmp = Xx[j];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	435 if (tmp != 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 if (ii == x_nz)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	438 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	439 // Resize the sparse matrix
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	440 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	441 sz = (sz > 10 ? sz : 10) + x_nz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	442 x.change_capacity (sz);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	443 x_nz = sz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	444 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	445 x.xdata(ii) = tmp;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	446 x.xridx(ii++) = j;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	447 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	448 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	449 x.xcidx(i+1) = ii;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	450 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	451 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	452 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	453 else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	454 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	455 SparseMatrix at = a.hermitian();
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	456 SparseQR q (at, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	457 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	458 return SparseMatrix();
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	459 x = SparseMatrix (nc, b_nc, b.nnz());
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	460 x.xcidx(0) = 0;
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	461 x_nz = b.nnz();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	462 ii = 0;
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	463 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	464 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	465 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	466 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
10314 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 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	469 for (octave_idx_type j = 0; j < b_nr; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	470 Xx[j] = b.xelem(j,i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	471 for (octave_idx_type j = nr; j < nbuf; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	472 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	473 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	474 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	475 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	476 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	477 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	478 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	479 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	480 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	481 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	482 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	483 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	484 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	485 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	486 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	487 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	488 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	489 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	490 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	491 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	492 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	493 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	494 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	495
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	496 for (octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	497 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	498 double tmp = Xx[j];
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	499 if (tmp != 0.0)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	500 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	501 if (ii == x_nz)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	502 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	503 // Resize the sparse matrix
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	504 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	505 sz = (sz > 10 ? sz : 10) + x_nz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	506 x.change_capacity (sz);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	507 x_nz = sz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	508 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	509 x.xdata(ii) = tmp;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	510 x.xridx(ii++) = j;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	511 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	512 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	513 x.xcidx(i+1) = ii;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	514 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	515 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	516 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	517
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	518 x.maybe_compress ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	519 return x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	520 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	521 return SparseMatrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	522 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	523 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	524
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	525 ComplexMatrix
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	526 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	527 {
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	528 info = -1;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	529 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	530 octave_idx_type nr = a.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	531 octave_idx_type nc = a.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	532 octave_idx_type b_nc = b.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	533 octave_idx_type b_nr = b.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	534 ComplexMatrix x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	535
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	536 if (nr < 0 \|\| nc < 0 \|\| nr != b_nr)
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	537 (*current_liboctave_error_handler)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	538 ("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	539 else if (nr == 0 \|\| nc == 0 \|\| b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	540 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	541 else if (nr >= nc)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	542 {
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	543 SparseQR q (a, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	544 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	545 return ComplexMatrix();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	546 x.resize(nc, b_nc);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	547 Complex *vec = x.fortran_vec();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	548 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	549 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	550 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	551 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	552 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	553 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	554 for (octave_idx_type j = 0; j < b_nr; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	555 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	556 Complex c = b.xelem (j,i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	557 Xx[j] = std::real (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	558 Xz[j] = std::imag (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	559 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	560 for (octave_idx_type j = nr; j < q.S()->m2; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	561 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	562 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	563 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	564 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	565 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	566 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	567 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	568 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	569 for (volatile octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	570 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	571 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	572 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	573 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	574 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	575 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	576 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	577 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	578 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	579 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	580 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	581 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	582 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	583 for (octave_idx_type j = nr; j < q.S()->m2; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	584 buf[j] = 0.;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	585 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	586 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	587 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	588 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	589 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	590 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	591 for (volatile octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	592 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	593 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	594 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	595 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	596 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	597 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	598 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	599 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	600 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	601 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	602 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	603 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	604 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	605 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	606 for (octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	607 vec[j+idx] = Complex (Xx[j], Xz[j]);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	608 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	609 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	610 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	611 else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	612 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	613 SparseMatrix at = a.hermitian();
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	614 SparseQR q (at, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	615 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	616 return ComplexMatrix();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	617 x.resize(nc, b_nc);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	618 Complex *vec = x.fortran_vec();
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	619 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	620 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	621 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	622 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	623 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	624 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	625 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	626 for (octave_idx_type j = 0; j < b_nr; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	627 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	628 Complex c = b.xelem (j,i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	629 Xx[j] = std::real (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	630 Xz[j] = std::imag (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	631 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	632 for (octave_idx_type j = nr; j < nbuf; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	633 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	634 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	635 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	636 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	637 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	638 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	639 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	640 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	641 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	642 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	643 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	644 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	645 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	646 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	647 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	648 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	649 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	650 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	651 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	652 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	653 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	654 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	655 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	656 for (octave_idx_type j = nr; j < nbuf; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	657 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	658 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	659 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	660 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	661 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	662 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	663 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	664 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	666 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	667 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	668 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	669 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	670 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	671 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	672 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	673 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	674 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	675 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	676 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	677 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	678 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	679 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	680 for (octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	681 vec[j+idx] = Complex (Xx[j], Xz[j]);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	682 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	683 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	684 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	685
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	686 return x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	687 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	688 return ComplexMatrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	689 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	690 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	691
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	692 SparseComplexMatrix
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	693 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	694 {
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	695 info = -1;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	696 #ifdef HAVE_CXSPARSE
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	697 octave_idx_type nr = a.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	698 octave_idx_type nc = a.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	699 octave_idx_type b_nr = b.rows();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	700 octave_idx_type b_nc = b.cols();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	701 SparseComplexMatrix x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	702 volatile octave_idx_type ii, x_nz;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	703
6924 be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	704 if (nr < 0 \|\| nc < 0 \|\| nr != b_nr)
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	705 (*current_liboctave_error_handler)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	706 ("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	707 else if (nr == 0 \|\| nc == 0 \|\| b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman] dbateman parents: 5797 diff changeset	708 x = SparseComplexMatrix (nc, b_nc);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	709 else if (nr >= nc)
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	710 {
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	711 SparseQR q (a, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	712 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	713 return SparseComplexMatrix();
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	714 x = SparseComplexMatrix (nc, b_nc, b.nnz());
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	715 x.xcidx(0) = 0;
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	716 x_nz = b.nnz();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	717 ii = 0;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	718 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	719 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	720 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	721 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	722 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	723 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	724 for (octave_idx_type j = 0; j < b_nr; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	725 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	726 Complex c = b.xelem (j,i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	727 Xx[j] = std::real (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	728 Xz[j] = std::imag (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	729 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	730 for (octave_idx_type j = nr; j < q.S()->m2; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	731 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	732 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	733 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	734 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	735 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	736 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	737 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	738 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	739 for (volatile octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	740 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	741 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	742 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	743 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	744 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	745 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	746 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	747 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	748 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	749 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	750 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	751 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	752 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	753 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	754 for (octave_idx_type j = nr; j < q.S()->m2; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	755 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	756 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	757 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	758 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	759 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	760 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	761 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	762 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	763 for (volatile octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	764 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	765 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	766 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	767 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	768 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	769 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	770 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	771 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	772 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	773 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	774 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	775 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	776 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	777 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	778
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	779 for (octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	780 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	781 Complex tmp = Complex (Xx[j], Xz[j]);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	782 if (tmp != 0.0)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	783 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	784 if (ii == x_nz)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	785 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	786 // Resize the sparse matrix
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	787 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	788 sz = (sz > 10 ? sz : 10) + x_nz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	789 x.change_capacity (sz);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	790 x_nz = sz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	791 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	792 x.xdata(ii) = tmp;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	793 x.xridx(ii++) = j;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	794 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	795 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	796 x.xcidx(i+1) = ii;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	797 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	798 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	799 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	800 else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	801 {
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	802 SparseMatrix at = a.hermitian();
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	803 SparseQR q (at, 3);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	804 if (! q.ok ())
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	805 return SparseComplexMatrix();
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	806 x = SparseComplexMatrix (nc, b_nc, b.nnz());
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	807 x.xcidx(0) = 0;
10506 bdf5d85cfc5e replace nzmax by nnz where appropriate in liboctave Jaroslav Hajek <highegg@gmail.com> parents: 10350 diff changeset	808 x_nz = b.nnz();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	809 ii = 0;
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	810 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	811 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	812 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
5681 233d98d95659 [project @ 2006-03-16 17:48:55 by dbateman] dbateman parents: 5648 diff changeset	813 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	814 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	815 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	816 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	817 for (octave_idx_type j = 0; j < b_nr; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	818 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	819 Complex c = b.xelem (j,i);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	820 Xx[j] = std::real (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	821 Xz[j] = std::imag (c);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	822 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	823 for (octave_idx_type j = nr; j < nbuf; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	824 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	825 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	826 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	827 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	828 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	829 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	830 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	831 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	832 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	833 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	834 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	835 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	836 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	837 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	838 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	839 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	840 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	841 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	842 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	843 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	844 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	845 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	846 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	847 for (octave_idx_type j = nr; j < nbuf; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	848 buf[j] = 0.;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	849 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	850 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	851 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	852 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	853 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	854 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	855 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	856 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	857 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	858 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	859 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	860 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	861 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	862 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	863 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	864 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	865 #if defined(CS_VER) && (CS_VER >= 2)
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	866 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	867 #else
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	868 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz);
5792 eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman] dbateman parents: 5775 diff changeset	869 #endif
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	870 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	871
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	872 for (octave_idx_type j = 0; j < nc; j++)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	873 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	874 Complex tmp = Complex (Xx[j], Xz[j]);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	875 if (tmp != 0.0)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	876 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	877 if (ii == x_nz)
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	878 {
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	879 // Resize the sparse matrix
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	880 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	881 sz = (sz > 10 ? sz : 10) + x_nz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	882 x.change_capacity (sz);
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	883 x_nz = sz;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	884 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	885 x.xdata(ii) = tmp;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	886 x.xridx(ii++) = j;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	887 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	888 }
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	889 x.xcidx(i+1) = ii;
07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	890 }
5797 11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman] dbateman parents: 5792 diff changeset	891 info = 0;
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	892 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	893
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	894 x.maybe_compress ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	895 return x;
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	896 #else
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	897 return SparseComplexMatrix ();
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	898 #endif
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	899 }
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	900
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	901 Matrix
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	902 qrsolve(const SparseMatrix &a, const MArray<double> &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	903 octave_idx_type &info)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	904 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	905 return qrsolve (a, Matrix (b), info);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	906 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	907
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	908 ComplexMatrix
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10314 diff changeset	909 qrsolve(const SparseMatrix &a, const MArray<Complex> &b,
10314 07ebe522dac2 untabify liboctave C++ sources John W. Eaton <jwe@octave.org> parents: 10158 diff changeset	910 octave_idx_type &info)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	911 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	912 return qrsolve (a, ComplexMatrix (b), info);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7505 diff changeset	913 }

Mercurial > octave

annotate liboctave/SparseQR.cc @ 11523:fd0a3ac60b0e