26#ifndef _CXSC_SRMATRIX_HPP_INCLUDED
27#define _CXSC_SRMATRIX_HPP_INCLUDED
31#include <srvector.hpp>
36#include <sparsedot.hpp>
37#include <sparsematrix.hpp>
59inline 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);
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];
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]));
924inline 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);
1341inline std::ostream& operator<<(std::ostream& os,
const srmatrix& A) {
1342 return sp_m_output<srmatrix,real>(os,A);
1351inline 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"
The Data Type cdotprecision.
int get_k() const
Get currently set precision for computation of dot products.
The Data Type cidotprecision.
int get_k() const
Get currently set precision for computation of dot products.
The Data Type dotprecision.
void set_k(unsigned int i)
Set precision for computation of dot products.
The Data Type idotprecision.
int get_k() const
Get currently set precision for computation of dot products.
The Data Type rmatrix_slice.
rmatrix_slice & operator=(const srmatrix &m)
Implementation of standard assigning operator.
rmatrix_slice & operator-=(const srmatrix &m)
Implementation of multiplication and allocation operation.
rmatrix_slice & operator*=(const srmatrix &m)
Implementation of multiplication and allocation operation.
rmatrix_slice & operator+=(const srmatrix &m)
Implementation of multiplication and allocation operation.
The Data Type rmatrix_subv.
rmatrix_subv & operator+=(const real &c) noexcept
Implementation of addition and allocation operation.
rmatrix_subv & operator=(const rmatrix_subv &rv) noexcept
Implementation of standard assigning operator.
rmatrix_subv & operator-=(const real &c) noexcept
Implementation of subtraction and allocation operation.
rmatrix & operator=(const real &r) noexcept
Implementation of standard assigning operator.
rmatrix & operator-=(const srmatrix &m)
Implementation of substraction and allocation operation.
rmatrix() noexcept
Constructor of class rmatrix.
rmatrix & operator*=(const srmatrix &m)
Implementation of multiplication and allocation operation.
rmatrix & operator+=(const srmatrix &m)
Implementation of addition and allocation operation.
The Data Type rvector_slice.
A slice of a sparse complex interval matrix.
Represents a row or column vector of a sparse matrix.
A sparse complex interval matrix.
A sparse complex interval vector.
A slice of a sparse complex matrix.
Represents a row or column vector of a sparse matrix.
A slice of a sparse real interval matrix.
Represents a row or column vector of a sparse matrix.
A sparse interval matrix.
A sparse interval vector.
A slice of a sparse real matrix.
srmatrix_slice & operator-=(const rmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
srmatrix_slice & operator+=(const rmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
srmatrix_subv operator[](const int)
Returns a row of the matrix.
real & element(const int i, const int j)
Returns a reference to the element (i,j) of the matrix.
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.
srmatrix_slice & operator+=(const rmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
srmatrix_slice & operator*=(const real &r)
Assigns the component wise product of the sparse slice and r to the slice.
srmatrix_slice & operator*=(const srmatrix &M)
Assigns the product of the sparse slice and M to the slice.
srmatrix_slice & operator-=(const rmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
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_slice & operator+=(const srmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
srmatrix_slice & operator=(const srmatrix &C)
Assing C to the slice.
srmatrix_slice & operator=(const real &C)
Assing C to all elements of the slice.
srmatrix_slice & operator*=(const rmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
srmatrix_slice & operator*=(const srmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
friend int ColLen(const srmatrix_slice &)
Returns the number of rows of the matrix slice.
srmatrix_slice & operator*=(const rmatrix &M)
Assigns the product of the sparse slice and M to the slice.
srmatrix_slice & operator=(const srmatrix_slice &C)
Assing C to the slice.
friend int RowLen(const srmatrix_slice &)
Returns the number columns of the matrix slice.
const real operator()(const int i, const int j) const
Returns a copy of the element (i,j) of the matrix.
srmatrix_slice & operator/=(const real &r)
Assigns the component wise division of the sparse slice and M to the slice.
srmatrix_slice & operator=(const rmatrix &C)
Assing C to the slice.
srmatrix_slice & operator=(const rmatrix_slice &C)
Assing C to the slice.
srmatrix_slice & operator-=(const srmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
srmatrix_slice & operator+=(const srmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
srmatrix_slice & operator-=(const srmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Represents a row or column vector of a sparse matrix.
srmatrix_subv & operator/=(const real &)
Assign the componentwise division of the subvector with a scalar to the subvector.
friend std::istream & operator>>(std::istream &, srmatrix_subv &)
Standard input operator for subvectors.
friend int VecLen(const srmatrix_subv &)
Returns the length of the subvector.
srmatrix_subv & operator=(const srvector &v)
Assigns a vector to a subvector.
srmatrix_subv & operator+=(const srvector &)
Assign the sum of the subvector with a vector to the subvector.
friend int Ub(const srmatrix_subv &)
Returns the upper index bound of the subvector.
srmatrix_subv & operator=(const srmatrix_subv &v)
Assigns a vector to a subvector.
friend int Lb(const srmatrix_subv &)
Returns the lower index bound of the subvector.
friend srvector operator-(const srmatrix_subv &)
Unary negation operator.
srmatrix_subv & operator-=(const srvector &)
Assign the difference of the subvector with a vector to the subvector.
real & operator[](const int i)
Returns a reference to the i-th element of the subvector.
srmatrix_subv & operator=(const real &v)
Assigns v to all elements of the subvector.
srmatrix_subv & operator=(const rvector_slice &v)
Assigns a vector to a subvector.
srmatrix_subv & operator=(const srvector_slice &v)
Assigns a vector to a subvector.
srmatrix_subv & operator=(const rvector &v)
Assigns a vector to a subvector.
const real operator[](const int i) const
Returns a copy of the i-th element of the subvector.
srmatrix_subv & operator*=(const real &)
Assign the componentwise product of the subvector with a scalar to the subvector.
srmatrix operator()(const intmatrix &P)
Performs a row permutation using the permutation matrix P. Faster than explicitly computing the produ...
srmatrix(const int ms, const int ns, const rmatrix &A)
Constructor for banded matrices.
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...
srmatrix & operator+=(const rmatrix &B)
Add B to the sparse matrix and assign the result to it.
friend srmatrix absmin(const simatrix &)
Returns the componentwise minimum absolute value.
real density() const
Returns the density (the number of non-zeros divided by the number of elements) of the matrix.
friend srmatrix Inf(const simatrix &)
Returns the Infimum of the matrix A.
std::vector< int > & column_pointers()
Returns a reference to the vector containing the column pointers (the array)
srmatrix & operator*=(const rmatrix_slice &B)
Multiply the sparse matrix by B and assign the result to it.
srmatrix operator()(const intvector &pervec)
Performs a row permutation using a permutation vector.
friend void SetUb(srmatrix &, const int, const int)
Sets the upper index bound of the rows (i==ROW) or columns (i==COL) to j.
srmatrix & operator-=(const srmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
friend srmatrix Re(const scmatrix &)
Returns the real part of the sparse matrix A.
void dropzeros()
Drops explicitly stored zeros from the data structure.
srmatrix & operator*=(const rmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
srmatrix & operator*=(const srmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
friend srmatrix absmax(const simatrix &)
Returns the componentwise maximum absolute value.
srmatrix & operator/=(const real &r)
Divide all elements of the sparse matrix by r and assign the result to it.
friend srmatrix Id(const srmatrix &)
Return a sparse unity matrix of the same dimension as A.
friend srmatrix InfIm(const scimatrix &)
Returns the imaginary part of the infimum of the matrix A.
srmatrix operator()(const intmatrix &P, const intmatrix &Q)
Performs row and column permutations using the two permutation matrices P and Q. Faster than explicit...
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...
friend srmatrix diam(const simatrix &)
Returns the componentwise diameter of A.
friend srmatrix SupRe(const scimatrix &)
Returns the real part of the supremum of the matrix A.
srmatrix & operator+=(const srmatrix &B)
Add B to the sparse matrix and assign the result to it.
srmatrix operator()(const intvector &pervec, const intvector &q)
Performs a row and column permutation using two permutation vectors.
srmatrix & operator-=(const rmatrix_slice &B)
Subtract B from the sparse matrix and assign the result to it.
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 & operator=(const rmatrix &A)
Assigns a dense matrix to the sparse matrix. Only the non zero entries of the dense matrix are used.
friend srmatrix CompMat(const simatrix &)
Returns Ostroswkis comparison matrix for A.
srmatrix(const rmatrix &A)
Creates a sparse matrix out of a dense matrix A. Only the non zero elements of A are stored explicitl...
real & element(int i, int j)
Returns a reference to the element (i,j) of the matrix.
friend srmatrix transp(const srmatrix &)
Returns the transpose of A.
const real operator()(int i, int j) const
Returns a copy of the element in row i and column j.
friend srmatrix InfRe(const scimatrix &)
Returns the real part of the infimum of the matrix A.
friend int Ub(const srmatrix &, int)
Returns the upper index bound for the rows or columns of A.
srmatrix_subv operator[](const cxscmatrix_column &)
Returns a column of the matrix as a sparse subvector object.
srmatrix & operator+=(const rmatrix_slice &B)
Add B to the sparse matrix and assign the result to it.
friend srmatrix mid(const simatrix &)
Returns the midpoint matrix for A.
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.
srmatrix()
Standard constructor, creates an empty matrix of dimension 0x0.
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...
srmatrix & operator=(const real &A)
Assigns a real value to all elements of the matrix (resulting in a dense matrix!)
std::vector< int > & row_indices()
Returns a reference to the vector containing the row indices (the array)
const std::vector< real > & values() const
Returns a constant reference to the vector containing the stored values (the array)
srmatrix & operator*=(const real &r)
Multiply all elements of the sparse matrix by r and assign the result to it.
const std::vector< int > & column_pointers() const
Returns a constant reference to the vector containing the column pointers (the array)
int get_nnz() const
Returns the number of non-zero entries (including explicitly stored zeros).
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.
std::vector< real > & values()
Returns a reference to the vector containing the stored values (the array)
friend srmatrix Sup(const simatrix &)
Returns the Supremum of the matrix A.
srmatrix & operator-=(const rmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
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,...
friend srmatrix Im(const scmatrix &)
Returns the imaginary part of the sparse matrix A.
friend int Lb(const srmatrix &, int)
Returns the lower index bound for the rows or columns of A.
friend int RowLen(const srmatrix &)
Returns the number of columns of the matrix.
friend srmatrix SupIm(const scimatrix &)
Returns the imaginary part of the supremum of the matrix A.
const std::vector< int > & row_indices() const
Returns a constant reference to the vector containing the row indices (the array)
friend srmatrix abs(const srmatrix &)
Returns the componentwise absolute value of A.
friend int ColLen(const srmatrix &)
Returns the number of rows of the matrix.
friend std::istream & operator>>(std::istream &, srmatrix_slice &)
Standard input operator for sparse matrix slice.
Helper class for slices of sparse vectors.
srvector()
Default constructor, creates an empty vector of size 0.
The namespace cxsc, providing all functionality of the class library C-XSC.
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
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...
int ColLen(const cimatrix &)
Returns the column dimension.
rmatrix CompMat(const cimatrix &A)
Returns Ostrowski's comparison matrix.
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
int Ub(const cimatrix &rm, const int &i) noexcept
Returns the upper bound index.
cimatrix & SetLb(cimatrix &m, const int &i, const int &j) noexcept
Sets the lower bound index.
int RowLen(const cimatrix &)
Returns the row dimension.
STORAGE_TYPE
Enumeration depicting the storage type of a sparse matrix (Triplet storage, Compressed column storage...
cimatrix Id(const cimatrix &A)
Returns the Identity matrix.
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
void Resize(cimatrix &A) noexcept
Resizes the matrix.
cimatrix & SetUb(cimatrix &m, const int &i, const int &j) noexcept
Sets the upper bound index.
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.