← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Numeric/GmshMatrix.cpp

  • 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:
Numeric/GmshMatrix.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 <complex>
7
 
#include <string.h>
8
 
#include "GmshConfig.h"
9
 
#include "GmshMatrix.h"
10
 
#include "GmshMessage.h"
11
 
 
12
 
//#if defined(_MSC_VER)
13
 
//#define F77NAME(x) (x)
14
 
//#endif
15
 
 
16
 
#if !defined(F77NAME)
17
 
#define F77NAME(x) (x##_)
18
 
#endif
19
 
 
20
 
#if defined(HAVE_BLAS)
21
 
 
22
 
extern "C" {
23
 
  void F77NAME(dgemm)(const char *transa, const char *transb, int *m, int *n, int *k, 
24
 
                      double *alpha, double *a, int *lda, 
25
 
                      double *b, int *ldb, double *beta, 
26
 
                      double *c, int *ldc);
27
 
  void F77NAME(zgemm)(const char *transa, const char *transb, int *m, int *n, int *k, 
28
 
                      std::complex<double> *alpha, std::complex<double> *a, int *lda, 
29
 
                      std::complex<double> *b, int *ldb, std::complex<double> *beta, 
30
 
                      std::complex<double> *c, int *ldc);
31
 
  void F77NAME(dgemv)(const char *trans, int *m, int *n, 
32
 
                      double *alpha, double *a, int *lda, 
33
 
                      double *x, int *incx, double *beta, 
34
 
                      double *y, int *incy);
35
 
  void F77NAME(zgemv)(const char *trans, int *m, int *n, 
36
 
                      std::complex<double> *alpha, std::complex<double> *a, int *lda, 
37
 
                      std::complex<double> *x, int *incx, std::complex<double> *beta, 
38
 
                      std::complex<double> *y, int *incy);
39
 
}
40
 
 
41
 
template<> 
42
 
void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c)
43
 
{
44
 
  int M = c.size1(), N = c.size2(), K = _c;
45
 
  int LDA = _r, LDB = b.size1(), LDC = c.size1();
46
 
  double alpha = 1., beta = 0.;
47
 
  F77NAME(dgemm)("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB, 
48
 
                 &beta, c._data, &LDC);
49
 
}
50
 
 
51
 
template<> 
52
 
void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<double> > &b, 
53
 
                                             gmshMatrix<std::complex<double> > &c)
54
 
{
55
 
  int M = c.size1(), N = c.size2(), K = _c;
56
 
  int LDA = _r, LDB = b.size1(), LDC = c.size1();
57
 
  std::complex<double> alpha = 1., beta = 0.;
58
 
  F77NAME(zgemm)("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB, 
59
 
                 &beta, c._data, &LDC);
60
 
}
61
 
 
62
 
template<> 
63
 
void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b, 
64
 
                              double alpha, double beta)
65
 
{
66
 
  int M = size1(), N = size2(), K = a.size2();
67
 
  int LDA = a.size1(), LDB = b.size1(), LDC = size1();
68
 
  F77NAME(dgemm)("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB, 
69
 
                 &beta, _data, &LDC);
70
 
}
71
 
 
72
 
template<> 
73
 
void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &a, 
74
 
                                             gmshMatrix<std::complex<double> > &b, 
75
 
                                             std::complex<double> alpha, 
76
 
                                             std::complex<double> beta)
77
 
{
78
 
  int M = size1(), N = size2(), K = a.size2();
79
 
  int LDA = a.size1(), LDB = b.size1(), LDC = size1();
80
 
  F77NAME(zgemm)("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB, 
81
 
                 &beta, _data, &LDC);
82
 
}
83
 
 
84
 
template<> 
85
 
void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y)
86
 
{
87
 
  int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
88
 
  double alpha = 1., beta = 0.;
89
 
  F77NAME(dgemv)("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
90
 
                 &beta, y._data, &INCY);
91
 
}
92
 
 
93
 
template<> 
94
 
void gmshMatrix<std::complex<double> >::mult(const gmshVector<std::complex<double> > &x, 
95
 
                                             gmshVector<std::complex<double> > &y)
96
 
