SparseLibrary  Version 1.6.0
HBICRS.hpp
1 /*
2  * Copyright (c) 2007-2014, A. N. Yzelman, Utrecht University 2007-2011;
3  * KU Leuven 2011-2014.
4  * R. H. Bisseling, Utrecht University 2007-2014.
5  *
6  * This file is part of the Sparse Library.
7  *
8  * This library was developed under supervision of Prof. dr. Rob H. Bisseling at
9  * Utrecht University, from 2007 until 2011. From 2011-2014, development continued
10  * at KU Leuven, where Prof. dr. Dirk Roose contributed significantly to the ideas
11  * behind the newer parts of the library code.
12  *
13  * The Sparse Library is free software: you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by the
15  * Free Software Foundation, either version 3 of the License, or (at your
16  * option) any later version.
17  *
18  * The Sparse Library is distributed in the hope that it will be useful, but
19  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
20  * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  * for more details.
22  *
23  * You should have received a copy of the GNU General Public License along
24  * with the Sparse Library. If not, see <http://www.gnu.org/licenses/>.
25  */
26 
27 
28 /*
29  * File created by:
30  * A. N. Yzelman, Dept. of Mathematics, Utrecht University, 2010.
31  */
32 
33 
34 #include "Matrix.hpp"
35 #include "SparseMatrix.hpp"
36 #include "TS.hpp"
37 #include "CRS.hpp"
38 #include "ICRS.hpp"
39 #include "ZZ_CRS.hpp"
40 #include "ZZ_ICRS.hpp"
41 #include "SVM.hpp"
42 #include "HTS.hpp"
43 #include "BICRS.hpp"
44 #include "CCSWrapper.hpp"
45 
46 
47 #include <assert.h>
48 #include <iostream>
49 
50 // Set when in debug mode
51 //#define _DEBUG
52 
53 #ifndef _H_HBICRS
54 #define _H_HBICRS
55 
77 template< typename _t_value >
78 class HBICRS: public SparseMatrix< _t_value, LI > {
79 
80  protected:
81 
83  size_t jumps;
84 
86  LI* r_inc;
87 
89  LI* c_inc;
90 
92  Matrix< _t_value > **vals;
93 
95  ULI ntt;
96 
97  public:
98 
100  ~HBICRS() {
101  delete [] r_inc;
102  delete [] c_inc;
103  for( LI i=0; i<this->nnz; i++ )
104  delete vals[ i ];
105  delete [] vals;
106  }
107 
109  HBICRS() {}
110 
115  HBICRS( std::string file, _t_value zero = 0 ) {
116  std::cerr << "Error: Hierarchical BICRS, loading in from MatrixMarket file not possible, since no way to infer a hierarchy!" << std::endl;
117  exit( 1 );
118  }
119 
123  HBICRS( std::vector< Triplet< _t_value > >& input, LI m, LI n, _t_value zero ) {
124  load( input, m, n, zero );
125  }
126 
130  HBICRS( std::vector< std::vector< Triplet< _t_value > > >& input, signed char* group_type, LI m, LI n, _t_value zero = 0 ) {
131  load( input, group_type, m, n, zero );
132  }
133 
140  virtual void load( std::vector< Triplet< _t_value > >& input, LI m, LI n, _t_value zero ) {
141  std::cerr << "Error: Hierarchical BICRS, loading in from plain triplets not possible, since no way to infer a hierarchy!" << std::endl;
142  exit( 1 );
143  }
144 
146  virtual void load( std::vector< std::vector< Triplet< _t_value > > >& input, signed char* group_type, LI m, LI n, _t_value zero ) {
147  LI* rrow = new LI[ input.size() ];
148  LI* ccol = new LI[ input.size() ];
149 
150  if( group_type == NULL ) {
151  char *group_type = new char[ input.size() ];
152  for( unsigned int i=0; i<input.size(); i++ ) group_type[i] = 2;
153  }
154 
155  vals = new Matrix< _t_value >*[ input.size() ];
156 
157  unsigned long int g = 0;
158  typename std::vector< std::vector< Triplet< _t_value > > >::iterator it = input.begin();
159  for( ; it!=input.end(); it++, g++ ) {
160  //skip empty blocks
161  if( it->size() == 0 ) {
162  g--;
163  continue;
164  }
165  //std::cout << "\tHBICRS: Adding a group type " << static_cast< int >( group_type[g] ) << std::endl;
166  switch( group_type[ g ] ) {
167  case -7:
168  vals[g] = new CCSWrapper< double, BICRS< _t_value >, ULI >( *it, m, n, zero );
169  rrow[g] = 0;
170  ccol[g] = 0;
171  break;
172  case -6:
173  vals[g] = new CCSWrapper< double, HTS< _t_value >, ULI >( *it, m, n, zero );
174  rrow[g] = 0;
175  ccol[g] = 0;
176  break;
177  case -5:
178  vals[g] = new CCSWrapper< double, SVM< _t_value >, ULI >( *it, m, n, zero );
179  rrow[g] = 0;
180  ccol[g] = 0;
181  break;
182  case -4:
183  vals[g] = new CCSWrapper< _t_value, ZZ_ICRS< _t_value >, LI >( *it, m, n, zero );
184  rrow[g] = 0;
185  ccol[g] = 0;
186  break;
187  case -3:
188  vals[g] = new CCSWrapper< double, ZZ_CRS< _t_value >, LI >( *it, m, n, zero );
189  rrow[g] = 0;
190  ccol[g] = 0;
191  break;
192  case -2:
193  vals[g] = new CCSWrapper< double, ICRS< _t_value >, ULI >( *it, m, n, zero );
194  rrow[g] = 0;
195  ccol[g] = 0;
196  break;
197  case -1:
198  vals[g] = new CCSWrapper< double, CRS< _t_value >, ULI >( *it, m, n, zero );
199  rrow[g] = 0;
200  ccol[g] = 0;
201  break;
202  case 0:
203  vals[g] = new TS< _t_value >( *it, m, n, zero );
204  rrow[g] = 0;
205  ccol[g] = 0;
206  break;
207  case 1:
208  vals[g] = new CRS< _t_value >( *it, m, n, zero );
209  rrow[g] = 0;
210  ccol[g] = 0;
211  break;
212  case 2:
213  vals[g] = new ICRS< _t_value >( *it, m, n, zero );
214  rrow[g] = 0;
215  ccol[g] = 0;
216  break;
217  case 3:
218  vals[g] = new ZZ_CRS< _t_value >( *it, m, n, zero );
219  rrow[g] = 0;
220  ccol[g] = 0;
221  break;
222  case 4:
223  vals[g] = new ZZ_ICRS< _t_value >( *it, m, n, zero );
224  rrow[g] = 0;
225  ccol[g] = 0;
226  break;
227  case 5:
228  vals[g] = new SVM< _t_value >( *it, m, n, zero );
229  rrow[g] = 0;
230  ccol[g] = 0;
231  break;
232  case 6:
233  vals[g] = new HTS< _t_value >( *it, m, n, zero );
234  rrow[g] = 0;
235  ccol[g] = 0;
236  break;
237  case 7:
238  vals[g] = new BICRS< _t_value >( *it, m, n, zero );
239  rrow[g] = 0;
240  ccol[g] = 0;
241  break;
242  default:
243  std::cerr << "Hierarchical sparse matrix: invalid subscheme ID (" << static_cast<int>(group_type[g]) << "), group " << g << "ignored!" << std::endl;
244  }
245 
246  }
247  this->zero_element = zero;
248  load( rrow, ccol, m, n, g );
249  delete [] rrow;
250  delete [] ccol;
251  }
252 
254  void load( LI* row, LI* col, LI m, LI n, LI nb ) {
255 #ifdef _DEBUG
256  std::cerr << "Warning: _DEBUG flag set." << std::endl;
257 #endif
258 
259  this->nnz = nb;
260  this->nor = m;
261  this->noc = n;
262  this->ntt=2*n;
263  jumps = 0;
264  int prevrow = row[ 0 ];
265  for( LI i=1; i<this->nnz; i++ ) {
266  if( row[ i ] != prevrow )
267  jumps++;
268  prevrow = row[ i ];
269  }
270 #ifdef _DEBUG
271  std::cout << jumps << " row jumps found." << std::endl;
272 #endif
273  r_inc = new LI[ jumps + 2 ];
274  c_inc = new LI[ this->nnz + 1 ];
275 
276  r_inc[ 0 ] = prevrow = row[ 0 ];
277  int prevcol = c_inc[ 0 ] = col[ 0 ];
278 #ifdef _DEBUG
279  std::cout << "c_inc: " << prevcol << std::endl;
280  std::cout << "r_inc: " << prevrow << std::endl;
281 #endif
282  int c = 1;
283  for( LI i=1; i<this->nnz; i++ ) {
284  this->c_inc[ i ] = col[ i ] - prevcol;
285  if( row[ i ] != prevrow ) {
286  this->c_inc[ i ] += ntt;
287  this->r_inc[ c++ ] = row[ i ] - prevrow;
288 #ifdef _DEBUG
289  std::cout << "c_inc: " << ntt << std::endl;
290  std::cout << "r_inc: " << row[ i ] - prevrow << std::endl;
291 #endif
292  prevrow = row[ i ];
293  }
294 #ifdef _DEBUG
295  else
296  std::cout << "c_inc: " << col[ i ] - prevcol << std::endl;
297 #endif
298  prevcol = col[ i ];
299  }
300  //overflow so to signal end of matrix
301  c_inc[ this->nnz ] = ntt;
302  //initialise last row jump to zero (prevent undefined jump)
303  r_inc[ c ] = 0;
304 #ifdef _DEBUG
305  std::cout << "Construction done." << std::endl;
306 #endif
307  }
308 
310  virtual void getFirstIndexPair( LI &i, LI &j ) {
311  i = r_inc[ 0 ];
312  j = c_inc[ 0 ];
313  }
314 
320  virtual void zxa( const _t_value*__restrict__ x_p, _t_value*__restrict__ y_p ) {
321  const _t_value * y = y_p;
322  const _t_value * y_end = y+this->noc;
323  LI *__restrict__ c_inc_p = c_inc;
324  LI *__restrict__ r_inc_p = r_inc;
325 
326 #ifndef NDEBUG
327  const _t_value * x = x_p;
328  const _t_value *__restrict__ x_end = x+this->nor;
329  const LI *__restrict__ c_inc_end = c_inc+this->nnz+1;
330 #endif
331 
332  Matrix< _t_value >**__restrict__ v_end = vals+this->nnz;
333  Matrix< _t_value >**__restrict__ v_p = vals;
334 #ifdef _DEBUG
335  //simulation trace:
336  int yc = c_inc[0];
337  int xc = r_inc[0];
338  std::cout << "( " << xc << " , " << yc << " )" << std::endl;
339  int t=1;
340  for( int i=1; i<this->nnz; i++ ) {
341  yc += c_inc[i];
342  if( yc > this->noc ) {
343  yc -= ntt;
344  xc += r_inc[t++];
345  }
346  assert( yc < this->noc );
347  assert( xc < this->nor );
348  assert( yc >=0 );
349  assert( xc >=0 );
350  std::cout << "( " << yc << " , " << xc << " )" << std::endl;
351  }
352 #endif
353  x_p += *r_inc_p++;
354  y_p += *c_inc_p++;
355  while( v_p < v_end ) {
356  assert( y_p >= y );
357  assert( y_p < y_end );
358  assert( v_p >= vals );
359  assert( v_p < v_end );
360  assert( x_p >= x );
361  assert( x_p < x_end );
362  assert( c_inc_p >= c_inc );
363  assert( c_inc_p < c_inc_end );
364  assert( r_inc_p >= r_inc );
365  while( y_p < y_end ) {
366  (*(v_p++))->zxa( x_p, y_p );
367  y_p += *c_inc_p++;
368  }
369  y_p -= ntt;
370  x_p += *r_inc_p++;
371  }
372  }
373 
379  virtual void zax( const _t_value*__restrict__ x_p, _t_value*__restrict__ y_p ) {
380  const _t_value * x_end = x_p+this->noc;
381  LI *__restrict__ c_inc_p = c_inc;
382  LI *__restrict__ r_inc_p = r_inc;
383 #ifndef NDEBUG
384  const _t_value * x = x_p;
385  const _t_value * const y = y_p;
386  const _t_value * y_end = y+this->nor;
387  const LI *__restrict__ c_inc_end = c_inc+this->nnz+1;
388 #endif
389 
390  Matrix< _t_value >**__restrict__ v_end = vals+this->nnz;
391  Matrix< _t_value >**__restrict__ v_p = vals;
392 
393 #ifdef _DEBUG
394  //simulation trace:
395  int xc = c_inc[0];
396  int yc = r_inc[0];
397  std::cout << "( " << xc << " , " << yc << " )" << std::endl;
398  int t=1;
399  for( int i=1; i<this->nnz; i++ ) {
400  xc += c_inc[i];
401  if( xc > this->noc ) {
402  xc -= ntt;
403  yc += r_inc[t++];
404  }
405  assert( xc < this->nor );
406  assert( yc < this->noc );
407  assert( xc >=0 );
408  assert( yc >=0 );
409  std::cout << "( " << xc << " , " << yc << " )" << std::endl;
410  }
411 #endif
412 
413  x_p += *c_inc_p++;
414  y_p += *r_inc_p++;
415  while( v_p < v_end ) {
416  assert( y_p >= y );
417  assert( y_p < y_end );
418  assert( v_p >= vals );
419  assert( v_p < v_end );
420  assert( x_p >= x );
421  assert( x_p < x_end );
422  assert( c_inc_p >= c_inc );
423  assert( c_inc_p < c_inc_end );
424  assert( r_inc_p >= r_inc );
425  while( x_p < x_end ) {
426  (*(v_p++))->zax( x_p, y_p );
427  x_p += *c_inc_p++;
428  }
429  x_p -= ntt;
430  y_p += *r_inc_p++;
431  }
432  }
433 
434  virtual size_t bytesUsed() {
435  //base case; indices + pointers to substructures
436  size_t ret = sizeof( LI ) * ( this->nnz + jumps + 2 ) + sizeof( void* ) * this->nnz;
437  //get leaf data-structure memory use
438  for( size_t i = 0; i < static_cast< size_t >(this->nnz); ++i )
439  ret += vals[ i ]->bytesUsed();
440  return ret;
441  }
442 
443 };
444 
445 #endif
446 
SparseMatrix< _t_value, LI >::nnz
LI nnz
Number of non-zeros.
Definition: SparseMatrix.hpp:58
Matrix.hpp
File created by: A.
CCSWrapper
Automatically transforms a row-major scheme into an column-major scheme.
Definition: CCSWrapper.hpp:49
ICRS
The incremental compressed row storage sparse matrix data structure.
Definition: ICRS.hpp:53
BICRS
Bi-directional Incremental Compressed Row Storage scheme.
Definition: BICRS.hpp:58
HBICRS::~HBICRS
~HBICRS()
Base deconstructor.
Definition: HBICRS.hpp:100
HBICRS
Hierarchical Bi-directional Incremental Compressed Row Storage scheme.
Definition: HBICRS.hpp:78
SparseMatrix< _t_value, LI >::m
virtual unsigned long int m()
Queries the number of rows this matrix contains.
Definition: SparseMatrix.hpp:107
HBICRS::zax
virtual void zax(const _t_value *__restrict__ x_p, _t_value *__restrict__ y_p)
Calculates y=Ax, but does not allocate y itself.
Definition: HBICRS.hpp:379
HBICRS::HBICRS
HBICRS(std::vector< std::vector< Triplet< _t_value > > > &input, signed char *group_type, LI m, LI n, _t_value zero=0)
Base constructor.
Definition: HBICRS.hpp:130
CRS
The compressed row storage sparse matrix data structure.
Definition: CRS.hpp:52
HBICRS::load
virtual void load(std::vector< Triplet< _t_value > > &input, LI m, LI n, _t_value zero)
This function will rewrite the std::vector< Triplet > structure to one suitable for the other load fu...
Definition: HBICRS.hpp:140
TS
The triplet scheme; a storage scheme for sparse matrices using triplets.
Definition: TS.hpp:44
HBICRS::load
virtual void load(std::vector< std::vector< Triplet< _t_value > > > &input, signed char *group_type, LI m, LI n, _t_value zero)
Constructs the hierarchical part.
Definition: HBICRS.hpp:146
HBICRS::HBICRS
HBICRS(std::vector< Triplet< _t_value > > &input, LI m, LI n, _t_value zero)
Base constructor.
Definition: HBICRS.hpp:123
HBICRS::vals
Matrix< _t_value > ** vals
Stores the values individual storage schemes.
Definition: HBICRS.hpp:92
SparseMatrix
Interface common to all sparse matrix storage schemes.
Definition: SparseMatrix.hpp:46
ICRS.hpp
File created by: A.
ZZ_ICRS
The zig-zag incremental compressed row storage sparse matrix data structure.
Definition: ZZ_ICRS.hpp:58
SparseMatrix< _t_value, LI >::noc
LI noc
Number of columns.
Definition: SparseMatrix.hpp:55
HBICRS::r_inc
LI * r_inc
Stores the row jumps; size is at maximum the number of nonzeros.
Definition: HBICRS.hpp:86
SVM
The sparse vector matrix format.
Definition: SVM.hpp:51
SparseMatrix< _t_value, LI >::nor
LI nor
Number of rows.
Definition: SparseMatrix.hpp:52
HBICRS::zxa
virtual void zxa(const _t_value *__restrict__ x_p, _t_value *__restrict__ y_p)
Calculates y=xA, but does not allocate y itself.
Definition: HBICRS.hpp:320
SparseMatrix< _t_value, LI >::zero_element
_t_value zero_element
The element considered to be zero.
Definition: SparseMatrix.hpp:63
ZZ_CRS
The zig-zag compressed row storage sparse matrix data structure.
Definition: ZZ_CRS.hpp:50
HTS
The Hilbert triplet scheme.
Definition: HTS.hpp:51
HBICRS::HBICRS
HBICRS()
Base constructor.
Definition: HBICRS.hpp:109
HBICRS::jumps
size_t jumps
The number of row jumps.
Definition: HBICRS.hpp:83
SparseMatrix< _t_value, LI >::n
virtual unsigned long int n()
Queries the number of columns this matrix contains.
Definition: SparseMatrix.hpp:115
Triplet
A single triplet value.
Definition: Triplet.hpp:52
HBICRS::load
void load(LI *row, LI *col, LI m, LI n, LI nb)
Builds the BICRS structure.
Definition: HBICRS.hpp:254
HBICRS::HBICRS
HBICRS(std::string file, _t_value zero=0)
Base constructor.
Definition: HBICRS.hpp:115
HBICRS::bytesUsed
virtual size_t bytesUsed()
Function to query the amount of storage required by this sparse matrix.
Definition: HBICRS.hpp:434
HBICRS::ntt
ULI ntt
Caches n times two.
Definition: HBICRS.hpp:95
HBICRS::c_inc
LI * c_inc
Stores the column jumps; size is exactly the number of nonzeros.
Definition: HBICRS.hpp:89
HBICRS::getFirstIndexPair
virtual void getFirstIndexPair(LI &i, LI &j)
Definition: HBICRS.hpp:310
Generated on Wed Nov 19 2014 14:47:33 for SparseLibrary by   doxygen 1.8.7