← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Solver/linearSystemCSR.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-27 17:36:40 UTC
  • mfrom: (1.4.2 upstream)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20090927173640-meoeywl4f5hq5qas
Tags: 2.4.2.dfsg-1
[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.cpp
 
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
#include <stdlib.h>
 
7
#include <stdio.h>
 
8
#include <string.h>
 
9
#include "GmshConfig.h"
 
10
#include "GmshMessage.h"
 
11
#include "linearSystemCSR.h"
 
12
 
 
13
#define SWAP(a,b)  temp=(a);(a)=(b);(b)=temp;
 
14
#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
 
15
 
 
16
static void *CSRMalloc(size_t size)
 
17
{
 
18
  void *ptr;
 
19
  if (!size) return(NULL);
 
20
  ptr = malloc(size);
 
21
  return(ptr);
 
22
}
 
23
 
 
24
static void *CSRRealloc(void *ptr, size_t size)
 
25
{
 
26
  if (!size) return(NULL);
 
27
  ptr = realloc(ptr,size);
 
28
  return(ptr);
 
29
}
 
30
 
 
31
static void CSRList_Realloc(CSRList_T *liste,int n)
 
32
{
 
33
  char* temp;
 
34
  if (n <= 0) return;
 
35
  if (liste->array == NULL) {
 
36
    liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
 
37
    liste->array = (char *)CSRMalloc(liste->nmax * liste->size);
 
38
  }
 
39
  else {
 
40
    if (n > liste->nmax) {
 
41
      liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
 
42
      temp = (char *)CSRRealloc(liste->array, liste->nmax * liste->size);
 
43
      liste->array = temp;
 
44
    }
 
45
  }
 
46
}
 
47
 
 
48
 
 
49
static CSRList_T *CSRList_Create(int n, int incr, int size)
 
50
{
 
51
  CSRList_T *liste;
 
52
  
 
53
  if (n <= 0)  n = 1 ;
 
54
  if (incr <= 0) incr = 1;
 
55
  
 
56
  liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
 
57
  
 
58
  liste->nmax    = 0;
 
59
  liste->incr    = incr;
 
60
  liste->size    = size;
 
61
  liste->n       = 0;
 
62
  liste->isorder = 0;
 
63
  liste->array   = NULL;
 
64
  
 
65
  CSRList_Realloc(liste,n);
 
66
  return(liste);
 
67
}
 
68
 
 
69
static void CSRList_Delete(CSRList_T *liste)
 
70
{
 
71
  if (liste != 0) {
 
72
    free(liste->array);
 
73
    free(liste);
 
74
  }
 
75
}
 
76
 
 
77
void CSRList_Add(CSRList_T *liste, void *data)
 
78
{
 
79
  liste->n++;
 
80
  
 
81
  CSRList_Realloc(liste,liste->n);
 
82
  liste->isorder = 0;
 
83
  memcpy(&liste->array[(liste->n - 1) * liste->size],data,liste->size);
 
84
}
 
85
 
 
86
int CSRList_Nbr(CSRList_T *liste)
 
87
{
 
88
  return(liste->n);
 
89
}
 
90
 
 
91
template<>
 
92
void linearSystemCSR<double>::allocate(int _nbRows)
 
93
{
 
94
  if(a_) {
 
95
    CSRList_Delete(a_);
 
96
    CSRList_Delete(ai_);
 
97
    CSRList_Delete(ptr_);
 
98
    CSRList_Delete(jptr_);
 
99
    delete _x;
 
100
    delete _b;
 
101
    delete something;
 
102
  }
 
103
  
 
104
  if(_nbRows == 0){
 
105
    a_ = 0; 
 
106
    ai_ = 0; 
 
107
    ptr_ = 0; 
 
108
    jptr_ = 0; 
 
109
    _b = 0;
 
110
    _x = 0;
 
111
    something = 0;
 
112
    return;
 
113
  }
 
114
  
 
115
  a_    = CSRList_Create (_nbRows, _nbRows, sizeof(double));
 
116
  ai_   = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
 
117
  ptr_  = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
 
118
  jptr_ = CSRList_Create (_nbRows+1, _nbRows, sizeof(INDEX_TYPE));
 
119
  
 
120
  something = new char[_nbRows];
 
121
  
 
122
  for (int i=0;i<_nbRows;i++)something[i] = 0;
 
123
  
 
124
  _b = new std::vector<double>(_nbRows);
 
125
  _x = new std::vector<double>(_nbRows);
 
126
}
 
127
 
 
128
const int NSTACK   = 50;
 
129
const unsigned int M_sort2  = 7;
 
130
 
 
131
static void free_ivector(int *v, long nl, long nh)
 
132
{
 
133
  // free an int vector allocated with ivector() 
 
134
  free((char*)(v+nl-1));
 
135
}
 
136
 
 
137
static int *ivector(long nl, long nh)
 
138
{
 
139
  // allocate an int vector with subscript range v[nl..nh] 
 
140
  int *v;  
 
141
  v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
 
142
  if (!v) fprintf(stderr, "allocation failure in ivector()\n");
 
143
  return v-nl+1;
 
144
}
 
145
 
 
146
static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj)
 
