linalg.cpp
1.61 KB
1
2
3
4
5
6
7
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
#include "linalg.h"
// This file contains a set of wrapper functions that are linked to the corresponding functions in CLAPACK.
extern "C" {
#include "f2c.h"
#include "clapack.h"
#include "cblas.h"
}
void LINALG_zgetrf(
int M,
int N,
std::complex<double>* A,
int LDA,
int* IPIV)
{
integer INFO;
zgetrf_((integer*)&M, (integer*)&N, (doublecomplex*)A, (integer*)&LDA, (integer*)IPIV, &INFO);
}
void LINALG_zgetri(
size_t N,
std::complex<double>* A,
int LDA,
int* IPIV)
{
integer LWORK = -1;
std::complex<double> WORK[1];
integer INFO;
zgetri_((integer*)&N, (doublecomplex*)A, (integer*)&LDA, (integer*)IPIV, (doublecomplex*)WORK, &LWORK, &INFO);
}
void LINALG_inverse(std::complex<double>* A, int N)
{
int* IPIV = new int[N + (size_t)1];
integer LWORK = N * N;
std::complex<double>* WORK = new std::complex<double>[LWORK];
integer INFO;
zgetrf_((integer*)&N, (integer*)&N, (doublecomplex*)A, (integer*)&N, (integer*)IPIV, &INFO);
zgetri_((integer*)&N, (doublecomplex*)A, (integer*)&N, (integer*)IPIV, (doublecomplex*)WORK, &LWORK, &INFO);
delete[] IPIV;
delete[] WORK;
}
void LINALG_zgemm(
const int M, //A(M*K) B(K*N)
const int N,
const int K,
std::complex<double>* A,
const int LDA, //=K
std::complex<double>* B,
const int LDB, //=N
std::complex<double>* C,
const int LDC) //=columns of C.
{
std::complex<double> alpha = 1;
std::complex<double> beta = 0;
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, (OPENBLAS_CONST blasint)M, (OPENBLAS_CONST blasint)N, (OPENBLAS_CONST blasint)K,
&alpha, A, (OPENBLAS_CONST blasint)LDA, B, (OPENBLAS_CONST blasint)LDB, &beta, C, (OPENBLAS_CONST blasint)LDC);
}