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