TBCI Numerical high perf. C++ Library: cscmatrix.h Source File

00001 
00005 //----------------------------------------------------------
00006 // version 1.0, 08.1998
00007 //
00008 // Andreas Ahland
00009 //
00010 // $Id: cscmatrix.h,v 1.18.2.29 2009/12/05 08:35:30 garloff Exp $
00011 //----------------------------------------------------------
00012 
00013 #ifndef TBCI_CSCMATRIX_H
00014 #define TBCI_CSCMATRIX_H
00015 
00016 #include "vector.h"
00017 #include "matrix_sig.h"
00018 #include "f_matrix.h"
00019 #include "band_matrix.h"
00020 
00021 // Avoid -fguiding-decls
00022 #if !defined(NO_GD) && !defined(AUTO_DECL)
00023 # include "cscmatrix_gd.h"
00024 #endif
00025 
00026 NAMESPACE_TBCI
00027 
00028 // exception class
00029 #ifdef EXCEPT
00030 //# include "except.h" is include in basics.h
00031 
00032 class CSCMatErr : public NumErr
00033 {
00034  public:
00035         CSCMatErr() 
00036         : NumErr("Error in CSCMatrix") {}
00037         CSCMatErr(const char* t, const long i = 0) 
00038         : NumErr(t, i) {}
00039         CSCMatErr(const CSCMatErr& ce)
00040         : NumErr(ce) {}
00041         virtual ~CSCMatErr() throw() {}
00042 };
00043 #endif
00044 
00045 #ifdef PRAGMA_I
00046 # pragma interface "cscmatrix.h"
00047 #endif
00048 
00049 template <typename T> class CSCMatrix;
00050 template <typename T> class F_TSMatrix;
00051 template <typename T> class F_TMatrix;
00052 template <typename T> class F_Matrix;
00053 
00054 //----------
00055 // CSCMatrix
00056 //----------
00064 template <typename T>
00065 class CSCMatrix : public Matrix_Sig<T>
00066 {
00067         friend class F_TSMatrix<T>;
00068 
00069 public:
00070         typedef T value_type;
00071         typedef T element_type;
00072         typedef T aligned_value_type TALIGN(MIN_ALIGN2);
00073 
00075         CSCMatrix() { allocate(1,1); }
00076         CSCMatrix(const T& val, const unsigned int rows, 
00077                           const unsigned int columns, const unsigned int nnzeros=1) 
00078         {   allocate(rows, columns, nnzeros); 
00079                 if (val != (T)0) fill(val); }
00080         CSCMatrix(const unsigned int rows, const unsigned int columns, 
00081                           const unsigned int nnzeros=1) 
00082         { allocate(rows, columns, nnzeros); }
00083         CSCMatrix(const CSCMatrix<T>& m) 
00084         { copy(m); }
00085         CSCMatrix(const F_Matrix<T>& m, const double tol=0)
00086         { allocate(m.rows(), m.columns());      fill(m,tol); }
00087         CSCMatrix(const BdMatrix<T>& m, const double tol=0)
00088         { allocate(m.rows(), m.columns());      fill(m,tol); }
00089         CSCMatrix(const   Matrix<T>& m, const double tol=0)
00090         { allocate(m.rows(), m.columns());      fill(m,tol); }
00091         ~CSCMatrix() 
00092         { destroy (); }
00093         
00094         operator F_TMatrix<T> () const;
00095 
00097         /*virtual*/ static const char* mat_info ()  { return ("CSCMatrix"); }
00098 
00099 
00101         unsigned int rows()    const  { return n_rows; }
00102         unsigned int columns() const  { return n_cols; }
00103         unsigned int size()    const  { return n_size; } // number of non-zero elements
00104 
00105 
00107         CSCMatrix<T>& resize(const unsigned int newRows, 
00108                              const unsigned int newColumns, 
00109                              const unsigned int nnzeros = 1);
00110 
00112         CSCMatrix<T>& clear ();
00114         CSCMatrix<T>& fill (const T&);
00116         CSCMatrix<T>& setunit (const T& = (T)1);
00117 
00119 #ifdef PREFER_COPY
00120         T operator() (unsigned int row, unsigned int column) const HOT;
00121         T get (unsigned int row, unsigned int column) const
00122         { return (*this) (row, column); }
00123 #else
00124         const T& operator() (unsigned int row, unsigned int column) const HOT;
00125         const T& get (unsigned int row, unsigned int column) const
00126         { return this->operator() (row, column); }
00127 #endif
00129         CSCMatrix<T>& setval(const T& z, unsigned int row, unsigned int column);
00130         T& setval(unsigned int row, unsigned int column);
00131         T& operator () (unsigned int row, unsigned int column)
00132         { return setval (row, column); }
00133         T* dataPointer()                { return comp; }
00134         unsigned int* columnPointer()   { return pcol; }
00135         unsigned int* rowIndexPointer() { return irow; }
00136 
00138         bool operator == (const CSCMatrix<T>& m) const;
00139         bool operator != (const CSCMatrix<T>& m) const    
00140         { return !(this->operator == (m)); }
00141         CSCMatrix<T>& operator= (const CSCMatrix<T>& m);
00142 
00144         //CSCMatrix<T>& do_import (const Matrix_Sig<T>& M)
00145         template <typename MatType>
00146         CSCMatrix<T>& do_import  /*FGD*/ (const MatType& M)
00147         {
00148                 destroy (); allocate(M.rows(), M.columns(), 20);
00149                 fill (M);
00150                 return *this;
00151         }
00152         
00154         //void do_export (Matrix_Sig<T>& M)
00155         template <typename MatType>
00156         void do_export /*FGD*/ (MatType& M)
00157         {
00158                 M.clear (); //Matrix should be cleared!
00159                 for (unsigned int i=0; i<n_cols; ++i)
00160                         //for (unsigned int j=pcol[i]; j<pcol[i+1]; ++j)  M.setval(comp[j], i,irow[j]);
00161                         for (unsigned int j=pcol[i]; j<pcol[i+1]; ++j)  
00162                                 M(i,irow[j]) = comp[j];
00163         }
00164 
00166         CSCMatrix<T> operator - () const;
00167 
00169         CSCMatrix<T> operator + (const CSCMatrix<T>&) const;
00170         CSCMatrix<T> operator - (const CSCMatrix<T>&) const;
00171 
00173         F_TMatrix<T> operator * (const CSCMatrix<T>&) const;
00174 
00176         //friend  F_TMatrix<T> operator* FGD (const F_Matrix<T>& m1, const CSCMatrix<T>& m2);
00177         F_TMatrix<T> multf (const F_Matrix<T>&) const;
00178         F_TMatrix<T> mult  (const F_Matrix<T>&) const;
00179         F_TMatrix<T> mult1 (const F_Matrix<T>&) const;
00180 
00182         TVector<T> operator * (const Vector<T>& v) const HOT;
00183         TVector<T> operator * (const TVector<T>& tv) const HOT;
00184         TVector<T> operator * (TSVector<T>& tsv) const HOT;
00186         void MatVecMult (Vector<T>& res, const Vector<T>& v) const HOT;
00187         void MatVecMult (T* v, T* res) HOT; // res <- A*v , for ARPACK++
00188 
00190         TVector<T> transMult(const Vector<T>& v) const HOT;
00191         TVector<T> transMult(const TVector<T>& tv) const HOT;
00192         TVector<T> transMult(const TSVector<T>& tsv) const HOT;
00193 
00194 
00195         CSCMatrix<T> operator* (const T& z) const;       // matrix-scalar multiplication
00197         CSCMatrix<T> mult (const T&) const;
00198 
00199         CSCMatrix<T>& operator*= (const T& z);         // fast in-place matrix-scalar multiplication
00200         CSCMatrix<T> operator/ (const T& z) const;       // matrix-scalar division
00201         CSCMatrix<T>& operator/= (const T& z);       // fast in-place matrix-scalar division
00202 
00203         CSCMatrix<T>& swap (CSCMatrix<T>&);
00205         CSCMatrix<T>  transposed_copy() const;
00206         CSCMatrix<T>& transpose();
00207 
00209         friend STD__ ostream& operator<< FGDT (STD__ ostream& stream, const CSCMatrix<T>& m);
00210 
00217 protected:
00218         unsigned int n_rows, n_cols, n_size, n_max_size;
00219         unsigned int* pcol;  // pointer to the beginning of each column
00220         unsigned int* irow;  // row index of all nonzero elements
00221         T* comp;             // storage vector
00222         mutable T dummy;
00223 
00224         void allocate (unsigned int rows, unsigned int columns, unsigned int nnzeros=1);
00225         void destroy  ();       // used by destructor and "=" operator
00226         void copy     (const CSCMatrix<T>& m);  // used by copy constr. and "="
00227         void insert   (const unsigned int column, const unsigned int pos); //dynamic inserting
00228 
00230         //void fill (const Matrix_Sig<T>& M, double tol=0)
00231         template <typename MatType>
00232         void fill /*FGD*/ (const MatType& M, const double tol = 0.0)
00233         {
00234                 unsigned int count = 0, i, j;
00235                 //count nonzeros
00236                 for (j = 0; j < n_cols; ++j)
00237                         for (i = 0; i < n_rows; ++i)    
00238                                 if (MATH__ fabs(M(i,j)) > tol) 
00239                                         ++count; 
00240                 STD__ cout << "Matrix conversion, NNZ=" << count << STD__ endl;
00241                 resize(M.rows(), M.columns(), count);
00242                 for (j = 0; j < n_cols; ++j)
00243                         for (i = 0; i < n_rows; ++i)  
00244                                 if (MATH__ fabs(M(i,j)) > tol) 
00245                                         setval(M(i,j), i,j); 
00246         }
00247 };
00248 
00249 
00250 template <typename T>
00251 STD__ ostream& operator<< (STD__ ostream& stream, const CSCMatrix<T>& m)
00252 {
00253         for (unsigned int i = 0; i < m.rows(); ++i) {
00254                 //stream << "|\t";
00255                 for (unsigned int j = 0; j < m.columns(); ++j)
00256                         stream << m(i,j) << " " /* << "\t"*/;
00257                 stream /*<< "|" */ << "\n";
00258         }
00259         return stream /* << STD__ flush*/;
00260 }
00261 
00262 //-------------------------------
00263 // CSCMatrix function definitions 
00264 //-------------------------------
00265 
00267 template <typename T>
00268 F_TMatrix<T> CSCMatrix<T>::multf (const F_Matrix<T>& m1) const
00269 {       
00270         F_TMatrix<T> res((T)0, m1.rows(), columns());
00271         BCHK(rows() != m1.columns(), CSCMatErr, Sizes dont match in multf, m1.columns(), res);
00272         for (unsigned int i = 0; i < m1.rows(); ++i)
00273                 for (unsigned int j = 0; j < columns(); ++j) {
00274                         T val = 0;
00275                         for (unsigned int x = pcol[j]; x < pcol[j+1]; ++x)
00276                                 val += m1(i, irow[x] )*comp[x];
00277                         res.setval(val, i, j);
00278                 }
00279          return res;
00280 }
00281 
00282 INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const F_Matrix<T>&, const CSCMatrix<T>&);)
00283 template <typename T>
00284 inline F_TMatrix<T> operator * (const F_Matrix<T>& m1, const CSCMatrix<T>& m2)
00285 { return m2.multf (m1); }
00286 
00288 template <typename T>
00289 F_TMatrix<T> CSCMatrix<T>::mult (const F_Matrix<T>& m1) const
00290 {       
00291         F_TMatrix<T> res((T)0, rows(), m1.columns());
00292         BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
00293         for (unsigned int i = 0; i < rows(); ++i)
00294                 for (unsigned int j = 0; j < m1.columns(); ++j) {
00295                         T val = 0;
00296                         for (unsigned int x = 0; x < columns(); ++x) 
00297                                 val += (*this)(i, x) * m1(x, j);
00298                         res.setval(val, i, j);
00299                 }
00300         return res;
00301 }
00302 
00304 template <typename T>
00305 F_TMatrix<T> CSCMatrix<T>::mult1 (const F_Matrix<T>& m1) const
00306 {       
00307         F_TMatrix<T> res((T)0, rows(), m1.columns());
00308         BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
00309         for (unsigned int c = 0; c < columns(); ++c) 
00310                 for (unsigned int off = pcol[c]; off < pcol[c+1]; ++off) {
00311                         unsigned int r = irow[off];
00312                         const T el = comp[off]; // (*this)(r,c)
00313                         for (unsigned int j = 0; j < m1.columns(); ++j)
00314                                 res.setval(r, j) += el * m1(c, j);
00315                 }
00316         return res;
00317 }
00318 
00319 INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const CSCMatrix<T>&, const F_Matrix<T>&);)
00320 template <typename T>
00321 inline F_TMatrix<T> operator * (const CSCMatrix<T>& m1, const F_Matrix<T>& m2)
00322 { return m1.mult1 (m2); }
00323 
00325 template <typename T>
00326 F_TMatrix<T> CSCMatrix<T>::operator * (const CSCMatrix<T>& m1) const
00327 {
00328         F_TMatrix<T> res((T)0, rows(), m1.columns());
00329         BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
00330         for (unsigned int r0 = 0; r0 < rows(); ++r0)
00331                 for (unsigned int c1 = 0; c1 < m1.columns(); ++c1) {
00332                         T val = 0;
00333                         for (unsigned int off1 = m1.pcol[c1]; off1 < m1.pcol[c1+1]; ++off1)
00334                                 val += (*this)(r0, m1.irow[off1]) * m1.comp[off1];
00335                         res.setval(val, r0, c1);
00336                 }
00337         return res;
00338 }
00339 
00341 template <typename T>
00342 CSCMatrix<T> CSCMatrix<T>::operator + (const CSCMatrix<T>& m1) const
00343 {
00344         CSCMatrix<T> res(rows(), columns(), size());
00345         BCHK(columns() != m1.columns(), CSCMatErr, Columns dont match in add, columns(), res);
00346         BCHK(rows() != m1.rows(), CSCMatErr, Rows dont match in add, rows(), res);
00347         for (unsigned int c = 0; c < columns(); ++c) {
00348                 unsigned int off0 = pcol[c];
00349                 unsigned int off1 = m1.pcol[c];
00350                 while (off0 < pcol[c+1] || off1 < m1.pcol[c+1]) {
00351                         // only one is valid ...
00352                         if (off0 >= pcol[c+1]) {
00353                                 res(m1.irow[off1], c) = m1.comp[off1]; 
00354                                 ++off1;
00355                         } else if (off1 >= m1.pcol[c+1]) {
00356                                 res(irow[off0], c) = comp[off0]; 
00357                                 ++off0;
00358                         } else if (irow[off0] == m1.irow[off1]) {
00359                                 // both are valid
00360                                 res(irow[off0], c) = comp[off0] + m1.comp[off1];
00361                                 ++off0; ++off1;
00362                         } else if (irow[off0] < m1.irow[off1]) {
00363                                 res(irow[off0], c) = comp[off0];
00364                                 ++off0;
00365                         } else {
00366                                 res(m1.irow[off1], c) = m1.comp[off1];
00367                                 ++off1;
00368                         }
00369                 }
00370         }
00371         return res;
00372 }
00373 
00374 template <typename T>
00375 CSCMatrix<T> CSCMatrix<T>::operator - (const CSCMatrix<T>& m1) const
00376 {
00377         CSCMatrix<T> res(rows(), columns(), size());
00378         BCHK(columns() != m1.columns(), CSCMatErr, Columns dont match in sub, columns(), res);
00379         BCHK(rows() != m1.rows(), CSCMatErr, Rows dont match in sub, rows(), res);
00380         for (unsigned int c = 0; c < columns(); ++c) {
00381                 unsigned int off0 = pcol[c];
00382                 unsigned int off1 = m1.pcol[c];
00383                 while (off0 < pcol[c+1] || off1 < m1.pcol[c+1]) {
00384                         // only one is valid ...
00385                         if (off0 >= pcol[c+1]) {
00386                                 res(m1.irow[off1], c) = -m1.comp[off1]; 
00387                                 ++off1; continue;
00388                         }
00389                         if (off1 >= m1.pcol[c+1]) {
00390                                 res(irow[off0], c) = comp[off0]; 
00391                                 ++off0; continue;
00392                         }
00393                         // both are valid
00394                         if (irow[off0] == m1.irow[off1]) {
00395                                 res(irow[off0], c) = comp[off0] - m1.comp[off1];
00396                                 ++off0; ++off1; //continue;
00397                         } else if (irow[off0] < m1.irow[off1]) {
00398                                 res(irow[off0], c) = comp[off0];
00399                                 ++off0; //continue;
00400                         } else {
00401                                 res(m1.irow[off1], c) = -m1.comp[off1];
00402                                 ++off1; //continue;
00403                         }
00404                 }
00405         }
00406         return res;
00407 }
00408 
00409 
00410 
00412 template <typename T>
00413 #ifdef PREFER_COPY
00414 inline T CSCMatrix<T>::operator() (unsigned int row, unsigned int col) const
00415 #else
00416 inline const T& CSCMatrix<T>::operator() (unsigned int row, 
00417                                                                                   unsigned int col) const
00418 #endif
00419 {
00420         static T dummy = (T)0;
00421         BCHK(row >= n_rows, CSCMatErr, Illegal row index, row, dummy);
00422         BCHK(col >= n_cols, CSCMatErr, Illegal col index, col, dummy);
00423 /*
00424         // Version 1 
00425         for (unsigned int j=pcol[col]; j<pcol[col+1]; j++)
00426         {
00427                 if (irow[j] == row) return comp[j];
00428         }
00429 
00430         // Version 2
00431         unsigned int pos = pcol[col]; 
00432         // linear search ...
00433         while (pos < pcol[col+1] && irow[pos] < row) 
00434                 ++pos;
00435         if (pos < n_size && irow[pos] == row) 
00436                 return comp[pos];
00437         // Version 3
00438 
00439         unsigned int pos = pcol[col];
00440         // linear search ...
00441         while (pos < pcol[col+1]) {
00442                 if (irow[pos] < row) 
00443                         ++pos;
00444                 else if (irow[pos] > row)
00445                         break;
00446                 else
00447                         return comp[pos];
00448         }
00449         // Element not found
00450         return (dummy = 0);
00451 */
00452 
00453         // Version 4
00454         const unsigned int epos = pcol[col+1];
00455         for (unsigned int pos = pcol[col]; pos < epos; ++pos) {
00456                 if (irow[pos] == row)
00457                         return comp[pos];
00458                 else if (irow[pos] > row)
00459                         break;
00460         }
00461         return (dummy = 0);
00462 }
00463 
00465 template <typename T>
00466 CSCMatrix<T>::operator F_TMatrix<T> () const
00467 {
00468         F_TMatrix<T> res((T)0, rows(), columns());
00469         for (unsigned int c = 0; c < columns(); ++c) 
00470                 for (unsigned int off = pcol[c]; off < pcol[c+1]; ++off)
00471                         res.setval(comp[off], irow[off], c);
00472         return res;
00473 }
00474 
00475 
00477 template <typename T>
00478 CSCMatrix<T>& CSCMatrix<T>::resize(const unsigned int newRows, 
00479                                    const unsigned int newColumns,
00480                                    const unsigned int nnzeros)
00481 {
00482         if (UNLIKELY(newRows != n_rows && newColumns != n_cols)) {
00483                 destroy ();
00484                 allocate(newRows, newColumns, nnzeros);
00485         } else if (UNLIKELY(nnzeros != n_max_size)) {
00486                 // TODO: Handle size 0
00487                 unsigned int *newirow = new unsigned int[nnzeros];
00488                 T            *newcomp = new T           [nnzeros];
00489                 unsigned int i;
00490                 const unsigned int com_sz = MIN(nnzeros, n_max_size);
00491                 for (i = 0; i < com_sz; ++i) {
00492                         newirow[i] = irow[i];
00493                         newcomp[i] = comp[i];
00494                 }
00495                 for (; i < nnzeros; ++i)
00496                         newirow[i] = 0;
00497                 n_max_size = nnzeros;
00498                 delete[] irow;
00499                 delete[] comp;
00500                 irow = newirow;
00501                 comp = newcomp;
00502         } /* else do nothing */
00503         return *this;
00504 }
00505 
00506 
00507 template <typename T>
00508 void CSCMatrix<T>::insert(const unsigned int column, const unsigned int pos)
00509 {
00510         if (UNLIKELY(n_size >= n_max_size)) {
00511                 //not enough space left
00512                 n_max_size += 16;
00513                 unsigned int* new_irow = new unsigned int [n_max_size];
00514                 BCHKNR(new_irow==NULL, CSCMatErr, Out of memory, n_max_size);
00515                 T*            new_comp = new T            [n_max_size];
00516                 BCHKNR(new_comp==NULL, CSCMatErr, Out of memory, n_max_size);
00517 
00518                 // if data exist, copy data
00519                 if (n_size > 0) {
00520                         unsigned int i;
00521                         for (i = 0; i < pos; ++i)       { 
00522                                 // copy data up before
00523                                 new_irow[i] = irow[i];
00524                                 new_comp[i] = comp[i];
00525                         }
00526                         for (i = pos; i < n_size; ++i) {
00527                                 // copy data behind insertion
00528                                 new_irow[i+1] = irow[i];
00529                                 new_comp[i+1] = comp[i];
00530                         }
00531                         // Following elements are not yet used
00532                         for (i = n_size+1; i < n_max_size; ++i) 
00533                                 new_irow[i] = 0;
00534                 }
00535 
00536                 if (irow != NULL) 
00537                         delete[] irow;
00538                 if (comp != NULL) 
00539                         delete[] comp;
00540                 irow = new_irow;
00541                 comp = new_comp;
00542         } else if (n_size > 0) { 
00543                 // if data exist, copy data
00544                 for (int j = n_size-1; j >= (int)pos; --j) {
00545                          // copy data behind insertion
00546                         irow[j+1] = irow[j];
00547                         comp[j+1] = comp[j];
00548                 }
00549         }
00550         for (unsigned int i = column+1; i < n_cols+1; ++i)
00551                 ++pcol[i];
00552         ++n_size;
00553 }
00554 
00555 // element access (write)
00556 template <typename T>
00557 CSCMatrix<T>& CSCMatrix<T>::setval(const T& z, unsigned int row, unsigned int column)
00558 {
00559         setval(row, column) = z;
00560         return *this;
00561 }
00562 
00564 template <typename T>
00565 T& CSCMatrix<T>::setval(unsigned int row, unsigned int column)
00566 {
00567         BCHKNR(row >= n_rows, CSCMatErr, Illegal row index, row);
00568         BCHKNR(column >= n_cols, CSCMatErr, Illegal col index, column);
00569         unsigned int pos  = pcol[column]; 
00570         const unsigned int posn = pcol[column+1]; 
00571         while (pos < posn && irow[pos] < row) 
00572                 ++pos;
00573         if (pos < posn && irow[pos] == row) 
00574                 return comp[pos];
00575 
00576         // element not found -- dynamically growing required...
00577         insert(column, pos);
00578 
00579         // Write new element
00580         irow[pos] = row;
00581         comp[pos] = T(0);   // Init with zero!!
00582 
00583         return comp[pos];
00584 }
00585 
00586 
00587 // matrix-matrix assignment
00588 template <typename T>
00589 CSCMatrix<T>&   CSCMatrix<T>::operator= (const CSCMatrix<T>& m)
00590 {
00591         if (this == &m) 
00592                 return *this;
00593         destroy(); copy(m);
00594         return *this;
00595 }
00596 
00597 
00598 // KG 19980721
00599 template <typename T>
00600 bool CSCMatrix<T>::operator == (const CSCMatrix<T>& cm) const
00601 {
00602          if (n_rows != cm.n_rows || n_cols != cm.n_cols) 
00603                  return false;
00604          for (unsigned int i = 0; i < n_rows; ++i) {
00605                  for (unsigned int j = 0; j < n_cols; ++j) {
00606                          if (operator()(i,j) != cm(i,j)) 
00607                                  return false;
00608                  }
00609          }
00610          return true;
00611 }
00612 
00613 
00614 // matrix negation
00615 template <typename T>
00616 CSCMatrix<T> CSCMatrix<T>::operator - () const
00617 {
00618         CSCMatrix<T> res(n_rows, n_cols);
00619         for (unsigned int i = 0; i < n_size; ++i) 
00620                 res.comp[i] = -comp[i];
00621         return res; // avoid copy constructor soon?
00622 }
00623 
00624 
00625 // matrix-vector multiplication
00626 template <typename T>
00627 TVector<T> CSCMatrix<T>::operator * (const Vector<T>& v) const
00628 {
00629         TVector<T> res(0, n_rows);
00630         BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
00631         for (unsigned int i = 0; i < n_cols; ++i) {
00632                 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
00633                         res.setval(irow[j]) += comp[j]*v(i);  
00634         }
00635         return res;
00636 }
00637 
00638 template <typename T>
00639 inline TVector<T> CSCMatrix<T>::operator * (const TVector<T>& tv) const
00640 {
00641         Vector<T> v(tv);
00642         return (this->operator * (v));
00643 }
00644 
00645 template <typename T>
00646 inline TVector<T> CSCMatrix<T>::operator * (TSVector<T>& tsv) const
00647 {
00648         TVector<T> res(0, n_rows);
00649         BCHK(n_cols != tsv.size(), CSCMatErr, Dimension conflict, tsv.size(), res);
00650         for (unsigned int i = 0; i < n_cols; ++i) {
00651                 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
00652                         res.setval(irow[j]) += comp[j] * tsv.get(i);
00653         }
00654         tsv.destroy ();
00655         return res;
00656 }
00657 
00658 
00659 // matrix-vector multiplication
00660 template <typename T>
00661 void CSCMatrix<T>::MatVecMult(Vector<T>& res, const Vector<T>& v) const
00662 {
00663         BCHKNR(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size());
00664         res = T(0);
00665         for (unsigned int i = 0; i < n_cols; ++i) {
00666                 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
00667                         res.setval(irow[j]) += comp[j]*v(i);  
00668         }
00669 }
00670 
00671 INST(template <typename T> class CSCMatrix<T> friend void MatVecMult(Vector<T>&, const CSCMatrix<T>&, const Vector<T>&);)
00672 template <typename T>
00673 inline void MatVecMult(Vector<T>& res, const CSCMatrix<T>& m, const Vector<T>& v)
00674 { m.MatVecMult (res, v); }
00675 
00676 
00677 template <typename T>
00678 void CSCMatrix<T>::MatVecMult (T* v, T* res) // res <- A*v , for ARPACK++
00679 {
00680         STD__ cout << "CSCMatrix<T>::MatVecMult (T*, T*): Not tested!!" << STD__ endl;
00681         // No Boundcheck possible!!
00682         unsigned int i;
00683         for (i = 0; i < n_cols; ++i) 
00684                 res[i] = 0;
00685         for (i = 0; i < n_cols; ++i) {
00686                 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
00687                         res[irow[j]] += comp[j]*v[i];  
00688         }
00689 }
00690 
00691 
00692 // transpose-vector multiplication
00693 template <typename T>
00694 TVector<T> CSCMatrix<T>::transMult(const Vector<T>& v) const
00695 {
00696         TVector<T> res(n_cols);
00697         BCHK(n_rows != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
00698         for (unsigned int i = 0; i < n_cols; ++i) {
00699                 T val = 0;
00700                 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) 
00701                         val += comp[j]*v(irow[j]);
00702                 res.set(val, i);
00703         }
00704         return res;
00705 }
00706 
00707 
00708 template <typename T>
00709 TVector<T> CSCMatrix<T>::transMult(const TVector<T>& tv) const
00710 {
00711         Vector<T> v(tv);
00712         return(this->transMult(v));
00713 }
00714 
00715 
00716 template <typename T>
00717 TVector<T> CSCMatrix<T>::transMult(const TSVector<T>& tsv) const
00718 {
00719         Vector<T> v(tsv);
00720         return(this->transMult(v));
00721 }
00722 
00723 
00724 // matrix-scalar multiplication
00725 template <typename T>
00726 CSCMatrix<T>CSCMatrix<T>:: operator* (const T& z)const
00727 {
00728         CSCMatrix<T> res(*this);
00729         for (unsigned int i = 0; i < n_size; ++i) 
00730                 res.comp[i] = comp[i]*z;
00731         return res; // JENS: avoid copy constructor soon!
00732 }
00733 
00734 
00735 // matrix-scalar multiplication
00736 template <typename T>
00737 CSCMatrix<T> CSCMatrix<T>::mult (const T& z) const
00738 {
00739         CSCMatrix<T> res(*this);
00740         for (unsigned int i = 0; i < n_size; ++i) 
00741                 res.comp[i] = z*comp[i];
00742         return res;
00743 }
00744 
00745 
00746 INST(template <typename T> class CSCMatrix friend CSCMatrix<T> operator* (const T&, const CSCMatrix<T>&);)
00747 template <typename T>
00748 inline CSCMatrix<T> operator* (const T& z, const CSCMatrix<T>& m)
00749 { return m.mult (z); }
00750 
00751 // fast in-place matrix-scalar multiplication
00752 template <typename T>
00753 CSCMatrix<T>& CSCMatrix<T>::operator*= (const T& z)
00754 {
00755         for (unsigned int i = 0; i < n_size; ++i) 
00756                 comp[i] *= z;
00757         return *this;
00758 }
00759 
00760 // matrix-scalar division
00761 template <typename T>
00762 CSCMatrix<T> CSCMatrix<T>::operator/ (const T& z) const
00763 {
00764         CSCMatrix<T> res(*this);
00765         for (unsigned int i = 0; i < n_size; ++i) 
00766                 res.comp[i] = comp[i]/z;
00767         return res; // JENS: avoid copy constructor soon!
00768 }
00769 
00770 
00771 // fast in-place matrix-scalar division
00772 template <typename T>
00773 CSCMatrix<T>& CSCMatrix<T>::operator/= (const T& z)
00774 {
00775         for (unsigned int i = 0; i < n_size; ++i) 
00776                 comp[i] /= z;
00777         return *this;
00778 }
00779 
00780 
00781 // allocate memory space
00782 // used by constructor and resize()
00783 template <typename T>
00784 void CSCMatrix<T>::allocate(unsigned int rows, unsigned int columns, 
00785                             unsigned int nnzeros)
00786 {
00787         BCHKNR(rows <= 0 || columns <= 0, CSCMatErr, Zero space allocating, 0);
00788         n_rows = rows;
00789         n_cols = columns;
00790         pcol = new unsigned int[n_cols+1];
00791         BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_cols);
00792         irow = new unsigned int[nnzeros];
00793         BCHKNR(irow==NULL, CSCMatErr, Out of memory, nnzeros);
00794         comp = new T[nnzeros];
00795         BCHKNR(comp==NULL, CSCMatErr, Out of memory, nnzeros);
00796 
00797         unsigned int i;
00798         for (i = 0; i < n_cols+1; ++i)
00799                 pcol[i] = 0;
00800         for (i = 0; i < nnzeros;  ++i)
00801                 irow[i] = 0;
00802 
00803         n_size = 0;
00804         n_max_size = nnzeros;
00805 }
00806 
00807 
00808 // delete memory space
00809 // used by "=" operator and destructor
00810 template <typename T>
00811 void CSCMatrix<T>::destroy ()
00812 {
00813         if (irow != NULL) 
00814                 delete[] irow;
00815         if (comp != NULL) 
00816                 delete[] comp;
00817         if (pcol != NULL) 
00818                 delete[] pcol;
00819 
00820         n_size = 0;
00821         n_max_size = 0;
00822 }
00823 
00824 
00825 // copy to new memory space
00826 // used by "=" operator and copy constructor
00827 template <typename T>
00828 void CSCMatrix<T>::copy (const CSCMatrix<T>& m)
00829 {
00830         n_rows = m.n_rows;
00831         n_cols = m.n_cols;
00832         n_size = m.n_size;
00833         n_max_size = m.n_size;
00834         pcol = new unsigned int[n_cols+1];
00835         BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_cols);
00836         irow = new unsigned int[n_size];
00837         BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_size);
00838         comp = new T[n_size];
00839         BCHKNR(comp==NULL, CSCMatErr, Out of memory, n_size);
00840 
00841         for (unsigned int j = 0; j < n_cols+1; ++j)
00842                 pcol[j] = m.pcol[j];
00843         for (unsigned int i = 0; i < n_size;   ++i) {
00844                 irow[i] = m.irow[i];
00845                 comp[i] = m.comp[i];
00846         }
00847 }
00848 
00849 
00850 
00851 // set all elements to a certain value
00852 template <typename T>
00853 inline CSCMatrix<T>& CSCMatrix<T>::fill (const T& val)
00854 {
00855         if (val == (T)0)
00856                 return this->clear();
00857         //for (unsigned int i = 0; i < n_size; ++i)
00858         //      comp[i] = val;
00859         for (unsigned c = 0; c < columns(); ++c)
00860                 for (unsigned r = 0; r < rows(); ++r)
00861                         (*this)(r, c) = val;
00862         return *this;
00863 }
00864 
00865 // clear all elements
00866 template <typename T>
00867 inline CSCMatrix<T>& CSCMatrix<T>::clear()
00868 { CSTD__ memset (comp, 0, sizeof(T) * n_size); return *this; }
00869 
00870 template <typename T>
00871 CSCMatrix<T>& CSCMatrix<T>::setunit (const T& val)
00872 {
00873         BCHK (n_rows != n_cols, CSCMatErr, setunit is meaningless for non-square matrices, n_cols, *this);
00874         clear ();
00875         for (unsigned int i = 0; i < n_cols; i++)
00876                 setval (val, i, i);
00877         return *this;
00878 }
00879 
00880 
00881 template <typename T>
00882 CSCMatrix<T>& CSCMatrix<T>::swap(CSCMatrix<T>& m)
00883 {       
00884         SWAP(n_rows, m.n_rows); SWAP(n_cols, m.n_cols);
00885         SWAP(n_size, m.n_size); SWAP(n_max_size, m.n_max_size);
00886         SWAP(pcol, m.pcol); SWAP(irow, m.irow);
00887         SWAP(comp, m.comp);
00888         return *this;
00889 }
00890 
00891 template <typename T>
00892 CSCMatrix<T> CSCMatrix<T>::transposed_copy() const
00893 {
00894         CSCMatrix<T> m(columns(), rows());
00895         for (unsigned r = 0; r < rows(); ++r)
00896                 for (unsigned c = 0; c < columns(); ++c)
00897                         m(c,r) = (*this)(r,c);
00898         return m;
00899 }
00900 
00901 template <typename T>
00902 inline CSCMatrix<T>& CSCMatrix<T>::transpose()
00903 {
00904         CSCMatrix<T> tr(this->transposed_copy());
00905         return swap(tr);
00906 }
00907 
00908 INST(template <typename T> class CSCMatrix friend CSCMatrix<T> transpose(const CSCMatrix<T>& cscm);)
00909 template <typename T>
00910 inline CSCMatrix<T> transpose(const CSCMatrix<T>& cscm)
00911 {
00912         return cscm.transposed_copy();
00913 }
00914 
00915 NAMESPACE_END
00916 
00917 #endif  /* TBCI_CSCMATRIX_H */