147
{
 
148
  if(ai<bi)return -1;
 
149
  if(ai>bi)return 1;
 
150
  if(aj<bj)return -1;
 
151
  if(aj>bj)return 1;
 
152
  return 0;
 
153
}
 
154
 
 
155
template <class scalar>
 
156
static void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
 
157
{
 
158
  unsigned long i,ir=n,j,k,l=1;
 
159
  int *istack,jstack=0;
 
160
  INDEX_TYPE tempi;
 
161
  scalar a,temp;
 
162
  int    b,c;
 
163
  
 
164
  istack=ivector(1,NSTACK);
 
165
  for (;;) {
 
166
    if (ir-l < M_sort2) {
 
167
      for (j=l+1;j<=ir;j++) {
 
168
        a=arr[j -1];
 
169
        b=ai[j -1];
 
170
        c=aj[j -1];
 
171
        for (i=j-1;i>=1;i--) {
 
172
          if (cmpij(ai[i -1],aj[i -1],b,c) <= 0) break;
 
173
          arr[i+1 -1]=arr[i -1];
 
174
          ai[i+1 -1]=ai[i -1];
 
175
          aj[i+1 -1]=aj[i -1];
 
176
        }
 
177
        arr[i+1 -1]=a;
 
178
        ai[i+1 -1]=b;
 
179
        aj[i+1 -1]=c;
 
180
      }
 
181
      if (!jstack) {
 
182
        free_ivector(istack,1,NSTACK);
 
183
        return;
 
184
      }
 
185
      ir=istack[jstack];
 
186
      l=istack[jstack-1];
 
187
      jstack -= 2;
 
188
    } 
 
189
    else {
 
190
      k=(l+ir) >> 1;
 
191
      SWAP(arr[k -1],arr[l+1 -1])
 
192
        SWAPI(ai[k -1],ai[l+1 -1])
 
193
        SWAPI(aj[k -1],aj[l+1 -1])
 
194
        if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[ir -1],aj[ir -1])>0){
 
195
          SWAP(arr[l+1 -1],arr[ir -1])
 
196
            SWAPI(ai[l+1 -1],ai[ir -1])
 
197
            SWAPI(aj[l+1 -1],aj[ir -1])
 
198
            }
 
199
      if (cmpij(ai[l -1],aj[l -1],ai[ir -1],aj[ir -1])>0){
 
200
        SWAP(arr[l -1],arr[ir -1])
 
201
          SWAPI(ai[l -1],ai[ir -1])
 
202
          SWAPI(aj[l -1],aj[ir -1])
 
203
          }
 
204
      if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[l -1],aj[l -1])>0){
 
205
        SWAP(arr[l+1 -1],arr[l -1])
 
206
          SWAPI(ai[l+1 -1],ai[l -1])
 
207
          SWAPI(aj[l+1 -1],aj[l -1])
 
208
          }
 
209
      i=l+1;
 
210
      j=ir;
 
211
      a=arr[l -1];
 
212
      b=ai[l -1];
 
213
      c=aj[l -1];
 
