Newer
Older
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include <complex>
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#include "GmshConfig.h"
#include "GmshMatrix.h"
#include "GmshMessage.h"
#if defined(HAVE_BLAS)
extern "C" {
void dgemm_(const char *transa, const char *transb, int *m, int *n, int *k,
double *alpha, double *a, int *lda,
double *b, int *ldb, double *beta,
double *c, int *ldc);
void zgemm_(const char *transa, const char *transb, int *m, int *n, int *k,
std::complex<double> *alpha, std::complex<double> *a, int *lda,
std::complex<double> *b, int *ldb, std::complex<double> *beta,
std::complex<double> *c, int *ldc);
void dgemv_(const char *trans, int *m, int *n,
double *alpha, double *a, int *lda,
double *x, int *incx, double *beta,
double *y, int *incy);
void zgemv_(const char *trans, int *m, int *n,
std::complex<double> *alpha, std::complex<double> *a, int *lda,
std::complex<double> *x, int *incx, std::complex<double> *beta,
std::complex<double> *y, int *incy);
}
template<>
void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c)
{
int M = c.size1(), N = c.size2(), K = _c;
int LDA = _r, LDB = b.size1(), LDC = c.size1();
double alpha = 1., beta = 0.;
dgemm_("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB,
&beta, c._data, &LDC);
}
template<>
void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<double> > &b,
gmshMatrix<std::complex<double> > &c)
{
int M = c.size1(), N = c.size2(), K = _c;
int LDA = _r, LDB = b.size1(), LDC = c.size1();
std::complex<double> alpha = 1., beta = 0.;
zgemm_("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB,
&beta, c._data, &LDC);
}
template<>
void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b,
double alpha, double beta)
{
int M = size1(), N = size2(), K = a.size2();
int LDA = a.size1(), LDB = b.size1(), LDC = size1();
dgemm_("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB,
&beta, _data, &LDC);
}
template<>
void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &a,
gmshMatrix<std::complex<double> > &b,
std::complex<double> alpha,
std::complex<double> beta)
{
int M = size1(), N = size2(), K = a.size2();
int LDA = a.size1(), LDB = b.size1(), LDC = size1();
zgemm_("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB,
&beta, _data, &LDC);
}
template<>
void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y)
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
double alpha = 1., beta = 0.;
dgemv_("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
&beta, y._data, &INCY);
}
template<>
void gmshMatrix<std::complex<double> >::mult(const gmshVector<std::complex<double> > &x,
gmshVector<std::complex<double> > &y)
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
std::complex<double> alpha = 1., beta = 0.;
zgemv_("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
&beta, y._data, &INCY);
}
#endif
#if defined(HAVE_LAPACK)
extern "C" {
void dgesv_(int *N, int *nrhs, double *A, int *lda, int *ipiv,
double *b, int *ldb, int *info);
void dgetrf_(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
void dgetri_(int *M, double *A, int *lda, int *ipiv, double *work, int *lwork,
int *info);
void dgesvd_(const char* jobu, const char *jobvt, int *M, int *N,
double *A, int *lda, double *S, double* U, int *ldu,
double *VT, int *ldvt, double *work, int *lwork, int *info);
void dgeev_(const char *jobvl, const char *jobvr,
int *n, double *a, int *lda,
double *wr, double *wi,
double *vl, int *ldvl,
double *vr, int *ldvr,
double *work, int *lwork,
int *info);
template<>
bool gmshMatrix<double>::invertInPlace()
{
int N = size1(), nrhs = N, lda = N, ldb = N, info;
int *ipiv = new int[N];
double * invA = new double[N*N];
for (size_t i = 0; i < N * N; i++) invA[i] = 0.;
for (size_t i = 0; i < N; i++) invA[i * N + i] = 1.;
dgesv_(&N, &nrhs, _data, &lda, ipiv, invA, &ldb, &info);
delete [] invA;
delete [] ipiv;
if(info == 0) return true;
if(info > 0)
Msg::Error("U(%d,%d)=0 in matrix inversion", info, info);
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
else
Msg::Error("Wrong %d-th argument in matrix inversion", -info);
return false;
}
template<>
bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, // left eigenvectors
gmshVector<double> &DR, // Real part of eigenvalues
gmshVector<double> &DI, // Im part of eigenvalues
gmshMatrix<double> &VR )
{
int N = size1(), info;
int LWORK = 10*N;
double * work = new double[LWORK];
dgeev_("V","V",
&N,_data,
&N,DR._data,DI._data,
VL._data,&N,
VR._data,&N,
work,&LWORK,&info);
delete [] work;
if(info == 0) return true;
if(info > 0)
Msg::Error("QR Algorithm failed to compute all the eigenvalues", info, info);
else
Msg::Error("Wrong %d-th argument in eig", -info);
return false;
}
template<>
bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<double> &result)
{
int N = size1(), nrhs = 1, lda = N, ldb = N, info;
int *ipiv = new int[N];
for(int i = 0; i < N; i++) result(i) = rhs(i);
dgesv_(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
delete [] ipiv;
if(info == 0) return true;
if(info > 0)
Msg::Error("U(%d,%d)=0 in LU decomposition", info, info);
else
Msg::Error("Wrong %d-th argument in LU decomposition", -info);
template<>
bool gmshMatrix<double>::invert(gmshMatrix<double> &result)
{
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
result = *this;
dgetrf_(&M, &N, result._data, &lda, ipiv, &info);
if(info == 0){
int lwork = M * 4;
double *work = new double[lwork];
dgetri_(&M, result._data, &lda, ipiv, work, &lwork, &info);
delete [] work;
}
delete [] ipiv;
if(info == 0) return true;
else if(info > 0)
else
Msg::Error("Wrong %d-th argument in matrix inversion", -info);
return false;
template<>
double gmshMatrix<double>::determinant() const
{
gmshMatrix<double> tmp(*this);
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
dgetrf_(&M, &N, tmp._data, &lda, ipiv, &info);
double det = 1.;
if(info == 0){
for(int i = 0; i < size1(); i++){
det *= tmp(i, i);
if(ipiv[i] != i + 1) det = -det;
}
Msg::Error("Wrong %d-th argument in matrix factorization", -info);
return det;
}
template<>
bool gmshMatrix<double>::svd(gmshMatrix<double> &V, gmshVector<double> &S)
{
gmshMatrix<double> VT(V.size2(), V.size1());
int M = size1(), N = size2(), LDA = size1(), LDVT = VT.size1(), info;
int LWORK = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
gmshVector<double> WORK(LWORK);
dgesvd_("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
VT._data, &LDVT, WORK._data, &LWORK, &info);
V = VT.transpose();
if(info == 0) return true;
if(info > 0)
Msg::Error("SVD did not converge");
else
Msg::Error("Wrong %d-th argument in SVD decomposition", -info);