octave-nkf: liboctave/SparseQR.cc annotate

annotate liboctave/SparseQR.cc @ 10506:bdf5d85cfc5e

replace nzmax by nnz where appropriate in liboctave

author	Jaroslav Hajek <highegg@gmail.com>
date	Sat, 10 Apr 2010 21:04:55 +0200
parents	12884915a8e4
children	367bfee35ba0

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

Mercurial > hg > octave-nkf

annotate liboctave/SparseQR.cc @ 10506:bdf5d85cfc5e