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