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>
BdMatrix< T > transpose(BdMatrix< T > &mat)
TVector< T > transMult(const Vector_Sig< T > &) const
#define BCHKNR(cond, exc, txt, ind)
const Vector< T > const Vector< T > & x
CSCMatrix< T > & setval(const T &z, unsigned int row, unsigned int column)
element access (write)
T & setval(const unsigned long i) const
CSCMatErr(const char *t, const long i=0)
bool operator!=(const CSCMatrix< T > &m) const
unsigned int columns() const
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
CSCMatrix(const BdMatrix< T > &m, const double tol=0)
bool operator==(const CSCMatrix< T > &m) const
matrix-matrix assignment and comparison
F_TMatrix< T > mult(const F_Matrix< T > &) const
Calculate *this * m1 (dumb version)
TVector< T > transMult(const Vector< T > &v) const HOT
transpose-vector multiplication
CSCMatrix(const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
exception base class for the TBCI NumLib
unsigned int * columnPointer()
unsigned int rows() const
Common interface definition (signature) for all Matrices.
F_TMatrix< T > multf(const F_Matrix< T > &) const
matrix-matrix multiplication
TVector< T > cscm_vec_mul_exact(const Vector< T > &V) const
unsigned int columns() const
number of columns
#define BCHK(cond, exc, txt, ind, rtval)
CSCMatErr(const CSCMatErr &ce)
CSCMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
CSCMatrix< T > & operator*=(const T &z)
unsigned int rows() const
query matrix dimensions
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.
unsigned int size() const
void MatVecMult(Vector< T > &res, const Vector< T > &v) const HOT
for friend void MatVecMult FGD (Vector<T>& res, const CSCMatrix<T>& m, const Vector<T>& v); ...
CSCMatrix< T > & operator=(const CSCMatrix< T > &m)
T & set(const T &val, const unsigned long i) const
static const char * mat_info()
allow instantiation (Matrix_Sig)
Temporary object for scaled matrices.
Matrix_Sig< T > & clear()
F_TMatrix< T > operator*(const CSCMatrix< T > &) const
CSCMatrix * CSCMatrix.
CSCMatrix(const CSCMatrix< T > &m)
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > & z
CSCMatrix< T > & transpose()
CSCMatrix< T > operator/(const T &z) const
void setval(const T val, const unsigned int r, const unsigned int c)
T get(const unsigned long i) const HOT
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...
CSCMatrix< T > & clear()
set all elements defined to zero
unsigned int columns() const
number of columns
unsigned int * rowIndexPointer()
void insert(const unsigned int column, const unsigned int pos)
long int Vector< T > & index
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...
CSCMatrix(const F_Matrix< T > &m, const double tol=0)
CSCMatrix< T > operator+(const CSCMatrix< T > &) const
Addition.
void allocate(unsigned int rows, unsigned int columns, unsigned int nnzeros=1)
CSCMatrix< T > & operator/=(const T &z)
void copy(const CSCMatrix< T > &m)
CSCMatrix(const Matrix< T > &m, const double tol=0)
unsigned int columns() const
unsigned int rows() const
Matrix_Sig< T > & setunit(const T &=(T) 1)
CSCMatrix< T > & fill(const T &)
set all defined element to a val
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 = ...
tbci_traits< T >::const_refval_type operator()(unsigned int row, unsigned int column) const HOT
element access (read)
unsigned long size() const
#define EXPCHK(cond, exc, txt, ind, rtval)
Temporary Base Class (non referable!) (acc.
tbci_traits< T >::const_refval_type get(unsigned int row, unsigned int column) const HOT
element access (read)
exception class: Use MatErr from matrix.h
Temporary Base Class Idiom: Class TVector is used for temporary variables.
unsigned long size() const
void MatVecMult(Vector< T > &res, const CRMatrix< T > &m, const Vector< T > &v)
const Vector< T > Vector< T > Vector< T > Vector< T > & y
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int res
CSCMatrix< T > operator-() const
matrix negation
void do_export(MatType &M)
Export operation.
unsigned int do_exactsum2()
CSCMatrix< T > & resize(const unsigned int newRows, const unsigned int newColumns, const unsigned int nnzeros=1)
change matrix dimensions
CSCMatrix< T > & swap(CSCMatrix< T > &)
unsigned int columns() const
#define EXPCHKNR(cond, exc, txt, ind)
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int int char v
< find minimun of func on grid with resolution res
T aligned_value_type TALIGN(MIN_ALIGN2)
CSCMatrix(const T &val, const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
CSCMatrix< T > & setunit(const T &=(T) 1)
set CSCMatrix to val times the unit matrix
unsigned int rows() const
number of rows