26 #ifndef _CXSC_SRMATRIX_HPP_INCLUDED 27 #define _CXSC_SRMATRIX_HPP_INCLUDED 30 #include <rmatrix.hpp> 31 #include <srvector.hpp> 36 #include <sparsedot.hpp> 37 #include <sparsematrix.hpp> 54 class scimatrix_slice;
56 class scivector_slice;
59 inline bool comp_pair(std::pair<int,real> p1, std::pair<int,real> p2) {
60 return p1.first < p2.first;
115 const std::vector<real>&
values()
const {
123 lb1 = lb2 = ub1 = ub2 = 0;
127 srmatrix(
const int r,
const int c) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
128 p = std::vector<int>((n>0) ? n+1 : 1, 0);
129 ind.reserve(2*(m+n));
136 srmatrix(
const int r,
const int c,
const int e) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
137 p = std::vector<int>((n>0) ? n+1 : 1, 0);
155 p = std::vector<int>(n+1,0);
161 std::vector<triplet_store<real> > work;
164 for(
int k=0 ; k<nnz ; k++) {
165 work.push_back(triplet_store<real>(rows[
Lb(rows)+k],cols[
Lb(cols)+k],
values[
Lb(
values)+k]));
168 sort(work.begin(), work.end());
172 for(
int j=0 ; j<n ; j++) {
174 while((
unsigned int)i < work.size() && work[i].col == j ) {
175 ind.push_back(work[i].row);
176 x.push_back(work[i].val);
183 }
else if(t == compressed_row) {
187 p = std::vector<int>(n+1,0);
193 for(
int i=0 ; i<n+1 ; i++)
194 p[i] = rows[
Lb(rows)+i];
196 std::vector<triplet_store<real> > work;
199 for(
int j=0 ; j<n ; j++) {
200 for(
int k=p[j] ; k<p[j+1] ; k++) {
201 work.push_back(triplet_store<real>(j,cols[
Lb(cols)+k],
values[
Lb(
values)+k]));
205 sort(work.begin(), work.end());
209 for(
int j=0 ; j<n ; j++) {
211 while((
unsigned int)i < work.size() && work[i].col == j ) {
212 ind.push_back(work[i].row);
213 x.push_back(work[i].val);
220 }
else if(t == compressed_column) {
223 p = std::vector<int>(n+1,0);
229 for(
int i=0 ; i<n+1 ; i++)
230 p[i] = rows[
Lb(rows)+i];
232 std::vector<std::pair<int,real> > work;
235 for(
int j=0 ; j<n ; j++) {
238 for(
int k=p[j] ; k<p[j+1] ; k++) {
242 std::sort(work.begin(),work.end(),comp_pair);
244 for(
unsigned int i=0 ; i<work.size() ; i++) {
245 ind.push_back(work[i].first);
246 x.push_back(work[i].second);
267 p = std::vector<int>(n+1,0);
273 std::vector<triplet_store<real> > work;
276 for(
int k=0 ; k<nnz ; k++) {
277 work.push_back(triplet_store<real>(rows[k],cols[k],
values[k]));
280 sort(work.begin(), work.end());
284 for(
int j=0 ; j<n ; j++) {
286 while((
unsigned int)i < work.size() && work[i].col == j ) {
287 ind.push_back(work[i].row);
288 x.push_back(work[i].val);
295 }
else if(t == compressed_row) {
299 p = std::vector<int>(n+1,0);
305 for(
int i=0 ; i<n+1 ; i++)
308 std::vector<triplet_store<real> > work;
311 for(
int j=0 ; j<n ; j++) {
312 for(
int k=p[j] ; k<p[j+1] ; k++) {
313 work.push_back(triplet_store<real>(j,cols[k],
values[k]));
317 sort(work.begin(), work.end());
321 for(
int j=0 ; j<n ; j++) {
323 while((
unsigned int)i < work.size() && work[i].col == j ) {
324 ind.push_back(work[i].row);
325 x.push_back(work[i].val);
332 }
else if(t == compressed_column) {
335 p = std::vector<int>(n+1,0);
341 for(
int i=0 ; i<n+1 ; i++)
344 std::vector<std::pair<int,real> > work;
347 for(
int j=0 ; j<n ; j++) {
350 for(
int k=p[j] ; k<p[j+1] ; k++) {
351 work.push_back(std::make_pair(cols[k],
values[k]));
354 std::sort(work.begin(),work.end(),comp_pair);
356 for(
unsigned int i=0 ; i<work.size() ; i++) {
357 ind.push_back(work[i].first);
358 x.push_back(work[i].second);
368 p = std::vector<int>((n>0) ? n+1 : 1, 0);
369 ind.reserve((m*n*0.1 < 2*m) ? (
int)(m*n*0.1) : 2*m);
370 x.reserve((m*n*0.1 < 2*m) ? (
int)(m*n*0.1) : 2*m);
375 for(
int j=0 ; j<n ; j++) {
376 for(
int i=0 ; i<m ; i++) {
377 if(A[i+lb1][j+lb2] != 0.0) {
379 x.push_back(A[i+lb1][j+lb2]);
394 srmatrix(
const int ms,
const int ns,
const rmatrix& A) : m(ms), n(ns), lb1(1), ub1(ms), lb2(1), ub2(ns) {
397 p = std::vector<int>((n>0) ? n+1 : 1, 0);
401 std::vector<triplet_store<real> > work;
405 for(
int i=0 ; i<
ColLen(A) ; i++) {
406 for(
int j=
Lb(A,2) ; j<=
Ub(A,2) ; j++) {
407 if(i+j >=0 && i+j < n) {
408 work.push_back(triplet_store<real>(i,i+j,A[i+
Lb(A,1)][j]));
413 sort(work.begin(), work.end());
417 for(
int j=0 ; j<n ; j++) {
419 while((
unsigned int)i < work.size() && work[i].col == j ) {
420 ind.push_back(work[i].row);
421 x.push_back(work[i].val);
437 for(
int j=0 ; j<n ; j++) {
438 for(
int k=p[j] ; k<p[j+1] ; k++) {
439 A[ind[k]+lb1][j+lb2] = x[k];
450 std::vector<int> pnew(n+1,0);
451 std::vector<int> indnew;
452 std::vector<real> xnew;
455 for(
int j=0 ; j<n ; j++) {
456 for(
int k=p[j] ; k<p[j+1] ; k++) {
458 xnew.push_back(x[k]);
459 indnew.push_back(ind[k]);
474 return sp_ms_assign<srmatrix,real,real>(*
this,A);
479 return spf_mm_assign<srmatrix,rmatrix,real>(*
this,A);
484 return spf_mm_assign<srmatrix,rmatrix,real>(*
this,A);
505 #if(CXSC_INDEX_CHECK) 506 if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
507 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator()(int, int)"));
510 for(
int k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
511 if(ind[k] == i-lb1) r = x[k];
526 #if(CXSC_INDEX_CHECK) 527 if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
528 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::element(int, int)"));
531 for(k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
532 if(ind[k] == i-lb1)
return x[k];
536 std::vector<int>::iterator ind_it = ind.begin() + k;
537 std::vector<real>::iterator x_it = x.begin() + k;
538 ind.insert(ind_it, i-lb1);
539 x_it = x.insert(x_it, 0.0);
540 for(k=j-lb2+1 ; k<(int)p.size() ; k++)
566 for(
int k=0 ; k<n ; k++) {
569 std::map<int,real> work;
570 for(
int j=p[q[
Lb(q)+k]] ; j<p[q[
Lb(q)+k]+1] ; j++)
571 work.insert(std::make_pair(per[
Lb(per)+ind[j]], x[j]));
573 for(std::map<int,real>::iterator it = work.begin() ; it != work.end() ; it++) {
574 A.ind.push_back(it->first);
575 A.x.push_back(it->second);
592 for(
int k=0 ; k<n ; k++) {
595 std::map<int,real> work;
596 for(
int j=p[k] ; j<p[k+1] ; j++)
597 work.insert(std::make_pair(per[
Lb(per)+ind[j]], x[j]));
599 for(std::map<int,real>::iterator it = work.begin() ; it != work.end() ; it++) {
600 A.ind.push_back(it->first);
601 A.x.push_back(it->second);
626 return p[n]/((double)m*n);
636 return spf_mm_addassign<srmatrix,rmatrix,rmatrix>(*
this,B);
641 return spf_mm_addassign<srmatrix,rmatrix_slice,rmatrix>(*
this,B);
646 return spsp_mm_addassign<srmatrix,srmatrix,real>(*
this,B);
651 return spf_mm_subassign<srmatrix,rmatrix,rmatrix>(*
this,B);
656 return spf_mm_subassign<srmatrix,rmatrix_slice,rmatrix>(*
this,B);
661 return spsp_mm_subassign<srmatrix,srmatrix,real>(*
this,B);
666 return spf_mm_multassign<srmatrix,rmatrix,sparse_dot,rmatrix>(*
this,B);
671 return spf_mm_multassign<srmatrix,rmatrix,sparse_dot,rmatrix>(*
this,B);
676 return spsp_mm_multassign<srmatrix,srmatrix,sparse_dot,real>(*
this,B);
681 return sp_ms_multassign(*
this,r);
686 return sp_ms_divassign(*
this,r);
738 #include "matrix_friend_declarations.inl" 742 dat =
new real[A.m*A.n];
743 lb1 = A.lb1; lb2 = A.lb2; ub1 = A.ub1; ub2 = A.ub2;
748 for(
int j=0 ; j<A.n ; j++) {
749 for(
int k=A.p[j] ; k<A.p[j+1] ; k++) {
750 dat[A.ind[k]*A.n+j] = A.x[k];
757 srmatrix I(A.m, A.n, (A.m>A.n) ? A.m : A.n);
758 I.lb1 = A.lb1; I.lb2 = A.lb2;
759 I.ub1 = A.ub1; I.ub2 = A.ub2;
762 for(
int i=0 ; i<A.m ; i++) {
763 I.p[i+1] = I.p[i] + 1;
768 for(
int i=0 ; i<A.n ; i++) {
769 I.p[i+1] = I.p[i] + 1;
783 std::vector<int> w(A.m,0);
784 for(
unsigned int i=0 ; i<A.ind.size() ; i++)
790 for(
unsigned int i=1 ; i<B.p.size() ; i++)
791 B.p[i] = w[i-1] + B.p[i-1];
794 w.insert(w.begin(), 0);
795 for(
unsigned int i=1 ; i<w.size() ; i++) {
803 for(
int j=0 ; j<A.n ; j++) {
804 for(
int k=A.p[j] ; k<A.p[j+1] ; k++) {
817 for(
unsigned int i=0 ; i<ret.x.size() ; i++)
818 ret.x[i] =
abs(ret.x[i]);
831 res.p[A.n] = A.p[A.n];
833 for(
int j=0 ; j<res.n ; j++) {
834 for(
int k=A.p[j] ; k<A.p[j+1] ; k++) {
836 res.x.push_back(
abs(A.x[k]));
838 res.x.push_back(-
abs(A.x[k]));
924 inline void Resize(
srmatrix& A,
const int l1,
const int u1,
const int l2,
const int u2) {
925 sp_m_resize(A,l1,u1,l2,u2);
936 return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(A,B);
947 return spf_mm_mult<srmatrix,rmatrix,rmatrix,sparse_dot>(A,B);
958 return fsp_mm_mult<rmatrix_slice,srmatrix,rmatrix,sparse_dot>(A,B);
969 return spf_mm_mult<srmatrix,rmatrix_slice,rmatrix,sparse_dot>(A,B);
980 return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(A,B);
985 return sp_ms_mult<srmatrix,real,srmatrix>(A,r);
990 return sp_ms_div<srmatrix,real,srmatrix>(A,r);
995 return sp_sm_mult<real,srmatrix,srmatrix>(r,A);
1006 return spf_mv_mult<srmatrix,rvector,rvector,sparse_dot>(A,v);
1017 return spf_mv_mult<srmatrix,rvector_slice,rvector,sparse_dot>(A,v);
1028 return spsp_mv_mult<srmatrix,srvector,srvector,sparse_dot,real>(A,v);
1039 return spsl_mv_mult<srmatrix,srvector_slice,srvector,sparse_dot,real>(A,v);
1050 return fsp_mv_mult<rmatrix,srvector,rvector,sparse_dot>(A,v);
1061 return fsp_mv_mult<rmatrix_slice,srvector,rvector,sparse_dot>(A,v);
1072 return fsl_mv_mult<rmatrix,srvector_slice,rvector,sparse_dot>(A,v);
1083 return fsl_mv_mult<rmatrix_slice,srvector_slice,rvector,sparse_dot>(A,v);
1088 return fsp_mm_add<rmatrix,srmatrix,rmatrix>(A,B);
1093 return spf_mm_add<srmatrix,rmatrix,rmatrix>(A,B);
1098 return fsp_mm_add<rmatrix_slice,srmatrix,rmatrix>(A,B);
1103 return spf_mm_add<srmatrix,rmatrix_slice,rmatrix>(A,B);
1108 return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(A,B);
1113 return fsp_mm_sub<rmatrix,srmatrix,rmatrix>(A,B);
1118 return spf_mm_sub<srmatrix,rmatrix,rmatrix>(A,B);
1123 return fsp_mm_sub<rmatrix_slice,srmatrix,rmatrix>(A,B);
1128 return spf_mm_sub<srmatrix,rmatrix_slice,rmatrix>(A,B);
1133 return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(A,B);
1138 return sp_m_negative<srmatrix,srmatrix>(M);
1157 return fsp_mm_addassign(*
this,B);
1161 return fsp_mm_addassign(*
this,B);
1165 return fsp_mm_subassign(*
this,B);
1169 return fsp_mm_subassign(*
this,B);
1173 return fsp_mm_multassign<rmatrix,srmatrix,sparse_dot,rmatrix>(*
this,B);
1177 return fsp_mm_multassign<rmatrix_slice,srmatrix,sparse_dot,rmatrix>(*
this,B);
1182 return spsp_mm_comp(A,B);
1187 return spf_mm_comp(A,B);
1192 return fsp_mm_comp(A,B);
1197 return fsp_mm_comp(A,B);
1202 return spf_mm_comp(A,B);
1207 return !spsp_mm_comp(A,B);
1212 return !spf_mm_comp(A,B);
1217 return !fsp_mm_comp(A,B);
1222 return !fsp_mm_comp(A,B);
1227 return !spf_mm_comp(A,B);
1232 return spsp_mm_less<srmatrix,srmatrix,real>(A,B);
1237 return spf_mm_less<srmatrix,rmatrix,real>(A,B);
1242 return fsp_mm_less<rmatrix,srmatrix,real>(A,B);
1247 return fsp_mm_less<rmatrix_slice,srmatrix,real>(A,B);
1252 return spf_mm_less<srmatrix,rmatrix_slice,real>(A,B);
1257 return spsp_mm_leq<srmatrix,srmatrix,real>(A,B);
1262 return spf_mm_leq<srmatrix,rmatrix,real>(A,B);
1267 return fsp_mm_leq<rmatrix,srmatrix,real>(A,B);
1272 return fsp_mm_leq<rmatrix_slice,srmatrix,real>(A,B);
1277 return spf_mm_leq<srmatrix,rmatrix_slice,real>(A,B);
1282 return spsp_mm_greater<srmatrix,srmatrix,real>(A,B);
1287 return spf_mm_greater<srmatrix,rmatrix,real>(A,B);
1292 return fsp_mm_greater<rmatrix,srmatrix,real>(A,B);
1297 return fsp_mm_greater<rmatrix_slice,srmatrix,real>(A,B);
1302 return spf_mm_greater<srmatrix,rmatrix_slice,real>(A,B);
1307 return spsp_mm_geq<srmatrix,srmatrix,real>(A,B);
1312 return spf_mm_geq<srmatrix,rmatrix,real>(A,B);
1317 return fsp_mm_geq<rmatrix,srmatrix,real>(A,B);
1322 return fsp_mm_geq<rmatrix_slice,srmatrix,real>(A,B);
1327 return spf_mm_geq<srmatrix,rmatrix_slice,real>(A,B);
1341 inline std::ostream& operator<<(std::ostream& os,
const srmatrix& A) {
1342 return sp_m_output<srmatrix,real>(os,A);
1351 inline std::istream& operator>>(std::istream& is,
srmatrix& A) {
1352 return sp_m_input<srmatrix,real>(is,A);
1375 A.p = std::vector<int>(A.n+1, 0);
1376 A.ind.reserve(A.m + A.n);
1377 A.x.reserve(A.m + A.n);
1379 for(
int i=0 ; i<A.n ; i++) {
1381 for(
int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
1382 if(Mat.ind[j] >= sl1l-Mat.lb1 && Mat.ind[j] <= sl1u-Mat.lb1) {
1383 A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
1384 A.x.push_back(Mat.x[j]);
1403 A.p = std::vector<int>(A.n+1, 0);
1404 A.ind.reserve(A.m + A.n);
1405 A.x.reserve(A.m + A.n);
1407 for(
int i=0 ; i<A.n ; i++) {
1409 for(
int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
1410 if(Mat.ind[j] >= sl1l-Mat.lb1 && Mat.ind[j] <= sl1u-Mat.lb1) {
1411 A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
1412 A.x.push_back(Mat.x[j]);
1425 return sl_ms_assign<srmatrix_slice, real, std::vector<real>::iterator,
real>(*
this,C);
1430 return slsp_mm_assign<srmatrix_slice, srmatrix, std::vector<real>::iterator>(*
this,C);
1435 return slf_mm_assign<srmatrix_slice, rmatrix, std::vector<real>::iterator,
real>(*
this,C);
1440 return slf_mm_assign<srmatrix_slice, rmatrix_slice, std::vector<real>::iterator,
real>(*
this,C);
1539 #if(CXSC_INDEX_CHECK) 1540 if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
1541 cxscthrow(ELEMENT_NOT_IN_VEC(
"srmatrix_slice::operator()(int, int)"));
1553 #if(CXSC_INDEX_CHECK) 1554 if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
1555 cxscthrow(ELEMENT_NOT_IN_VEC(
"srmatrix::element(int, int)"));
1594 #include "matrix_friend_declarations.inl" 1598 dat =
new real[A.A.m*A.A.n];
1599 lb1 = A.A.lb1; lb2 = A.A.lb2; ub1 = A.A.ub1; ub2 = A.A.ub2;
1603 for(
int j=0 ; j<A.A.n ; j++) {
1604 for(
int k=A.A.p[j] ; k<A.A.p[j+1] ; k++) {
1605 dat[A.A.ind[k]*A.A.n+j] = A.A.x[k];
1631 #if(CXSC_INDEX_CHECK) 1632 if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
1633 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator()(int, int, int, int)"));
1639 #if(CXSC_INDEX_CHECK) 1640 if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
1641 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator()(int, int, int, int) const"));
1663 return sp_m_negative<srmatrix,srmatrix>(M.A);
1679 return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1.A,M2.A);
1690 return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1.A,M2);
1701 return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1,M2.A);
1712 return spf_mm_mult<srmatrix,rmatrix,rmatrix,sparse_dot>(M1.A,M2);
1723 return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(M1,M2.A);
1734 return spf_mm_mult<srmatrix,rmatrix_slice,rmatrix,sparse_dot>(M1.A,M2);
1745 return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(M1,M2.A);
1756 return spsp_mv_mult<srmatrix,srvector,srvector,sparse_dot,real>(M.A,v);
1767 return spsl_mv_mult<srmatrix,srvector_slice,srvector,sparse_dot,real>(M.A,v);
1778 return spf_mv_mult<srmatrix,rvector,rvector,sparse_dot>(M.A,v);
1789 return spf_mv_mult<srmatrix,rvector_slice,rvector,sparse_dot>(M.A,v);
1794 return sp_ms_mult<srmatrix,real,srmatrix>(M.A,r);
1799 return sp_ms_div<srmatrix,real,srmatrix>(M.A,r);
1804 return sp_sm_mult<real,srmatrix,srmatrix>(r,M.A);
1809 return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1.A,M2.A);
1814 return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1.A,M2);
1819 return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1,M2.A);
1824 return spf_mm_add<srmatrix,rmatrix,rmatrix>(M1.A,M2);
1829 return fsp_mm_add<rmatrix,srmatrix,rmatrix>(M1,M2.A);
1834 return spf_mm_add<srmatrix,rmatrix_slice,rmatrix>(M1.A,M2);
1839 return fsp_mm_add<rmatrix_slice,srmatrix,rmatrix>(M1,M2.A);
1844 return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1.A,M2.A);
1849 return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1.A,M2);
1854 return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1,M2.A);
1859 return spf_mm_sub<srmatrix,rmatrix,rmatrix>(M1.A,M2);
1864 return fsp_mm_sub<rmatrix,srmatrix,rmatrix>(M1,M2.A);
1869 return spf_mm_sub<srmatrix,rmatrix_slice,rmatrix>(M1.A,M2);
1874 return fsp_mm_sub<rmatrix_slice,srmatrix,rmatrix>(M1,M2.A);
1914 return spsp_mm_comp(M1.A,M2.A);
1919 return spsp_mm_comp(M1.A,M2);
1924 return spsp_mm_comp(M1,M2.A);
1929 return spf_mm_comp(M1.A,M2);
1934 return fsp_mm_comp(M1,M2.A);
1939 return fsp_mm_comp(M1,M2.A);
1944 return spf_mm_comp(M1.A,M2);
1949 return !spsp_mm_comp(M1.A,M2.A);
1954 return !spsp_mm_comp(M1.A,M2);
1959 return !spsp_mm_comp(M1,M2.A);
1964 return !spf_mm_comp(M1.A,M2);
1969 return !fsp_mm_comp(M1,M2.A);
1974 return !fsp_mm_comp(M1,M2.A);
1979 return !spf_mm_comp(M1.A,M2);
1984 return spsp_mm_less<srmatrix,srmatrix,real>(M1.A,M2.A);
1989 return spsp_mm_less<srmatrix,srmatrix,real>(M1.A,M2);
1994 return spsp_mm_less<srmatrix,srmatrix,real>(M1,M2.A);
1999 return spf_mm_less<srmatrix,rmatrix,real>(M1.A,M2);
2004 return fsp_mm_less<rmatrix,srmatrix,real>(M1,M2.A);
2009 return fsp_mm_less<rmatrix_slice,srmatrix,real>(M1,M2.A);
2014 return spf_mm_less<srmatrix,rmatrix_slice,real>(M1.A,M2);
2019 return spsp_mm_leq<srmatrix,srmatrix,real>(M1.A,M2.A);
2024 return spsp_mm_leq<srmatrix,srmatrix,real>(M1.A,M2);
2029 return spsp_mm_leq<srmatrix,srmatrix,real>(M1,M2.A);
2034 return spf_mm_leq<srmatrix,rmatrix,real>(M1.A,M2);
2039 return fsp_mm_leq<rmatrix,srmatrix,real>(M1,M2.A);
2044 return fsp_mm_leq<rmatrix_slice,srmatrix,real>(M1,M2.A);
2049 return spf_mm_leq<srmatrix,rmatrix_slice,real>(M1.A,M2);
2054 return spsp_mm_greater<srmatrix,srmatrix,real>(M1.A,M2.A);
2059 return spsp_mm_greater<srmatrix,srmatrix,real>(M1.A,M2);
2064 return spsp_mm_greater<srmatrix,srmatrix,real>(M1,M2.A);
2069 return spf_mm_greater<srmatrix,rmatrix,real>(M1.A,M2);
2074 return fsp_mm_greater<rmatrix,srmatrix,real>(M1,M2.A);
2079 return fsp_mm_greater<rmatrix_slice,srmatrix,real>(M1,M2.A);
2084 return spf_mm_greater<srmatrix,rmatrix_slice,real>(M1.A,M2);
2089 return spsp_mm_geq<srmatrix,srmatrix,real>(M1.A,M2.A);
2094 return spsp_mm_geq<srmatrix,srmatrix,real>(M1.A,M2);
2099 return spsp_mm_geq<srmatrix,srmatrix,real>(M1,M2.A);
2104 return spf_mm_geq<srmatrix,rmatrix,real>(M1.A,M2);
2109 return fsp_mm_geq<rmatrix,srmatrix,real>(M1,M2.A);
2114 return fsp_mm_geq<rmatrix_slice,srmatrix,real>(M1,M2.A);
2119 return spf_mm_geq<srmatrix,rmatrix_slice,real>(M1.A,M2);
2124 return sp_m_not(M.A);
2134 return sp_m_output<srmatrix,real>(os, M.A);
2145 sp_m_input<srmatrix,real>(is, tmp);
2164 srmatrix_subv(
srmatrix& A,
bool r,
int i,
int j,
int k,
int l) : dat(A,i,j,k,l), row(r) {
2165 if(row) index=i;
else index=k;
2168 srmatrix_subv(
const srmatrix& A,
bool r,
int i,
int j,
int k,
int l) : dat(A,i,j,k,l), row(r){
2169 if(row) index=i;
else index=k;
2181 #if(CXSC_INDEX_CHECK) 2182 if(i<dat.A.lb2 || i>dat.A.ub2)
2183 cxscthrow(ELEMENT_NOT_IN_VEC(
"srmatrix_subv::operator[](int)"));
2187 #if(CXSC_INDEX_CHECK) 2188 if(i<dat.A.lb1 || i>dat.A.ub1)
2189 cxscthrow(ELEMENT_NOT_IN_VEC(
"srmatrix_subv::operator[](int)"));
2201 #if(CXSC_INDEX_CHECK) 2202 if(i<dat.A.lb2 || i>dat.A.ub2)
2203 cxscthrow(ELEMENT_NOT_IN_VEC(
"srmatrix_subv::operator[](int)"));
2205 return dat(index,i);
2207 #if(CXSC_INDEX_CHECK) 2208 if(i<dat.A.lb1 || i>dat.A.ub1)
2209 cxscthrow(ELEMENT_NOT_IN_VEC(
"srmatrix_subv::operator[](int)"));
2211 return dat(i,index);
2217 return svsp_vv_assign(*
this,
srvector(v));
2222 return sv_vs_assign(*
this,v);
2227 return svsp_vv_assign(*
this,v);
2232 return svsl_vv_assign(*
this,v);
2237 return svf_vv_assign(*
this,v);
2242 return svf_vv_assign(*
this,v);
2287 #include "vector_friend_declarations.inl" 2293 return Lb(S.dat, 2);
2295 return Lb(S.dat, 1);
2301 return Ub(S.dat, 2);
2303 return Ub(S.dat, 1);
2308 return Ub(S)-
Lb(S)+1;
2320 if(v.row) n=v.dat.A.n;
else n=v.dat.A.m;
2328 #if(CXSC_INDEX_CHECK) 2329 if(c.col()<lb2 || c.col()>ub2)
2330 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator[](const cxscmatrix_column&)"));
2332 return srmatrix_subv(*
this,
false, lb1, ub1, c.col(), c.col());
2336 #if(CXSC_INDEX_CHECK) 2338 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator[](const int)"));
2344 #if(CXSC_INDEX_CHECK) 2345 if(c.col()<lb2 || c.col()>ub2)
2346 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator[](const cxscmatrix_column&)"));
2348 return srmatrix_subv(*
this,
false, lb1, ub1, c.col(), c.col());
2352 #if(CXSC_INDEX_CHECK) 2354 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix::operator[](const int)"));
2360 #if(CXSC_INDEX_CHECK) 2361 if(i<A.lb1 || i>A.ub1)
2362 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix_slice::operator[](const int)"));
2368 #if(CXSC_INDEX_CHECK) 2369 if(c.col()<A.lb2 || c.col()>A.ub2)
2370 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix_slice::operator[](const cxscmatrix_column&)"));
2372 return srmatrix_subv(*M,
false, A.lb1, A.ub1, c.col(), c.col());
2376 #if(CXSC_INDEX_CHECK) 2377 if(i<A.lb1 || i>A.ub1)
2378 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix_slice::operator[](const int)"));
2384 #if(CXSC_INDEX_CHECK) 2385 if(c.col()<A.lb2 || c.col()>A.ub2)
2386 cxscthrow(ROW_OR_COL_NOT_IN_MAT(
"srmatrix_slice::operator[](const cxscmatrix_column&)"));
2388 return srmatrix_subv(*M,
false, A.lb1, A.ub1, c.col(), c.col());
2400 for(
int j=0 ; j<n ; j++) {
2401 for(
int k=A.dat.A.p[j] ; k<A.dat.A.p[j+1] ; k++) {
2403 x.push_back(A.dat.A.x[k]);
2412 for(
unsigned int k=0 ; k<A.dat.A.ind.size() ; k++) {
2413 p.push_back(A.dat.A.ind[k]);
2414 x.push_back(A.dat.A.x[k]);
2943 spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot,
srvector(v1),
srvector(v2));
2951 spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot,
srvector(v1), v2);
2959 spsl_vv_accu<dotprecision,srvector,srvector_slice,sparse_dot>(dot,
srvector(v1), v2);
2967 spf_vv_accu<dotprecision,srvector,rvector,sparse_dot>(dot,
srvector(v1), v2);
2975 spf_vv_accu<dotprecision,srvector,rvector_slice,sparse_dot>(dot,
srvector(v1), v2);
2983 spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, v1,
srvector(v2));
2991 slsp_vv_accu<dotprecision,srvector_slice,srvector,sparse_dot>(dot, v1,
srvector(v2));
2999 fsp_vv_accu<dotprecision,rvector,srvector,sparse_dot>(dot, v1,
srvector(v2));
3007 fsp_vv_accu<dotprecision,rvector_slice,srvector,sparse_dot>(dot, v1,
srvector(v2));
3017 spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot,
srvector(v1),
srvector(v2));
3027 spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot,
srvector(v1), v2);
3037 spsl_vv_accuapprox<dotprecision,srvector,srvector_slice,sparse_dot>(dot,
srvector(v1), v2);
3047 spf_vv_accuapprox<dotprecision,srvector,rvector,sparse_dot>(dot,
srvector(v1), v2);
3057 spf_vv_accuapprox<dotprecision,srvector,rvector_slice,sparse_dot>(dot,
srvector(v1), v2);
3067 spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, v1,
srvector(v2));
3077 slsp_vv_accuapprox<dotprecision,srvector_slice,srvector,sparse_dot>(dot, v1,
srvector(v2));
3087 fsp_vv_accuapprox<dotprecision,rvector,srvector,sparse_dot>(dot, v1,
srvector(v2));
3097 fsp_vv_accuapprox<dotprecision,rvector_slice,srvector,sparse_dot>(dot, v1,
srvector(v2));
3503 #include "sparsematrix.inl" rmatrix_subv & operator=(const rmatrix_subv &rv)
Implementation of standard assigning operator.
const real operator[](const int i) const
Returns a copy of the i-th element of the subvector.
friend srmatrix Re(const scmatrix &)
Returns the real part of the sparse matrix A.
rmatrix & operator-=(const srmatrix &m)
Implementation of substraction and allocation operation.
void set_k(unsigned int i)
Set precision for computation of dot products.
srmatrix & operator*=(const rmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
srmatrix & operator=(const rmatrix_slice &A)
Assigns a dense matrix slice to the sparse matrix. Only the non zero entries of the dense matrix slic...
The Data Type rmatrix_slice.
srmatrix_slice & operator=(const rmatrix &C)
Assing C to the slice.
srmatrix_subv & operator*=(const real &)
Assign the componentwise product of the subvector with a scalar to the subvector. ...
srmatrix_slice & operator-=(const srmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
friend void SetLb(srmatrix &, const int, const int)
Sets the lower index bound of the rows (i==ROW) or columns (i==COL) to j.
srmatrix(const int r, const int c)
Creates an empty matrix with r rows and c columns, pre-reserving space for 2*(r+c) elements...
srmatrix_slice & operator-=(const rmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
The Data Type idotprecision.
friend void SetUb(srmatrix &, const int, const int)
Sets the upper index bound of the rows (i==ROW) or columns (i==COL) to j.
The Data Type dotprecision.
srvector()
Default constructor, creates an empty vector of size 0.
friend srmatrix Im(const scmatrix &)
Returns the imaginary part of the sparse matrix A.
srmatrix_slice & operator+=(const srmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
real & element(const int i, const int j)
Returns a reference to the element (i,j) of the matrix.
friend int ColLen(const srmatrix_slice &)
Returns the number of rows of the matrix slice.
srmatrix & operator*=(const rmatrix_slice &B)
Multiply the sparse matrix by B and assign the result to it.
Represents a row or column vector of a sparse matrix.
friend srmatrix Sup(const simatrix &)
Returns the Supremum of the matrix A.
friend srmatrix absmin(const simatrix &)
Returns the componentwise minimum absolute value.
srmatrix_subv & operator=(const real &v)
Assigns v to all elements of the subvector.
srmatrix & operator+=(const rmatrix &B)
Add B to the sparse matrix and assign the result to it.
friend int Ub(const srmatrix_slice &, const int)
Returns the upper index bound of the rows (if i==ROW) or columns (if i==COL) of the slice...
int Lb(const cimatrix &rm, const int &i)
Returns the lower bound index.
srmatrix_subv & operator=(const srmatrix_subv &v)
Assigns a vector to a subvector.
The namespace cxsc, providing all functionality of the class library C-XSC.
civector operator*(const cimatrix_subv &rv, const cinterval &s)
Implementation of multiplication operation.
srmatrix(const int m, const int n, const int nnz, const intvector &rows, const intvector &cols, const rvector &values, const enum STORAGE_TYPE t=triplet)
Creates a sparse matrix out of three vectors (arrays) forming a matrix stored in triplet, compressed row or compressed column storage.
const real operator()(int i, int j) const
Returns a copy of the element in row i and column j.
std::vector< int > & row_indices()
Returns a reference to the vector containing the row indices (the array)
srmatrix(const int r, const int c, const int e)
Creates an empty matrix with r rows and c columns, pre-reserving space for e elements.
A sparse interval vector.
srmatrix & operator/=(const real &r)
Divide all elements of the sparse matrix by r and assign the result to it.
real & operator[](const int i)
Returns a reference to the i-th element of the subvector.
cimatrix & SetLb(cimatrix &m, const int &i, const int &j)
Sets the lower bound index.
srmatrix_slice & operator*=(const srmatrix &M)
Assigns the product of the sparse slice and M to the slice.
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
srmatrix operator()(const intvector &pervec, const intvector &q)
Performs a row and column permutation using two permutation vectors.
real & element(int i, int j)
Returns a reference to the element (i,j) of the matrix.
int get_k() const
Get currently set precision for computation of dot products.
srmatrix_subv operator[](const cxscmatrix_column &)
Returns a column of the matrix as a sparse subvector object.
srmatrix_slice & operator+=(const rmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
A sparse complex interval matrix.
void accumulate_approx(cdotprecision &dp, const cmatrix_subv &rv1, const cmatrix_subv &rv2)
The accurate scalar product of the last two arguments added to the value of the first argument (witho...
std::vector< int > & column_pointers()
Returns a reference to the vector containing the column pointers (the array)
srmatrix()
Standard constructor, creates an empty matrix of dimension 0x0.
srmatrix_subv & operator=(const rvector &v)
Assigns a vector to a subvector.
srmatrix_slice & operator-=(const rmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
friend int Ub(const srmatrix &, int)
Returns the upper index bound for the rows or columns of A.
friend srmatrix SupIm(const scimatrix &)
Returns the imaginary part of the supremum of the matrix A.
friend int VecLen(const srmatrix_subv &)
Returns the length of the subvector.
rmatrix & operator+=(const srmatrix &m)
Implementation of addition and allocation operation.
friend int RowLen(const srmatrix &)
Returns the number of columns of the matrix.
srmatrix(const rmatrix &A)
Creates a sparse matrix out of a dense matrix A. Only the non zero elements of A are stored explicitl...
rmatrix_slice & operator*=(const srmatrix &m)
Implementation of multiplication and allocation operation.
A sparse complex interval vector.
const std::vector< int > & row_indices() const
Returns a constant reference to the vector containing the row indices (the array) ...
srmatrix & operator*=(const real &r)
Multiply all elements of the sparse matrix by r and assign the result to it.
srmatrix & operator-=(const rmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
srmatrix_slice & operator=(const srmatrix_slice &C)
Assing C to the slice.
srmatrix_slice & operator+=(const srmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
The Data Type cidotprecision.
cimatrix Id(const cimatrix &A)
Returns the Identity matrix.
friend srmatrix mid(const simatrix &)
Returns the midpoint matrix for A.
The Data Type rvector_slice.
friend srmatrix transp(const srmatrix &)
Returns the transpose of A.
A slice of a sparse real interval matrix.
srmatrix(const int m, const int n, const int nnz, const int *rows, const int *cols, const real *values, const enum STORAGE_TYPE t=triplet)
Creates a sparse matrix out of three arrays forming a matrix stored in triplet, compressed row or com...
friend srmatrix Id(const srmatrix &)
Return a sparse unity matrix of the same dimension as A.
cimatrix & SetUb(cimatrix &m, const int &i, const int &j)
Sets the upper bound index.
srmatrix_slice & operator*=(const real &r)
Assigns the component wise product of the sparse slice and r to the slice.
The Data Type rmatrix_subv.
srmatrix & operator-=(const srmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
srmatrix_subv & operator=(const srvector &v)
Assigns a vector to a subvector.
std::vector< real > & values()
Returns a reference to the vector containing the stored values (the array)
srmatrix operator()(const intmatrix &P)
Performs a row permutation using the permutation matrix P. Faster than explicitly computing the produ...
friend int Lb(const srmatrix &, int)
Returns the lower index bound for the rows or columns of A.
int get_nnz() const
Returns the number of non-zero entries (including explicitly stored zeros).
A slice of a sparse real matrix.
rmatrix_slice & operator+=(const srmatrix &m)
Implementation of multiplication and allocation operation.
The Data Type cdotprecision.
srmatrix_slice & operator=(const rmatrix_slice &C)
Assing C to the slice.
A slice of a sparse complex matrix.
void Resize(cimatrix &A)
Resizes the matrix.
srmatrix & operator=(const real &A)
Assigns a real value to all elements of the matrix (resulting in a dense matrix!) ...
srmatrix_subv & operator+=(const srvector &)
Assign the sum of the subvector with a vector to the subvector.
int get_k() const
Get currently set precision for computation of dot products.
const std::vector< int > & column_pointers() const
Returns a constant reference to the vector containing the column pointers (the array) ...
rmatrix_subv & operator-=(const real &c)
Implementation of subtraction and allocation operation.
srmatrix_subv & operator=(const rvector_slice &v)
Assigns a vector to a subvector.
rmatrix_slice & operator-=(const srmatrix &m)
Implementation of multiplication and allocation operation.
friend srmatrix abs(const srmatrix &)
Returns the componentwise absolute value of A.
srmatrix & operator-=(const rmatrix_slice &B)
Subtract B from the sparse matrix and assign the result to it.
srmatrix & operator+=(const rmatrix_slice &B)
Add B to the sparse matrix and assign the result to it.
Represents a row or column vector of a sparse matrix.
friend srmatrix absmax(const simatrix &)
Returns the componentwise maximum absolute value.
int ColLen(const cimatrix &)
Returns the column dimension.
srmatrix & operator*=(const srmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
friend srmatrix InfRe(const scimatrix &)
Returns the real part of the infimum of the matrix A.
friend srmatrix CompMat(const simatrix &)
Returns Ostroswkis comparison matrix for A.
rmatrix_subv & operator+=(const real &c)
Implementation of addition and allocation operation.
srmatrix_subv operator[](const int)
Returns a row of the matrix.
rmatrix_slice & operator=(const srmatrix &m)
Implementation of standard assigning operator.
srmatrix_slice & operator=(const real &C)
Assing C to all elements of the slice.
void full(rmatrix &A) const
Creates a full matrix out of the sparse matrix and stores it in A. This should normally be done using...
int get_k() const
Get currently set precision for computation of dot products.
srmatrix_slice & operator*=(const rmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
int RowLen(const cimatrix &)
Returns the row dimension.
srmatrix_slice & operator-=(const srmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
rmatrix & operator*=(const srmatrix &m)
Implementation of multiplication and allocation operation.
friend srmatrix SupRe(const scimatrix &)
Returns the real part of the supremum of the matrix A.
rmatrix CompMat(const cimatrix &A)
Returns Ostrowski's comparison matrix.
srmatrix_subv & operator-=(const srvector &)
Assign the difference of the subvector with a vector to the subvector.
srmatrix_subv & operator=(const srvector_slice &v)
Assigns a vector to a subvector.
friend std::istream & operator>>(std::istream &, srmatrix_slice &)
Standard input operator for sparse matrix slice.
srmatrix_slice & operator*=(const srmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
srmatrix_subv & operator/=(const real &)
Assign the componentwise division of the subvector with a scalar to the subvector.
friend srmatrix diam(const simatrix &)
Returns the componentwise diameter of A.
srmatrix_slice & operator*=(const rmatrix &M)
Assigns the product of the sparse slice and M to the slice.
friend int RowLen(const srmatrix_slice &)
Returns the number columns of the matrix slice.
void dropzeros()
Drops explicitly stored zeros from the data structure.
srmatrix_slice & operator+=(const rmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Represents a row or column vector of a sparse matrix.
friend srmatrix InfIm(const scimatrix &)
Returns the imaginary part of the infimum of the matrix A.
const real operator()(const int i, const int j) const
Returns a copy of the element (i,j) of the matrix.
STORAGE_TYPE
Enumeration depicting the storage type of a sparse matrix (Triplet storage, Compressed column storage...
friend srmatrix Inf(const simatrix &)
Returns the Infimum of the matrix A.
Helper class for slices of sparse vectors.
real density() const
Returns the density (the number of non-zeros divided by the number of elements) of the matrix...
int Ub(const cimatrix &rm, const int &i)
Returns the upper bound index.
friend int ColLen(const srmatrix &)
Returns the number of rows of the matrix.
srmatrix_slice & operator/=(const real &r)
Assigns the component wise division of the sparse slice and M to the slice.
Represents a row or column vector of a sparse matrix.
rmatrix()
Constructor of class rmatrix.
rmatrix & operator=(const real &r)
Implementation of standard assigning operator.
civector operator/(const cimatrix_subv &rv, const cinterval &s)
Implementation of division operation.
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
A slice of a sparse complex interval matrix.
srmatrix_slice & operator=(const srmatrix &C)
Assing C to the slice.
srmatrix & operator=(const rmatrix &A)
Assigns a dense matrix to the sparse matrix. Only the non zero entries of the dense matrix are used...
srmatrix operator()(const intmatrix &P, const intmatrix &Q)
Performs row and column permutations using the two permutation matrices P and Q. Faster than explicit...
friend int Ub(const srmatrix_subv &)
Returns the upper index bound of the subvector.
srmatrix & operator+=(const srmatrix &B)
Add B to the sparse matrix and assign the result to it.
friend int Lb(const srmatrix_slice &, const int)
Returns the lower index bound of the rows (if i==ROW) or columns (if i==COL) of the slice...
srmatrix operator()(const intvector &pervec)
Performs a row permutation using a permutation vector.
friend srvector operator-(const srmatrix_subv &)
Unary negation operator.
A sparse interval matrix.
friend std::istream & operator>>(std::istream &, srmatrix_subv &)
Standard input operator for subvectors.
const std::vector< real > & values() const
Returns a constant reference to the vector containing the stored values (the array) ...
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
srmatrix(const int ms, const int ns, const rmatrix &A)
Constructor for banded matrices.
friend int Lb(const srmatrix_subv &)
Returns the lower index bound of the subvector.