24 #if defined (HAVE_CONFIG_H)
37 template <
typename SPARSE_T>
48 #if defined (HAVE_CXSPARSE)
52 typedef void symbolic_type;
53 typedef void numeric_type;
62 #if defined (HAVE_CXSPARSE)
66 typedef void symbolic_type;
67 typedef void numeric_type;
71 template <
typename SPARSE_T>
82 #if defined (HAVE_CXSPARSE)
89 SPARSE_T
V (
void)
const;
95 SPARSE_T
R (
bool econ)
const;
97 typename SPARSE_T::dense_matrix_type
98 C (
const typename SPARSE_T::dense_matrix_type&
b)
const;
100 typename SPARSE_T::dense_matrix_type
111 template <
typename RHS_T,
typename RET_T>
115 template <
typename RHS_T,
typename RET_T>
128 template <
typename SPARSE_T>
132 #if defined (HAVE_CXSPARSE)
148 template <
typename SPARSE_T>
152 #if defined (HAVE_CXSPARSE)
175 : count (1), nrows (a.
rows ()), ncols (a.
columns ())
180 #if defined (HAVE_CXSPARSE)
191 A.x =
const_cast<double *
>(a.
data ());
200 (*current_liboctave_error_handler)
201 (
"sparse_qr: sparse matrix QR factorization filled");
207 octave_unused_parameter (order);
209 (*current_liboctave_error_handler)
210 (
"sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
218 #if defined (HAVE_CXSPARSE)
228 #if defined (HAVE_CXSPARSE)
246 ret.
xcidx (j) =
N->L->p[j];
250 ret.
xridx (j) =
N->L->i[j];
251 ret.
xdata (j) =
N->L->x[j];
267 #if defined (HAVE_CXSPARSE)
283 SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
286 ret.
xcidx (j) =
N->U->p[j];
290 ret.
xridx (j) =
N->U->i[j];
291 ret.
xdata (j) =
N->U->x[j];
298 octave_unused_parameter (econ);
309 #if defined (HAVE_CXSPARSE)
322 if (nr < 0 || nc < 0 || nr != b_nr)
323 (*current_liboctave_error_handler) (
"matrix dimension mismatch");
325 if (nr == 0 || nc == 0 || b_nc == 0)
326 ret =
Matrix (nc, b_nc, 0.0);
362 octave_unused_parameter (b);
373 #if defined (HAVE_CXSPARSE)
379 if (nr < 0 || nc < 0)
380 (*current_liboctave_error_handler) (
"matrix dimension mismatch");
382 if (nr == 0 || nc == 0)
383 ret =
Matrix (nc, nr, 0.0);
440 #if defined (HAVE_CXSPARSE)
448 const double *bvec = b.
data ();
456 i++, idx+=nc, bidx+=
b_nr)
488 octave_unused_parameter (b);
503 #if defined (HAVE_CXSPARSE)
514 const double *bvec = b.
data ();
524 i++, idx+=nc, bidx+=
b_nr)
556 octave_unused_parameter (b);
571 #if defined (HAVE_CXSPARSE)
579 SparseMatrix
x (nc, b_nc, b.
nnz ());
626 sz = (sz > 10 ? sz : 10) + x_nz;
627 x.change_capacity (sz);
645 octave_unused_parameter (b);
647 return SparseMatrix ();
660 #if defined (HAVE_CXSPARSE)
671 SparseMatrix
x (nc, b_nc, b.
nnz ());
719 sz = (sz > 10 ? sz : 10) + x_nz;
720 x.change_capacity (sz);
740 octave_unused_parameter (b);
742 return SparseMatrix ();
755 #if defined (HAVE_CXSPARSE)
822 vec[j+idx] =
Complex (Xx[j], Xz[j]);
831 octave_unused_parameter (b);
846 #if defined (HAVE_CXSPARSE)
920 vec[j+idx] =
Complex (Xx[j], Xz[j]);
929 octave_unused_parameter (b);
941 : count (1), nrows (a.
rows ()), ncols (a.
columns ())
946 #if defined (HAVE_CXSPARSE)
957 A.x =
const_cast<cs_complex_t *
> (
958 reinterpret_cast<const cs_complex_t *
> (a.
data ()));
967 (*current_liboctave_error_handler)
968 (
"sparse_qr: sparse matrix QR factorization filled");
974 octave_unused_parameter (order);
976 (*current_liboctave_error_handler)
977 (
"sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
985 #if defined (HAVE_CXSPARSE)
995 #if defined (HAVE_CXSPARSE)
1012 ret.
xcidx (j) =
N->L->p[j];
1016 ret.
xridx (j) =
N->L->i[j];
1033 #if defined (HAVE_CXSPARSE)
1052 ret.
xcidx (j) =
N->U->p[j];
1056 ret.
xridx (j) =
N->U->i[j];
1064 octave_unused_parameter (econ);
1075 #if defined (HAVE_CXSPARSE)
1080 const cs_complex_t *bvec =
reinterpret_cast<const cs_complex_t *
> (b.
fortran_vec ());
1084 if (nr < 0 || nc < 0 || nr != b_nr)
1085 (*current_liboctave_error_handler) (
"matrix dimension mismatch");
1087 if (nr == 0 || nc == 0 || b_nc == 0)
1103 reinterpret_cast<cs_complex_t *
> (buf),
1113 reinterpret_cast<cs_complex_t *
> (buf));
1118 vec[
i+idx] = buf[
i];
1126 octave_unused_parameter (b);
1137 #if defined (HAVE_CXSPARSE)
1143 if (nr < 0 || nc < 0)
1144 (*current_liboctave_error_handler) (
"matrix dimension mismatch");
1146 if (nr == 0 || nc == 0)
1153 bvec[
i] = cs_complex_t (0.0, 0.0);
1161 bvec[j] = cs_complex_t (1.0, 0.0);
1167 reinterpret_cast<cs_complex_t *
> (buf),
1177 reinterpret_cast<cs_complex_t *
> (buf));
1182 vec[
i+idx] = buf[
i];
1184 bvec[j] = cs_complex_t (0.0, 0.0);
1201 SparseComplexMatrix>
1206 #if defined (HAVE_CXSPARSE)
1214 SparseComplexMatrix
x (nc, b_nc, b.
nnz ());
1287 sz = (sz > 10 ? sz : 10) + x_nz;
1288 x.change_capacity (sz);
1306 octave_unused_parameter (b);
1308 return SparseComplexMatrix ();
1317 SparseComplexMatrix>
1322 #if defined (HAVE_CXSPARSE)
1333 SparseComplexMatrix
x (nc, b_nc, b.
nnz ());
1407 sz = (sz > 10 ? sz : 10) + x_nz;
1408 x.change_capacity (sz);
1422 x.maybe_compress ();
1428 octave_unused_parameter (b);
1430 return SparseComplexMatrix ();
1444 #if defined (HAVE_CXSPARSE)
1453 cs_complex_t *vec =
reinterpret_cast<cs_complex_t *
> (x.
fortran_vec ());
1466 buf[j] = cs_complex_t (0.0, 0.0);
1470 reinterpret_cast<cs_complex_t *
>(Xx),
1495 octave_unused_parameter (b);
1511 #if defined (HAVE_CXSPARSE)
1523 cs_complex_t *vec =
reinterpret_cast<cs_complex_t *
> (x.
fortran_vec ());
1542 buf[j] = cs_complex_t (0.0, 0.0);
1545 CXSPARSE_ZNAME (_pvec) (S->q,
reinterpret_cast<cs_complex_t *
>(Xx),
1570 octave_unused_parameter (b);
1581 SparseComplexMatrix>
1586 #if defined (HAVE_CXSPARSE)
1594 SparseComplexMatrix
x (nc, b_nc, b.
nnz ());
1611 buf[j] = cs_complex_t (0.0, 0.0);
1615 reinterpret_cast<cs_complex_t *
>(Xx),
1631 reinterpret_cast<cs_complex_t *
> (Xx),
1645 sz = (sz > 10 ? sz : 10) + x_nz;
1646 x.change_capacity (sz);
1660 x.maybe_compress ();
1666 octave_unused_parameter (b);
1668 return SparseComplexMatrix ();
1677 SparseComplexMatrix>
1682 #if defined (HAVE_CXSPARSE)
1693 SparseComplexMatrix
x (nc, b_nc, b.
nnz ());
1715 buf[j] = cs_complex_t (0.0, 0.0);
1719 reinterpret_cast<cs_complex_t *
>(Xx),
1735 reinterpret_cast<cs_complex_t *
> (Xx),
1749 sz = (sz > 10 ? sz : 10) + x_nz;
1750 x.change_capacity (sz);
1764 x.maybe_compress ();
1770 octave_unused_parameter (b);
1772 return SparseComplexMatrix ();
1786 #if defined (HAVE_CXSPARSE)
1794 const cs_complex_t *bvec =
reinterpret_cast<const cs_complex_t *
>
1798 cs_complex_t *vec =
reinterpret_cast<cs_complex_t *
>
1804 i++, idx+=nc, bidx+=
b_nr)
1809 buf[j] = cs_complex_t (0.0, 0.0);
1836 octave_unused_parameter (b);
1852 #if defined (HAVE_CXSPARSE)
1863 const cs_complex_t *bvec =
reinterpret_cast<const cs_complex_t *
>
1867 cs_complex_t *vec =
reinterpret_cast<cs_complex_t *
> (x.
fortran_vec ());
1878 i++, idx+=nc, bidx+=
b_nr)
1883 buf[j] = cs_complex_t (0.0, 0.0);
1910 octave_unused_parameter (b);
1925 #if defined (HAVE_CXSPARSE)
1933 SparseComplexMatrix
x (nc, b_nc, b.
nnz ());
1950 buf[j] = cs_complex_t (0.0, 0.0);
1954 reinterpret_cast<cs_complex_t *
>(Xx),
1970 reinterpret_cast<cs_complex_t *
> (Xx),
1984 sz = (sz > 10 ? sz : 10) + x_nz;
1985 x.change_capacity (sz);
1999 x.maybe_compress ();
2005 octave_unused_parameter (b);
2007 return SparseComplexMatrix ();
2020 #if defined (HAVE_CXSPARSE)
2031 SparseComplexMatrix
x (nc, b_nc, b.
nnz ());
2053 buf[j] = cs_complex_t (0.0, 0.0);
2056 CXSPARSE_ZNAME (_pvec) (S->q,
reinterpret_cast<cs_complex_t *
>(Xx),
2072 reinterpret_cast<cs_complex_t *
>(Xx), nc);
2085 sz = (sz > 10 ? sz : 10) + x_nz;
2086 x.change_capacity (sz);
2100 x.maybe_compress ();
2106 octave_unused_parameter (b);
2108 return SparseComplexMatrix ();
2113 template <
typename SPARSE_T>
2118 template <
typename SPARSE_T>
2123 template <
typename SPARSE_T>
2130 template <
typename SPARSE_T>
2133 if (--rep->count == 0)
2137 template <
typename SPARSE_T>
2143 if (--rep->count == 0)
2153 template <
typename SPARSE_T>
2160 template <
typename SPARSE_T>
2167 template <
typename SPARSE_T>
2174 template <
typename SPARSE_T>
2181 template <
typename SPARSE_T>
2185 return rep->R (econ);
2188 template <
typename SPARSE_T>
2189 typename SPARSE_T::dense_matrix_type
2195 template <
typename SPARSE_T>
2196 typename SPARSE_T::dense_matrix_type
2207 template <
typename SPARSE_T>
2212 enum { order = -1 };
2231 template <
typename SPARSE_T>
2232 template <
typename RHS_T,
typename RET_T>
2247 if (nr < 0 || nc < 0 || nr != b_nr)
2248 (*current_liboctave_error_handler)
2249 (
"matrix dimension mismatch in solution of minimum norm problem");
2251 if (nr == 0 || nc == 0 || b_nc == 0)
2255 return RET_T (nc, b_nc, 0.0);
2261 return q.
ok () ? q.
tall_solve<RHS_T, RET_T> (
b, info) : RET_T ();
2267 return q.
ok () ? q.wide_solve<RHS_T, RET_T> (
b, info) : RET_T ();
2271 template <
typename SPARSE_T>
2272 template <
typename RHS_T,
typename RET_T>
2276 return rep->template tall_solve<RHS_T, RET_T> (
b, info);
2279 template <
typename SPARSE_T>
2280 template <
typename RHS_T,
typename RET_T>
2284 return rep->template wide_solve<RHS_T, RET_T> (
b, info);
2288 qrsolve (
const SparseMatrix& a,
const MArray<double>& b,
2296 qrsolve (
const SparseMatrix& a,
const SparseMatrix& b,
2304 qrsolve (
const SparseMatrix& a,
const MArray<Complex>& b,
2312 qrsolve (
const SparseMatrix& a,
const SparseComplexMatrix& b,
2316 SparseComplexMatrix> (
a,
b, info);
2320 qrsolve (
const SparseComplexMatrix& a,
const MArray<double>& b,
2328 qrsolve (
const SparseComplexMatrix& a,
const SparseMatrix& b,
2332 SparseComplexMatrix>
2337 qrsolve (
const SparseComplexMatrix& a,
const MArray<Complex>& b,
2345 qrsolve (
const SparseComplexMatrix& a,
const SparseComplexMatrix& b,
2349 SparseComplexMatrix>
2355 template class sparse_qr<SparseMatrix>;
2357 template class sparse_qr<SparseComplexMatrix>;
octave_idx_type * xridx(void)
octave_idx_type cols(void) const
Octave interface to the compression and uncompression libraries.
for(octave_idx_type n=0;n< hcv.numel();n++)
SPARSE_T R(bool econ) const
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
octave_idx_type rows(void) const
octave_refcount< int > count
T & xelem(octave_idx_type n)
SPARSE_T::dense_matrix_type Q(void) const
octave_idx_type * xcidx(void)
sparse_qr & operator=(const sparse_qr &a)
octave_idx_type * cidx(void)
octave_idx_type columns(void) const
octave_idx_type rows(void) const
octave_idx_type nnz(void) const
Actual number of nonzero terms.
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
SPARSE_T::dense_matrix_type Q(void) const
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
OCTAVE_EXPORT octave_value_list isdir nd deftypefn *std::string nm
ComplexMatrix hermitian(void) const
ColumnVector P(void) const
const T * data(void) const
F77_RET_T const F77_INT & N
Matrix transpose(void) const
SPARSE_T R(bool econ=false) const
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
T & xelem(octave_idx_type n)
cxsparse_types< SPARSE_T >::numeric_type * N
#define CXSPARSE_DNAME(name)
defaults to zero A value of zero computes the digamma a value the trigamma and so on The digamma function is defined
octave_idx_type * ridx(void)
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
static RET_T solve(const SPARSE_T &a, const RHS_T &b, octave_idx_type &info)
=val(i)}if ode{val(i)}occurs in table i
RET_T tall_solve(const RHS_T &b, octave_idx_type &info) const
cxsparse_types< SPARSE_T >::symbolic_type * S
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
ColumnVector Pinv(void) const
sparse_qr_rep(const SPARSE_T &a, int order)
ColumnVector Pinv(void) const
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
ColumnVector P(void) const
std::complex< double > Complex
const T * fortran_vec(void) const
SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type &b) const
octave_idx_type cols(void) const
SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type &b) const
#define CXSPARSE_ZNAME(name)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
RET_T wide_solve(const RHS_T &b, octave_idx_type &info) const
F77_RET_T const F77_INT F77_CMPLX * A