{
97
 
  int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
98
 
  std::complex<double> alpha = 1., beta = 0.;
99
 
  F77NAME(zgemv)("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
100
 
                 &beta, y._data, &INCY);
101
 
}
102
 
 
103
 
#endif
104
 
 
105
 
#if defined(HAVE_LAPACK)
106
 
 
107
 
extern "C" {
108
 
  void F77NAME(dgesv)(int *N, int *nrhs, double *A, int *lda, int *ipiv, 
109
 
                      double *b, int *ldb, int *info);
110
 
  void F77NAME(dgetrf)(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
111
 
  void F77NAME(dgetri)(int *M, double *A, int *lda, int *ipiv, double *work, int *lwork, 
112
 
                       int *info);
113
 
  void F77NAME(dgesvd)(const char* jobu, const char *jobvt, int *M, int *N,
114
 
                       double *A, int *lda, double *S, double* U, int *ldu,
115
 
                       double *VT, int *ldvt, double *work, int *lwork, int *info);
116
 
  void F77NAME(dgeev)(const char *jobvl, const char *jobvr, 
117
 
                      int *n, double *a, int *lda, 
118
 
                      double *wr, double *wi, 
119
 
                      double *vl, int *ldvl, 
120
 
                      double *vr, int *ldvr, 
121
 
                      double *work, int *lwork,
122
 
                      int *info); 
123
 
}
124
 
 
125
 
template<> 
126
 
bool gmshMatrix<double>::invertInPlace()
127
 
{
128
 
  int N = size1(), nrhs = N, lda = N, ldb = N, info;
129
 
  int *ipiv = new int[N];
130
 
  double * invA = new double[N*N];
131
 
 
132
 
  for (int i = 0; i < N * N; i++) invA[i] = 0.;
133
 
  for (int i = 0; i < N; i++) invA[i * N + i] = 1.;
134
 
 
135
 
  F77NAME(dgesv)(&N, &nrhs, _data, &lda, ipiv, invA, &ldb, &info);
136
 
  memcpy(_data, invA, N * N * sizeof(double));
137
 
 
138
 
  delete [] invA;
139
 
  delete [] ipiv;
140
 
 
141
 
  if(info == 0) return true;
142
 
  if(info > 0)
143
 
    Msg::Error("U(%d,%d)=0 in matrix inversion", info, info);
144
 
  else
145
 
    Msg::Error("Wrong %d-th argument in matrix inversion", -info);
146
 
  return false;
147
 
}
148
 
 
149
 
 
150
 
template<> 
151
 
bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, gmshVector<double> &DR,
152
 
                             gmshVector<double> &DI, gmshMatrix<double> &VR,
153
 
                             bool sortRealPart)
154
 
{
155
 
  int N = size1(), info;
156
 
  int LWORK = 10 * N;
157
 
  double *work = new double[LWORK];
158
 
 
159
 
  F77NAME(dgeev)("V", "V", &N, _data, &N, DR._data, DI._data,
160
 
                 VL._data, &N, VR._data, &N, work, &LWORK, &info);
161
 
  
162
 
  delete [] work;
163
 
 
164
 
  if(info > 0)
165
 
    Msg::Error("QR Algorithm failed to compute all the eigenvalues", info, info);
166
 
  else if(info < 0)
167
 
    Msg::Error("Wrong %d-th argument in eig", -info);
168
 
 
169
 
  if (sortRealPart) {
170
 
    double tmp[8];
171
 
    // do permutations
172
 
    for (int i = 0; i < size1() - 1; i++) {
173
 
      int minR = i;
174
 
      for (int j = i + 1; j < size1(); j++) 
175
 
        if (fabs(DR(j)) < fabs(DR(minR))) minR = j;
176
 
      if ( minR != i ){
177
 
        tmp[0] = DR(i); tmp[1] = DI(i);
178
 
        tmp[2] = VL(0,i); tmp[3] = VL(1,i); tmp[4] = VL(2,i);
179
 
        tmp[5] = VR(0,i); tmp[6] = VR(1,i); tmp[7] = VR(2,i);
180
 
        
181
 
        DR(i) = DR(minR); DI(i) = DI(minR);
182
 
        VL(0,i) = VL(0,minR); VL(1,i) = VL(1,minR); VL(2,i) = VL(2,minR);
183
 
        VR(0,i) = VR(0,minR); VR(1,i) = VR(1,minR); VR(2,i) = VR(2,minR);
184
 
        
185
 
        DR(minR) = tmp[0]; DI(minR) = tmp[1];
186
 
        VL(0,minR) = tmp[2]; VL(1,minR) = tmp[3]; VL(2,minR) = tmp[4];
187
 
        VR(0,minR) = tmp[5]; VR(1,minR) = tmp[6]; VR(2,minR) = tmp[7];
188
 
      }
189
 
    }
190
 
  }
191
 
  
192
 
  return true;
193
 
}
194
 
 
195
 
