13 #ifndef TBCI_CSCMATRIX_H
14 #define TBCI_CSCMATRIX_H
16 #include "tbci/vector.h"
17 #include "tbci/matrix_sig.h"
18 #include "tbci/f_matrix.h"
19 #include "tbci/band_matrix.h"
22 #if !defined(NO_GD) && !defined(AUTO_DECL)
23 # include "tbci/cscmatrix_gd.h"
29 #ifndef TBCI_DISABLE_EXCEPT
36 :
NumErr(
"Error in CSCMatrix") {}
46 # pragma interface "cscmatrix.h"
52 template <
typename T>
class F_Matrix;
77 const unsigned int columns,
const unsigned int nnzeros=1)
79 if (val != (
T)0)
fill(val); }
81 const unsigned int nnzeros=1)
82 {
allocate(rows, columns, nnzeros); }
97 static const char*
mat_info () {
return (
"CSCMatrix"); }
108 const unsigned int newColumns,
109 const unsigned int nnzeros = 1);
119 typename tbci_traits<T>::const_refval_type
121 typename tbci_traits<T>::const_refval_type
122 get (
unsigned int row,
unsigned int column)
const HOT;
126 T&
setval(
unsigned int row,
unsigned int column);
128 {
return setval (row, column); }
141 template <
typename MatType>
151 template <
typename MatType>
157 for (
unsigned int j=
pcol[
i]; j<
pcol[
i+1]; ++j)
206 friend STD__ ostream& operator<< FGDT (STD__ ostream& stream, const CSCMatrix<T>& m);
224 void insert (
const unsigned int column,
const unsigned int pos);
228 template <
typename MatType>
229 void fill (
const MatType& M,
const double tol = 0.0)
231 unsigned int count = 0,
i, j;
233 for (j = 0; j <
n_cols; ++j)
237 STD__ cout <<
"Matrix conversion, NNZ=" << count <<
STD__ endl;
238 resize(M.rows(), M.columns(), count);
239 for (j = 0; j <
n_cols; ++j)
247 template <
typename T>
248 STD__ ostream& operator<< (STD__ ostream& stream, const CSCMatrix<T>& m)
250 for (
unsigned int i = 0;
i < m.rows(); ++
i) {
252 for (
unsigned int j = 0; j < m.columns(); ++j)
253 stream << m(
i,j) <<
" " ;
264 template <
typename T>
270 for (
unsigned int i = 0;
i < m1.
rows(); ++
i)
271 for (
unsigned int j = 0; j <
columns(); ++j) {
274 for (
unsigned int x = pcol[j];
x < pcol[j+1]; ++
x) {
280 res.setval(val-corr,
i, j);
283 for (
unsigned int i = 0;
i < m1.
rows(); ++
i)
284 for (
unsigned int j = 0; j <
columns(); ++j) {
286 for (
unsigned int x = pcol[j];
x < pcol[j+1]; ++
x)
287 val += m1(
i, irow[
x] )*comp[
x];
288 res.setval(val,
i, j);
295 template <typename T>
297 {
return m2.
multf (m1); }
300 template <
typename T>
306 for (
unsigned int i = 0;
i <
rows(); ++
i)
307 for (
unsigned int j = 0; j < m1.
columns(); ++j) {
310 for (
unsigned int x = 0;
x <
columns(); ++
x) {
316 res.setval(val-corr,
i, j);
319 for (
unsigned int i = 0;
i <
rows(); ++
i)
320 for (
unsigned int j = 0; j < m1.
columns(); ++j) {
323 val += (*
this)(
i,
x) * m1(
x, j);
324 res.setval(val,
i, j);
331 template <
typename T>
337 for (
unsigned int off = pcol[
c]; off < pcol[
c+1]; ++off) {
339 const unsigned int r = irow[off];
340 const T el = comp[off];
341 for (
unsigned int j = 0; j < m1.
columns(); ++j)
342 res.setval(r, j) += el * m1(
c, j);
348 template <typename T>
350 {
return m1.
mult1 (m2); }
353 template <
typename T>
359 for (
unsigned int r0 = 0; r0 <
rows(); ++r0)
360 for (
unsigned int c1 = 0; c1 < m1.
columns(); ++c1) {
363 for (
unsigned int off1 = m1.
pcol[c1]; off1 < m1.
pcol[c1+1]; ++off1) {
369 res.setval(val-corr, r0, c1);
372 for (
unsigned int r0 = 0; r0 <
rows(); ++r0)
373 for (
unsigned int c1 = 0; c1 < m1.
columns(); ++c1) {
375 for (
unsigned int off1 = m1.
pcol[c1]; off1 < m1.
pcol[c1+1]; ++off1)
376 val += (*this)(r0, m1.
irow[off1]) * m1.
comp[off1];
377 res.setval(val, r0, c1);
384 template <
typename T>
390 for (
unsigned int c = 0;
c <
columns(); ++
c) {
391 unsigned int off0 = pcol[
c];
392 unsigned int off1 = m1.
pcol[
c];
393 while (off0 < pcol[
c+1] || off1 < m1.
pcol[
c+1]) {
395 if (off0 >= pcol[
c+1]) {
398 }
else if (off1 >= m1.
pcol[
c+1]) {
399 res(irow[off0],
c) = comp[off0];
401 }
else if (irow[off0] == m1.
irow[off1]) {
403 res(irow[off0],
c) = comp[off0] + m1.
comp[off1];
405 }
else if (irow[off0] < m1.
irow[off1]) {
406 res(irow[off0],
c) = comp[off0];
417 template <
typename T>
423 for (
unsigned int c = 0;
c <
columns(); ++
c) {
424 unsigned int off0 = pcol[
c];
425 unsigned int off1 = m1.
pcol[
c];
426 while (off0 < pcol[
c+1] || off1 < m1.
pcol[
c+1]) {
428 if (off0 >= pcol[
c+1]) {
432 if (off1 >= m1.
pcol[
c+1]) {
433 res(irow[off0],
c) = comp[off0];
437 if (irow[off0] == m1.
irow[off1]) {
438 res(irow[off0],
c) = comp[off0] - m1.
comp[off1];
440 }
else if (irow[off0] < m1.
irow[off1]) {
441 res(irow[off0],
c) = comp[off0];
453 template <
typename T>
454 inline typename tbci_traits<T>::const_refval_type
457 static T dummy = (
T)0;
500 const unsigned int spos = pcol[
col], epos = pcol[col+1];
503 if (irow[spos] > row || irow[epos-1] < row)
505 if (irow[spos] == row)
507 if (irow[epos-1] == row)
510 if ((epos - spos) >= 32) {
511 const unsigned long int off =
_bin_search(irow, row, spos, epos);
512 if (off != (
unsigned long)(-1))
518 if (row <= (irow[spos] + irow[epos-1] + 1) / 2) {
519 for (
REGISTER unsigned int pos = spos; pos < epos; ++pos) {
520 if (irow[pos] == row)
522 else if (irow[pos] > row)
526 for (
REGISTER unsigned int pos = epos; pos > spos; --pos) {
527 const unsigned po1 = pos-1;
528 if (irow[po1] == row)
530 else if (irow[po1] < row)
538 template <
typename T>
539 inline typename tbci_traits<T>::const_refval_type
548 template <
typename T>
553 for (
unsigned int off = pcol[
c]; off < pcol[
c+1]; ++off)
554 res.
setval(comp[off], irow[off],
c);
560 template <
typename T>
562 const unsigned int newColumns,
563 const unsigned int nnzeros)
565 if (
UNLIKELY(newRows != n_rows && newColumns != n_cols)) {
567 allocate(newRows, newColumns, nnzeros);
568 }
else if (
UNLIKELY(nnzeros != n_max_size)) {
570 unsigned int *newirow =
new unsigned int[nnzeros];
571 T *newcomp =
new T [nnzeros];
573 const unsigned int com_sz =
MIN(nnzeros, n_max_size);
574 for (i = 0; i < com_sz; ++
i) {
575 newirow[
i] = irow[
i];
576 newcomp[
i] = comp[
i];
578 for (; i < nnzeros; ++
i)
580 n_max_size = nnzeros;
590 template <
typename T>
593 if (
UNLIKELY(n_size >= n_max_size)) {
596 unsigned int* new_irow =
new unsigned int [n_max_size];
598 T* new_comp =
new T [n_max_size];
604 for (i = 0; i < pos; ++
i) {
606 new_irow[
i] = irow[
i];
607 new_comp[
i] = comp[
i];
609 for (i = pos; i < n_size; ++
i) {
611 new_irow[i+1] = irow[
i];
612 new_comp[i+1] = comp[
i];
615 for (i = n_size+1; i < n_max_size; ++
i)
625 }
else if (n_size > 0) {
627 for (
int j = n_size-1; j >= (int)pos; --j) {
633 for (
unsigned int i = column+1;
i < n_cols+1; ++
i)
639 template <
typename T>
642 setval(row, column) = z;
647 template <
typename T>
652 unsigned int pos = pcol[column];
653 const unsigned int posn = pcol[column+1];
654 while (pos < posn && irow[pos] < row)
656 if (pos < posn && irow[pos] == row)
671 template <
typename T>
682 template <
typename T>
687 for (
unsigned int i = 0;
i < n_rows; ++
i) {
688 for (
unsigned int j = 0; j < n_cols; ++j) {
689 if (
operator()(
i,j) != cm(
i,j))
698 template <
typename T>
702 for (
unsigned int i = 0;
i < n_size; ++
i)
709 template <
typename T>
715 for (
unsigned row = 0;
row <= n_rows; ++
row) {
718 for (
unsigned int c = 0;
c < n_cols; ++
c) {
727 res.setval(
row) = val - corr;
732 template <
typename T>
737 for (
unsigned int i = 0;
i < n_cols; ++
i) {
738 for (
unsigned int j = pcol[
i]; j < pcol[
i+1]; ++j) {
739 const unsigned row = irow[j];
742 corr.
setval(row) += (t-
res.get(row)) - y;
751 template <
typename T>
755 return this->cscm_vec_mul_exact(v);
759 for (
unsigned int i = 0;
i < n_cols; ++
i) {
760 for (
unsigned int j = pcol[
i]; j < pcol[
i+1]; ++j)
761 res.
setval(irow[j]) += comp[j]*v(
i);
766 template <
typename T>
770 return (this->
operator * (v));
774 template <
typename T>
780 for (
unsigned int i = 0;
i < n_cols; ++
i) {
781 for (
unsigned int j = pcol[
i]; j < pcol[
i+1]; ++j)
790 template <
typename T>
796 for (
unsigned int i = 0;
i < n_cols; ++
i) {
797 for (
unsigned int j = pcol[
i]; j < pcol[
i+1]; ++j)
798 res.setval(irow[j]) += comp[j]*v(
i);
803 template <typename T>
808 template <
typename T>
811 STD__ cout <<
"CSCMatrix<T>::MatVecMult (T*, T*): Not tested!!" <<
STD__ endl;
814 for (i = 0; i < n_cols; ++
i)
817 for (i = 0; i < n_cols; ++
i) {
818 for (
unsigned int j = pcol[i]; j < pcol[i+1]; ++j)
819 res[irow[j]] += comp[j]*v[i];
825 template <
typename T>
831 for (
unsigned int i = 0;
i < n_cols; ++
i) {
833 for (
unsigned int j = pcol[
i]; j < pcol[
i+1]; ++j) {
839 res.
set(val-corr,
i);
842 for (
unsigned int i = 0;
i < n_cols; ++
i) {
844 for (
unsigned int j = pcol[
i]; j < pcol[
i+1]; ++j)
845 val += comp[j]*v(irow[j]);
853 template <
typename T>
861 template <
typename T>
870 template <
typename T>
874 for (
unsigned int i = 0;
i < n_size; ++
i)
881 template <
typename T>
885 for (
unsigned int i = 0;
i < n_size; ++
i)
892 template <typename T>
894 {
return m.
mult (z); }
897 template <
typename T>
900 for (
unsigned int i = 0;
i < n_size; ++
i)
906 template <
typename T>
910 for (
unsigned int i = 0;
i < n_size; ++
i)
917 template <
typename T>
920 for (
unsigned int i = 0;
i < n_size; ++
i)
928 template <
typename T>
930 unsigned int nnzeros)
932 BCHKNR(rows <= 0 || columns <= 0,
CSCMatErr, Zero space allocating, 0);
935 pcol =
new unsigned int[n_cols+1];
937 irow =
new unsigned int[nnzeros];
939 comp =
new T[nnzeros];
943 for (i = 0; i < n_cols+1; ++
i)
945 for (i = 0; i < nnzeros; ++
i)
949 n_max_size = nnzeros;
955 template <
typename T>
972 template <
typename T>
979 pcol =
new unsigned int[n_cols+1];
981 irow =
new unsigned int[n_size];
983 comp =
new T[n_size];
986 for (
unsigned int j = 0; j < n_cols+1; ++j)
988 for (
unsigned int i = 0;
i < n_size; ++
i) {
997 template <
typename T>
1001 return this->
clear();
1005 for (
unsigned r = 0; r <
rows(); ++r)
1006 (*
this)(r,
c) = val;
1011 template <
typename T>
1013 {
CSTD__ memset (comp, 0,
sizeof(
T) * n_size);
return *
this; }
1015 template <
typename T>
1018 BCHK (n_rows != n_cols,
CSCMatErr,
setunit is meaningless
for non-square matrices, n_cols, *
this);
1020 for (
unsigned int i = 0;
i < n_cols;
i++)
1026 template <
typename T>
1036 template <
typename T>
1040 for (
unsigned r = 0; r <
rows(); ++r)
1042 m(
c,r) = (*this)(r,
c);
1046 template <
typename T>
1054 template <typename T>
CSCMatrix(const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
TVector< T > transMult(const Vector_Sig< T > &) const
CSCMatrix(const F_Matrix< T > &m, const double tol=0)
void MatVecMult(Vector< T > &res, const CSCMatrix< T > &m, const Vector< T > &v)
#define BCHKNR(cond, exc, txt, ind)
unsigned int columns() const
unsigned int rows() const
CSCMatErr(const char *t, const long i=0)
unsigned long size() const
unsigned int columns() const
number of columns
unsigned int * rowIndexPointer()
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
CSCMatrix< T > & setval(const T &z, unsigned int row, unsigned int column)
element access (write)
F_TMatrix< T > mult1(const F_Matrix< T > &) const
Calculate *this * m1 (smart version)
CSCMatrix< T > & do_import(const MatType &M)
Import operation, automatic conversion to sparse matrix.
TVector< T > transMult(const Vector< T > &v) const HOT
transpose-vector multiplication
exception base class for the TBCI NumLib
F_TMatrix< T > operator*(const CSCMatrix< T > &) const
CSCMatrix * CSCMatrix.
unsigned int rows() const
Common interface definition (signature) for all Matrices.
const unsigned TMatrix< T > * res
CSCMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
unsigned int columns() const
number of columns
CSCMatrix< T > operator/(const T &z) const
#define BCHK(cond, exc, txt, ind, rtval)
CSCMatErr(const CSCMatErr &ce)
bool operator!=(const CSCMatrix< T > &m) const
TVector< T > cscm_vec_mul_exact(const Vector< T > &V) const
CSCMatrix(const CSCMatrix< T > &m)
const Vector< T > const Vector< T > & x
T get(const unsigned long i) const HOT
F_TMatrix< T > mult(const F_Matrix< T > &) const
Calculate *this * m1 (dumb version)
CSCMatrix< T > & fill(const T &)
set all defined element to a val
Matrix_Sig< T > & clear()
CSCMatrix< T > & clear()
set all elements defined to zero
bool operator==(const CSCMatrix< T > &m) const
matrix-matrix assignment and comparison
for(REGISTER T *p1=c.vec,*p2=b.vec;p1< c.endvec;p1++, p2++)*p1
unsigned long _bin_search(const T *vec, T el, unsigned long start, unsigned long end)
Search for an element el in a sorted vector between start and end-1, returns (unsigned long)-1 if ele...
exception class: Use MatErr from matrix.h
void copy(const CSCMatrix< T > &m)
unsigned long size() const
long int Vector< T > & index
Temporary Base Class Idiom: Class TVector is used for temporary variables.
static const char * mat_info()
allow instantiation (Matrix_Sig)
Temporary object for scaled matrices.
CSCMatrix< T > & operator=(const CSCMatrix< T > &m)
unsigned int rows() const
query matrix dimensions
unsigned int columns() const
CSCMatrix< T > operator-() const
matrix negation
T & setval(const unsigned long i) const
Matrix_Sig< T > & setunit(const T &=(T) 1)
CSCMatrix< T > & operator*=(const T &z)
unsigned int rows() const
number of rows
void SWAP(T &a, T &b)
SWAP function Note: We could implement a swap function without temporaries: a -= b b += a a -= b a = ...
CSCMatrix(const Matrix< T > &m, const double tol=0)
CSCMatrix< T > & setunit(const T &=(T) 1)
set CSCMatrix to val times the unit matrix
unsigned int rows() const
number of rows
Temporary Base Class (non referable!) (acc.
void allocate(unsigned int rows, unsigned int columns, unsigned int nnzeros=1)
#define EXPCHK(cond, exc, txt, ind, rtval)
CSCMatrix< T > & resize(const unsigned int newRows, const unsigned int newColumns, const unsigned int nnzeros=1)
change matrix dimensions
CSCMatrix(const T &val, const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
unsigned int * columnPointer()
unsigned int n_rows
Storage format: pcol holds the offsets of each column; the length of the column c is pcol[c+1]-pcol[c...
F_TMatrix< T > multf(const F_Matrix< T > &) const
matrix-matrix multiplication
CSCMatrix< T > & transpose()
tbci_traits< T >::const_refval_type get(unsigned int row, unsigned int column) const HOT
element access (read)
tbci_traits< T >::const_refval_type operator()(unsigned int row, unsigned int column) const HOT
element access (read)
void setval(const T val, const unsigned int r, const unsigned int c)
CSCMatrix< T > operator+(const CSCMatrix< T > &) const
Addition.
void MatVecMult(Vector< T > &res, const Vector< T > &v) const HOT
for friend void MatVecMult FGD (Vector& res, const CSCMatrix& m, const Vector& v); ...
unsigned int columns() const
void do_export(MatType &M)
Export operation.
unsigned int size() const
T aligned_value_type TALIGN(MIN_ALIGN2)
unsigned int do_exactsum2()
CSCMatrix(const BdMatrix< T > &m, const double tol=0)
T & set(const T &val, const unsigned long i) const
#define EXPCHKNR(cond, exc, txt, ind)
CSCMatrix< T > & operator/=(const T &z)
CSCMatrix< T > transpose(const CSCMatrix< T > &cscm)
CSCMatrix< T > & swap(CSCMatrix< T > &)
void insert(const unsigned int column, const unsigned int pos)