← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Solver/linearSystemCSR.h

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-27 17:36:40 UTC
  • mfrom: (1.2.9 upstream) (8.1.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090927173640-oxyhzt0eadjfrlwz
[Christophe Prud'homme]
* New upstream release
  + solver code refactoring
  + better IDE integration.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Solver/linearSystemCSR.h
 
1
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
 
2
//
 
3
// See the LICENSE.txt file for license information. Please report all
 
4
// bugs and problems to <gmsh@geuz.org>.
 
5
 
 
6
#ifndef _LINEAR_SYSTEM_CSR_H_
 
7
#define _LINEAR_SYSTEM_CSR_H_
 
8
 
 
9
#include <vector>
 
10
#include "GmshConfig.h"
 
11
#include "GmshMessage.h"
 
12
#include "linearSystem.h"
 
13
 
 
14
typedef int INDEX_TYPE ;  
 
15
typedef struct {
 
16
  int nmax;
 
17
  int size;
 
18
  int incr;
 
19
  int n;
 
20
  int isorder;
 
21
  char *array;
 
22
} CSRList_T;
 
23
  
 
24
void CSRList_Add(CSRList_T *liste, void *data);
 
25
int  CSRList_Nbr(CSRList_T *liste);
 
26
 
 
27
template <class scalar>
 
28
class linearSystemCSR : public linearSystem<scalar> {
 
29
 protected:
 
30
  bool sorted;
 
31
  char *something;
 
32
  CSRList_T *a_,*ai_,*ptr_,*jptr_; 
 
33
  std::vector<scalar> *_b, *_x;
 
34
 public:
 
35
  linearSystemCSR()
 
36
    : sorted(false), a_(0) {}
 
37
  virtual bool isAllocated() const { return a_ != 0; }
 
38
  virtual void allocate(int) ;
 
39
  virtual void clear()
 
40
  {
 
41
    allocate(0);
 
42
  }
 
43
  virtual ~linearSystemCSR()
 
44
  {
 
45
    allocate(0);
 
46
  }
 
47
  virtual void addToMatrix ( int il, int ic, double val) 
 
48
  {
 
49
    //    if (sorted)throw;
 
50
    
 
51
    INDEX_TYPE  *jptr  = (INDEX_TYPE*) jptr_->array;
 
52
    INDEX_TYPE  *ptr   = (INDEX_TYPE*) ptr_->array;
 
53
    INDEX_TYPE  *ai    = (INDEX_TYPE*) ai_->array;
 
54
    scalar      *a     = ( scalar * ) a_->array;
 
55
    
 
56
    INDEX_TYPE  position_ = jptr[il];
 
57
    
 
58
    if(something[il]) {
 
59
      while(1){
 
60
        if(ai[position_] == ic){
 
61
          a[position_] += val;
 
62
          //      if (il == 0)    printf("FOUND %d %d %d\n",il,ic,position_);
 
63
          return;
 
64
        }
 
65
        if (ptr[position_] == 0)break;
 
66
        position_ = ptr[position_];
 
67
      }
 
68
    }  
 
69
    
 
70
    INDEX_TYPE zero = 0;
 
71
    CSRList_Add (a_, &val);
 
72
    CSRList_Add (ai_, &ic);
 
73
    CSRList_Add (ptr_, &zero);
 
74
    // The pointers may have been modified
 
75
    // if there has been a reallocation in CSRList_Add  
 
76
    
 
77
    ptr = (INDEX_TYPE*) ptr_->array;
 
78
    ai  = (INDEX_TYPE*) ai_->array;
 
79
    a   = (scalar*) a_->array;
 
80
    
 
81
    INDEX_TYPE n = CSRList_Nbr(a_) - 1;
 
82
    
 
83
    if(!something[il]) {
 
84
      jptr[il] = n;
 
85
      something[il] = 1;      
 
86
    }
 
87
    else ptr[position_] = n;
 
88
  }
 
89
  
 
90
  virtual scalar getFromMatrix (int _row, int _col) const
 
91
  {
 
92
    throw;
 
93
  }
 
94
  virtual void addToRightHandSide(int _row, scalar _val) 
 
95
  {
 
96
    if(_val != 0.0) (*_b)[_row] += _val;
 
97
  }
 
98
  virtual scalar getFromRightHandSide(int _row) const 
 
99
  {
 
100
    return (*_b)[_row];
 
101
  }
 
102
  virtual scalar getFromSolution(int _row) const
 
103
  {
 
104
    return (*_x)[_row];
 
105
  }
 
106
  virtual void zeroMatrix()
 
107
  {
 
108
    int N=CSRList_Nbr(a_);
 
109
    scalar *a = (scalar*) a_->array;
 
110
    for (int i=0;i<N;i++)a[i]=0;
 
111
  }
 
112
  virtual void zeroRightHandSide() 
 
113
  {
 
114
    for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
 
115
  }
 
116
};
 
117
 
 
118
template <class scalar>
 
119
class linearSystemCSRGmm : public linearSystemCSR<scalar> {
 
120
 private:
 
121
  double _prec;
 
122
  int _noisy, _gmres;
 
123
 public:
 
124
  linearSystemCSRGmm()
 
125
    : _prec(1.e-8), _noisy(0), _gmres(0) {}
 
126
  virtual ~linearSystemCSRGmm()
 
127
  {}
 
128
  void setPrec(double p){ _prec = p; }
 
129
  void setNoisy(int n){ _noisy = n; }
 
130
  void setGmres(int n){ _gmres = n; }
 
131
  virtual int systemSolve() 
 
132
#if defined(HAVE_GMM)
 
133
    ;
 
134
#else
 
135
  {
 
136
    Msg::Error("Gmm++ is not available in this version of Gmsh");
 
137
    return 0;
 
138
  }
 
139
#endif
 
140
};
 
141
 
 
142
#endif