template<> 
196
 
bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<double> &result)
197
 
{
198
 
  int N = size1(), nrhs = 1, lda = N, ldb = N, info;
199
 
  int *ipiv = new int[N];
200
 
  for(int i = 0; i < N; i++) result(i) = rhs(i);
201
 
  F77NAME(dgesv)(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
202
 
  delete [] ipiv;
203
 
  if(info == 0) return true;
204
 
  if(info > 0)
205
 
    Msg::Error("U(%d,%d)=0 in LU decomposition", info, info);
206
 
  else
207
 
    Msg::Error("Wrong %d-th argument in LU decomposition", -info);
208
 
  return false;
209
 
}
210
 
 
211
 
template<>
212
 
bool gmshMatrix<double>::invert(gmshMatrix<double> &result)
213
 
{
214
 
  int M = size1(), N = size2(), lda = size1(), info;
215
 
  int *ipiv = new int[std::min(M, N)];
216
 
  result = *this;
217
 
  F77NAME(dgetrf)(&M, &N, result._data, &lda, ipiv, &info);
218
 
  if(info == 0){
219
 
    int lwork = M * 4;
220
 
    double *work = new double[lwork];
221
 
    F77NAME(dgetri)(&M, result._data, &lda, ipiv, work, &lwork, &info);
222
 
    delete [] work;
223
 
  }
224
 
  delete [] ipiv;
225
 
  if(info == 0) return true;
226
 
  else if(info > 0)
227
 
    Msg::Error("U(%d,%d)=0 in matrix inversion", info, info);
228
 
  else
229
 
    Msg::Error("Wrong %d-th argument in matrix inversion", -info);
230
 
  return false;
231
 
}
232
 
 
233
 
template<> 
234
 
double gmshMatrix<double>::determinant() const
235
 
{
236
 
  gmshMatrix<double> tmp(*this);
237
 
  int M = size1(), N = size2(), lda = size1(), info;
238
 
  int *ipiv = new int[std::min(M, N)];
239
 
  F77NAME(dgetrf)(&M, &N, tmp._data, &lda, ipiv, &info);
240
 
  double det = 1.;
241
 
  if(info == 0){
242
 
    for(int i = 0; i < size1(); i++){
243
 
      det *= tmp(i, i);
244
 
      if(ipiv[i] != i + 1) det = -det;
245
 
    }
246
 
  }
247
 
  else if(info > 0)
248
 
    det = 0.;
249
 
  else
250
 
    Msg::Error("Wrong %d-th argument in matrix factorization", -info);
251
 
  delete [] ipiv;  
252
 
  return det;
253
 
}
254
 
 
255
 
template<> 
256
 
bool gmshMatrix<double>::svd(gmshMatrix<double> &V, gmshVector<double> &S)
257
 
{
258
 
  gmshMatrix<double> VT(V.size2(), V.size1());
259
 
  int M = size1(), N = size2(), LDA = size1(), LDVT = VT.size1(), info;
260
 
  int LWORK = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
261
 
  gmshVector<double> WORK(LWORK);
262
 
  F77NAME(dgesvd)("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
263
 
                  VT._data, &LDVT, WORK._data, &LWORK, &info);
264
 
  V = VT.transpose();
265
 
  if(info == 0) return true;
266
 
  if(info > 0)
267
 
    Msg::Error("SVD did not converge");
268
 
  else
269
 
    Msg::Error("Wrong %d-th argument in SVD decomposition", -info);
270
 
  return false;
271
 
}
272
 
 
273
 
#endif