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"
77 const unsigned int columns,
const unsigned int nnzeros=1)
79 if (val != (
T)0)
fill(val); }
81 const unsigned int nnzeros=1)
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)
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)
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) <<
" " ;
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) {
288 res.setval(val,
i, j);
297{
return m2.
multf (m1); }
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);
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);
350{
return m1.
mult1 (m2); }
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);
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]) {
401 }
else if (
irow[off0] == m1.
irow[off1]) {
405 }
else if (
irow[off0] < m1.
irow[off1]) {
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]) {
440 }
else if (
irow[off0] < m1.
irow[off1]) {
454inline typename tbci_traits<T>::const_refval_type
510 if ((epos - spos) >= 32) {
511 const unsigned long int off = _bin_search(
irow,
row, spos, epos);
512 if (off != (
unsigned long)(-1))
519 for (
REGISTER unsigned int pos = spos; pos < epos; ++pos) {
526 for (
REGISTER unsigned int pos = epos; pos > spos; --pos) {
527 const unsigned po1 = pos-1;
539inline typename tbci_traits<T>::const_refval_type
553 for (
unsigned int off =
pcol[
c]; off <
pcol[
c+1]; ++off)
562 const unsigned int newColumns,
563 const unsigned int nnzeros)
567 allocate(newRows, newColumns, nnzeros);
570 unsigned int *newirow =
new unsigned int[nnzeros];
571 T *newcomp =
new T [nnzeros];
574 for (
i = 0;
i < com_sz; ++
i) {
578 for (;
i < nnzeros; ++
i)
596 unsigned int* new_irow =
new unsigned int [
n_max_size];
604 for (
i = 0;
i < pos; ++
i) {
627 for (
int j =
n_size-1; j >= (int)pos; --j) {
633 for (
unsigned int i = column+1;
i <
n_cols+1; ++
i)
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)
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))
702 for (
unsigned int i = 0;
i <
n_size; ++
i)
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;
737 for (
unsigned int i = 0;
i <
n_cols; ++
i) {
738 for (
unsigned int j =
pcol[
i]; j <
pcol[
i+1]; ++j) {
759 for (
unsigned int i = 0;
i <
n_cols; ++
i) {
760 for (
unsigned int j =
pcol[
i]; j <
pcol[
i+1]; ++j)
770 return (this->
operator * (v));
780 for (
unsigned int i = 0;
i <
n_cols; ++
i) {
781 for (
unsigned int j =
pcol[
i]; j <
pcol[
i+1]; ++j)
796 for (
unsigned int i = 0;
i <
n_cols; ++
i) {
797 for (
unsigned int j =
pcol[
i]; j <
pcol[
i+1]; ++j)
811 STD__ cout <<
"CSCMatrix<T>::MatVecMult (T*, T*): Not tested!!" <<
STD__ endl;
818 for (
unsigned int j =
pcol[
i]; j <
pcol[
i+1]; ++j)
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)
874 for (
unsigned int i = 0;
i <
n_size; ++
i)
885 for (
unsigned int i = 0;
i <
n_size; ++
i)
894{
return m.
mult (z); }
900 for (
unsigned int i = 0;
i <
n_size; ++
i)
910 for (
unsigned int i = 0;
i <
n_size; ++
i)
920 for (
unsigned int i = 0;
i <
n_size; ++
i)
930 unsigned int nnzeros)
937 irow =
new unsigned int[nnzeros];
939 comp =
new T[nnzeros];
945 for (
i = 0;
i < nnzeros; ++
i)
986 for (
unsigned int j = 0; j <
n_cols+1; ++j)
988 for (
unsigned int i = 0;
i <
n_size; ++
i) {
1001 return this->
clear();
1005 for (
unsigned r = 0; r <
rows(); ++r)
1006 (*
this)(r,
c) = val;
1011template <
typename T>
1015template <
typename T>
1020 for (
unsigned int i = 0;
i <
n_cols;
i++)
1026template <
typename T>
1036template <
typename T>
1040 for (
unsigned r = 0; r <
rows(); ++r)
1042 m(
c,r) = (*this)(r,
c);
1046template <
typename T>
1054template <typename T>
long int Vector< T > & index
const Vector< T > const Vector< T > & x
#define BCHK(cond, exc, txt, ind, rtval)
#define EXPCHK(cond, exc, txt, ind, rtval)
#define BCHKNR(cond, exc, txt, ind)
#define EXPCHKNR(cond, exc, txt, ind)
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
unsigned int rows() const
number of rows
unsigned int columns() const
number of columns
CSCMatErr(const CSCMatErr &ce)
CSCMatErr(const char *t, const long i=0)
exception class: Use MatErr from matrix.h
CSCMatrix(const Matrix< T > &m, const double tol=0)
CSCMatrix< T > & transpose()
CSCMatrix< T > & fill(const T &)
set all defined element to a val
CSCMatrix< T > operator+(const CSCMatrix< T > &) const
Addition.
F_TMatrix< T > operator*(const CSCMatrix< T > &) const
CSCMatrix * CSCMatrix.
CSCMatrix< T > operator-() const
matrix negation
CSCMatrix< T > & clear()
set all elements defined to zero
CSCMatrix(const F_Matrix< T > &m, const double tol=0)
F_TMatrix< T > mult1(const F_Matrix< T > &) const
Calculate *this * m1 (smart version).
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);
F_TMatrix< T > multf(const F_Matrix< T > &) const
matrix-matrix multiplication
void copy(const CSCMatrix< T > &m)
CSCMatrix(const T &val, const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
CSCMatrix< T > & resize(const unsigned int newRows, const unsigned int newColumns, const unsigned int nnzeros=1)
change matrix dimensions
CSCMatrix< T > & setval(const T &z, unsigned int row, unsigned int column)
element access (write)
CSCMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
bool operator==(const CSCMatrix< T > &m) const
matrix-matrix assignment and comparison
CSCMatrix< T > & operator/=(const T &z)
CSCMatrix< T > & setunit(const T &=(T) 1)
set CSCMatrix to val times the unit matrix
unsigned int size() const
void do_export(MatType &M)
Export operation.
CSCMatrix< T > & operator*=(const T &z)
CSCMatrix(const CSCMatrix< T > &m)
void insert(const unsigned int column, const unsigned int pos)
CSCMatrix< T > & swap(CSCMatrix< T > &)
TVector< T > transMult(const Vector< T > &v) const HOT
transpose-vector multiplication
tbci_traits< T >::const_refval_type get(unsigned int row, unsigned int column) const HOT
element access (read)
unsigned int columns() const
static const char * mat_info()
allow instantiation (Matrix_Sig)
CSCMatrix(const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
T aligned_value_type TALIGN(MIN_ALIGN2)
F_TMatrix< T > mult(const F_Matrix< T > &) const
Calculate *this * m1 (dumb version).
CSCMatrix(const BdMatrix< T > &m, const double tol=0)
CSCMatrix< T > & operator=(const CSCMatrix< T > &m)
TVector< T > cscm_vec_mul_exact(const Vector< T > &V) const
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...
bool operator!=(const CSCMatrix< T > &m) const
CSCMatrix< T > & do_import(const MatType &M)
Import operation, automatic conversion to sparse matrix.
CSCMatrix< T > operator/(const T &z) const
unsigned int rows() const
query matrix dimensions
unsigned int * rowIndexPointer()
void allocate(unsigned int rows, unsigned int columns, unsigned int nnzeros=1)
tbci_traits< T >::const_refval_type operator()(unsigned int row, unsigned int column) const HOT
element access (read)
unsigned int * columnPointer()
unsigned int columns() const
unsigned int rows() const
Temporary Base Class (non referable!) (acc.
Temporary object for scaled matrices.
T get(const unsigned r, const unsigned c) const
unsigned int rows() const
TVector< T > transMult(const Vector_Sig< T > &) const
Matrix_Sig< T > & clear()
unsigned int columns() const
unsigned int rows() const
number of rows
unsigned int columns() const
number of columns
T get(const unsigned long i) const HOT
unsigned long size() const
Temporary Base Class Idiom: Class TVector is used for temporary variables.
unsigned long size() const
T & setval(const unsigned long i) const
double fabs(const TBCI__ cplx< T > &c)
STD__ ostream & operator<<(STD__ ostream &stream, const CSCMatrix< T > &m)
CSCMatrix< T > transpose(const CSCMatrix< T > &cscm)
void MatVecMult(Vector< T > &res, const CSCMatrix< T > &m, const Vector< T > &v)
return F_TMatrix< T >(tm)
const unsigned TMatrix< T > * res
unsigned int do_exactsum2()