214
      for (;;) {
 
215
        do i++; while (cmpij(ai[i -1],aj[i -1],b,c) < 0);
 
216
        do j--; while (cmpij(ai[j -1],aj[j -1],b,c) > 0);
 
217
        if (j < i) break;
 
218
        SWAP(arr[i -1],arr[j -1])
 
219
          SWAPI(ai[i -1],ai[j -1])
 
220
          SWAPI(aj[i -1],aj[j -1])
 
221
          }
 
222
      arr[l -1]=arr[j -1];
 
223
      arr[j -1]=a;
 
224
      ai[l -1]=ai[j -1];
 
225
      ai[j -1]=b;
 
226
      aj[l -1]=aj[j -1];
 
227
      aj[j -1]=c;
 
228
      jstack += 2;
 
229
      if (jstack > NSTACK) {
 
230
        Msg::Fatal("NSTACK too small while sorting the columns of the matrix");
 
231
        throw;
 
232
      }
 
233
      if (ir-i+1 >= j-l) {
 
234
        istack[jstack]=ir;
 
235
        istack[jstack-1]=i;
 
236
        ir=j-1;
 
237
      } 
 
238
      else {
 
239
        istack[jstack]=j-1;
 
240
        istack[jstack-1]=l;
 
241
        l=i;
 
242
      }
 
243
    }
 
244
  }
 
245
}
 
246
 
 
247
template <class scalar>
 
248
static void sortColumns(int NbLines, 
 
249
                        int nnz, 
 
250
                        INDEX_TYPE *ptr, 
 
251
                        INDEX_TYPE *jptr, 
 
252
                        INDEX_TYPE *ai, 
 
253
                        scalar *a)
 
254
{
 
255
  // replace pointers by lines
 
256
  int *count = new int [NbLines];
 
257
  
 
258
  for(int i=0;i<NbLines;i++){
 
259
    count[i] = 0;
 
260
    INDEX_TYPE _position =  jptr[i];
 
261
    while(1){
 
262
      count[i]++;
 
263
      INDEX_TYPE _position_temp = _position;
 
264
      _position = ptr[_position];
 
265
      ptr[_position_temp] = i;
 
266
      if (_position == 0) break;
 
267
    }
 
268
  }   
 
269
  _sort2_xkws<double>(nnz,a,ptr,ai);   
 
270
  jptr[0] = 0;
 
271
  for(int i=1;i<=NbLines;i++){
 
272
    jptr[i] = jptr[i-1] + count[i-1];
 
273
  }
 
274
  
 
275
  for(int i=0;i<NbLines;i++){
 
276
    for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
 
277
      ptr[j] = j+1;
 
278
    }
 
279
    ptr[jptr[i+1]] = 0;
 
280
  }
 
281
  
 
282
  delete[] count;
 
283
}
 
284
 
 
285
#if defined(HAVE_GMM)
 
286
 
 
287
#include "gmm.h"
 
288
 
 
289
template<>
 
290
int linearSystemCSRGmm<double>::systemSolve()
 
291
{
 
292
  if (!sorted)
 
293
    sortColumns(_b->size(),
 
294
                CSRList_Nbr(a_),
 
295
                (INDEX_TYPE *) ptr_->array,
 
296
                (INDEX_TYPE *) jptr_->array, 
 
297
                (INDEX_TYPE *) ai_->array, 
 
298
                (double*) a_->array);
 
299
  sorted = true;
 
300
  
 
301
  gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>  
 
302
    ref((double*)a_->array, (INDEX_TYPE *) ai_->array,
 
303
        (INDEX_TYPE *)jptr_->array, _b->size(), _b->size());
 
304
  gmm::csr_matrix<double,0> M;
 
305
  M.init_with(ref);
 
306
  
 
307
  gmm::ildltt_precond<gmm::csr_matrix<double, 0> > P(M, 10, 1.e-10);
 
308
  gmm::iteration iter(_prec);
 
309
  iter.set_noisy(_noisy);
 
310
  if(_gmres) gmm::gmres(M, *_x, *_b, P, 100, iter);
 
311
  else gmm::cg(M, *_x, *_b, P, iter);
 
312
  return 1;
 
313
}
 
314
 
 
315
#endif
 
316