9 #include "GmshConfig.h"
14 #define F77NAME(x) (x##_)
19 #if defined(HAVE_BLAS) && !defined(HAVE_EIGEN)
22 void F77NAME(daxpy)(
int *n,
double *alpha,
double *x,
int *incx,
double *y,
24 void F77NAME(zaxpy)(
int *n, std::complex<double> *alpha,
25 std::complex<double> *x,
int *incx, std::complex<double> *y,
27 void F77NAME(dcopy)(
int *n,
double *a,
int *inca,
double *b,
int *incb);
28 void F77NAME(zcopy)(
int *n, std::complex<double> *a,
int *inca,
29 std::complex<double> *b,
int *incb);
30 void F77NAME(dgemm)(
const char *transa,
const char *transb,
int *m,
int *n,
31 int *k,
double *alpha,
double *a,
int *lda,
double *b,
32 int *ldb,
double *beta,
double *
c,
int *ldc);
33 void F77NAME(zgemm)(
const char *transa,
const char *transb,
int *m,
int *n,
34 int *k, std::complex<double> *alpha,
35 std::complex<double> *a,
int *lda, std::complex<double> *b,
36 int *ldb, std::complex<double> *beta,
37 std::complex<double> *
c,
int *ldc);
38 void F77NAME(dgemv)(
const char *trans,
int *m,
int *n,
double *alpha,
double *a,
39 int *lda,
double *x,
int *incx,
double *beta,
double *y,
41 void F77NAME(zgemv)(
const char *trans,
int *m,
int *n,
42 std::complex<double> *alpha, std::complex<double> *a,
43 int *lda, std::complex<double> *x,
int *incx,
44 std::complex<double> *beta, std::complex<double> *y,
46 void F77NAME(dscal)(
int *n,
double *alpha,
double *x,
int *incx);
47 void F77NAME(zscal)(
int *n, std::complex<double> *alpha,
48 std::complex<double> *x,
int *incx);
62 F77NAME(zcopy)(&_r, m._data, &stride, _data, &stride);
68 int M =
_r, INCX = 1, INCY = 1;
74 const fullVector<std::complex<double> > &x, std::complex<double> alpha)
76 int M = _r, INCX = 1, INCY = 1;
77 F77NAME(zaxpy)(&M, &alpha, x._data, &INCX, _data, &INCY);
82 if(_r != m.
_r || _c != m.
_c) {
87 for(
int i = 0; i < N; ++i) _data[i] = m.
_data[i];
105 if(_r != m.
_r || _c != m.
_c) {
126 std::complex<double> ss(s);
127 F77NAME(zscal)(&N, &ss, _data, &stride);
132 for(
int i = 0; i < _r * _c; ++i) _data[i] *= s;
139 int M =
c.size1(), N =
c.size2(), K =
_c;
140 int LDA =
_r, LDB = b.
size1(), LDC =
c.size1();
141 double alpha = 1., beta = 0.;
143 (
"N",
"N", &M, &N, &K, &alpha,
_data, &LDA, b.
_data, &LDB, &beta,
c._data,
152 int M =
c.size1(), N =
c.size2(), K = _c;
153 int LDA = _r, LDB = b.
size1(), LDC =
c.size1();
154 std::complex<double> alpha = 1., beta = 0.;
156 (
"N",
"N", &M, &N, &K, &alpha, _data, &LDA, b.
_data, &LDB, &beta,
c._data,
164 int M =
_r, N =
_c, LDA =
_r, INCX = 1, INCY = 1;
165 double alpha = 1., beta = 0.;
167 (
"N", &M, &N, &alpha,
_data, &LDA, x.
_data, &INCX, &beta, y.
_data, &INCY);
175 int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
176 std::complex<double> alpha = 1., beta = 0.;
178 (
"N", &M, &N, &alpha, _data, &LDA, x._data, &INCX, &beta, y._data, &INCY);
185 int M =
_r, N =
_c, LDA =
_r, INCX = 1, INCY = 1;
186 double alpha = 1., beta = 1.;
188 (
"N", &M, &N, &alpha,
_data, &LDA, x.
_data, &INCX, &beta, y.
_data, &INCY);
196 int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
197 std::complex<double> alpha = 1., beta = 1.;
199 (
"N", &M, &N, &alpha, _data, &LDA, x._data, &INCX, &beta, y._data, &INCY);
204 const int ncol,
const int fcol,
205 const int alpha_,
const int beta_,
208 int M = 1, N = ncol, K = b.
size1();
209 int LDA =
_r, LDB = b.
size1(), LDC = 1;
210 double alpha = alpha_, beta = beta_;
212 (
"N",
"N", &M, &N, &K, &alpha,
_data, &LDA, &(b.
_data[fcol * K]), &LDB, &beta,
213 &(
c._data[fcol]), &LDC);
218 double alpha,
double beta,
221 int M =
_r, N =
_c, LDA =
_r, INCX = 1, INCY = 1;
223 (
"T", &M, &N, &alpha,
_data, &LDA, x.
_data, &INCX, &beta, y.
_data, &INCY);
229 double beta,
bool transposeA,
bool transposeB)
234 (transposeA ?
"T" :
"N", transposeB ?
"T" :
"N", &M, &N, &K, &alpha, a.
_data,
241 const fullMatrix<std::complex<double> > &b, std::complex<double> alpha,
242 std::complex<double> beta,
bool transposeA,
bool transposeB)
244 int M = size1(), N = size2(), K = transposeA ? a.
size1() : a.
size2();
245 int LDA = a.
size1(), LDB = b.
size1(), LDC = size1();
247 (transposeA ?
"T" :
"N", transposeB ?
"T" :
"N", &M, &N, &K, &alpha, a.
_data,
248 &LDA, b.
_data, &LDB, &beta, _data, &LDC);
254 int M =
_r *
_c, INCX = 1, INCY = 1;
260 #if defined(HAVE_LAPACK) && !defined(HAVE_EIGEN)
263 void F77NAME(dgesv)(
int *N,
int *nrhs,
double *A,
int *lda,
int *ipiv,
264 double *b,
int *ldb,
int *info);
265 void F77NAME(dgetrf)(
int *M,
int *N,
double *A,
int *lda,
int *ipiv,
int *info);
266 void F77NAME(dgetrs)(
char *trans,
int *N,
int *nrhs,
double *A,
int *lda,
267 int *ipiv,
double *b,
int *ldb,
int *info);
268 void F77NAME(dgetri)(
int *M,
double *A,
int *lda,
int *ipiv,
double *work,
269 int *lwork,
int *info);
270 void F77NAME(dgesvd)(
const char *jobu,
const char *jobvt,
int *M,
int *N,
271 double *A,
int *lda,
double *
S,
double *U,
int *ldu,
272 double *VT,
int *ldvt,
double *work,
int *lwork,
274 void F77NAME(dgeev)(
const char *jobvl,
const char *jobvr,
int *n,
double *a,
275 int *lda,
double *wr,
double *wi,
double *vl,
int *ldvl,
276 double *vr,
int *ldvr,
double *work,
int *lwork,
int *info);
283 int N =
size1(), nrhs = 1, lda = N, ldb = N, info;
284 int *ipiv =
new int[N];
285 for(
int i = 0; i < N; i++) result(i) = rhs(i);
288 if(info == 0)
return true;
295 ipiv.
resize(std::min(M, N));
297 if(info == 0)
return true;
306 int N =
size1(), nrhs = 1, lda = N, ldb = N, info;
308 for(
int i = 0; i < N; i++) result(i) = rhs(i);
310 (&trans, &N, &nrhs,
_data, &lda, &ipiv(0), result.
_data, &ldb, &info);
311 if(info == 0)
return true;
318 int *ipiv =
new int[std::min(M, N)];
319 if(result.
size2() != M || result.
size1() != N) {
321 result.
resize(M, N,
false);
323 Msg::Error(
"FullMatrix: Bad dimension, I cannot write in proxy");
328 F77NAME(dgetrf)(&M, &N, result.
_data, &lda, ipiv, &info);
331 double *work =
new double[lwork];
332 F77NAME(dgetri)(&M, result.
_data, &lda, ipiv, work, &lwork, &info);
339 Msg::Warning(
"U(%d,%d)=0 in matrix inversion", info, info);
341 Msg::Error(
"Wrong %d-th argument in matrix inversion", -info);
347 int N =
size1(), nrhs = N, lda = N, ldb = N, info;
348 int *ipiv =
new int[N];
349 double *invA =
new double[N * N];
351 for(
int i = 0; i < N * N; i++) invA[i] = 0.;
352 for(
int i = 0; i < N; i++) invA[i * N + i] = 1.;
354 F77NAME(dgesv)(&N, &nrhs,
_data, &lda, ipiv, invA, &ldb, &info);
355 memcpy(
_data, invA, N * N *
sizeof(
double));
360 if(info == 0)
return true;
362 Msg::Warning(
"U(%d,%d)=0 in matrix in place inversion", info, info);
364 Msg::Error(
"Wrong %d-th argument in matrix inversion", -info);
372 int *ipiv =
new int[std::min(M, N)];
373 F77NAME(dgetrf)(&M, &N, tmp._data, &lda, ipiv, &info);
376 for(
int i = 0; i <
size1(); i++) {
378 if(ipiv[i] != i + 1) det = -det;
384 Msg::Error(
"Wrong %d-th argument in matrix factorization", -info);
394 int N =
size1(), info;
396 double *work =
new double[lwork];
398 (
"V",
"V", &N,
_data, &N, DR.
_data, DI.
_data, VL.
_data, &N, VR.
_data, &N,
399 work, &lwork, &info);
403 Msg::Error(
"QR Algorithm failed to compute all the eigenvalues", info,
406 Msg::Error(
"Wrong %d-th argument in eig", -info);
407 else if(sortRealPart)
418 int lwork = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
421 (
"O",
"A", &M, &N,
_data, &LDA,
S._data,
_data, &LDA, VT._data, &LDVT,
422 WORK._data, &lwork, &info);
424 if(info == 0)
return true;
428 Msg::Error(
"Wrong %d-th argument in SVD decomposition", -info);
446 for(
int i = 0; i < _r; ++i)
447 n += _data[i].real() * _data[i].real() + _data[i].imag() * _data[i].imag();
448 return std::complex<double>(sqrt(n), 0.);
453 const std::string format)
const
455 std::string rformat = (format ==
"") ?
"%12.5E " : format;
456 printf(
"double %s[%d]=\n", name.c_str(),
size());
458 for(
int I = 0; I <
size(); I++) { printf(rformat.c_str(), (*
this)(I)); }
464 const std::string format)
const
466 std::string rformat = (format ==
"") ?
"%12d " : format;
467 printf(
"double %s[%d]=\n", name.c_str(), size());
469 for(
int I = 0; I < size(); I++) { printf(rformat.c_str(), (*
this)(I)); }
475 const std::string &format)
const
477 std::string rformat = (format ==
"") ?
"%12.5E " : format;
480 printf(
"double %s [ %d ][ %d ]= { \n", name.c_str(), ni, nj);
481 for(
int I = 0; I < ni; I++) {
483 for(
int J = 0; J < nj; J++) {
484 printf(rformat.c_str(), (*
this)(I, J));
485 if(J != nj - 1) printf(
",");
497 const std::string &format)
const
499 std::string rformat = (format ==
"") ?
"%12d " : format;
502 printf(
"int %s [ %d ][ %d ]= { \n", name.c_str(), ni, nj);
503 for(
int I = 0; I < ni; I++) {
505 for(
int J = 0; J < nj; J++) {
506 printf(rformat.c_str(), (*
this)(I, J));
507 if(J != nj - 1) printf(
",");