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