// LAPACK: functions for computing eigenvalues and eigenvectors.
package scientific.linear
import std.unittest.*
import std.unittest.testmacro.*
import scientific.numbers.*
/* Compute eigenvalue and eigenvectors of a general matrix. */
foreign func LAPACKE_dgeev(
matrix_layout: IntNative, // row or column major
jobvl: UInt8, // whether to compute left eigenvectors
jobvr: UInt8, // whether to compute right eigenvectors
n: IntNative, // size of A
a: CPointer<Unit>, // matrix A (over-written, Float64)
lda: IntNative, // leading dimension of A
wr: CPointer<Unit>, // real part of eigenvalues, Float64
wi: CPointer<Unit>, // imaginary part of eigenvalues, Float64
vl: CPointer<Unit>, // left eigenvectors, Float64
ldvl: IntNative, // leading dimension of vl
vr: CPointer<Unit>, // right eigenvectors, Float64
ldvr: IntNative // leading dimension of vr
): Int64
/* Float32 version */
foreign func LAPACKE_sgeev(
matrix_layout: IntNative, // row or column major
jobvl: UInt8, // whether to compute left eigenvectors
jobvr: UInt8, // whether to compute right eigenvectors
n: IntNative, // size of A
a: CPointer<Unit>, // matrix A (over-written, Float32)
lda: IntNative, // leading dimension of A
wr: CPointer<Unit>, // real part of eigenvalues, Float32
wi: CPointer<Unit>, // imaginary part of eigenvalues, Float32
vl: CPointer<Unit>, // left eigenvectors, Float32
ldvl: IntNative, // leading dimension of vl
vr: CPointer<Unit>, // right eigenvectors, Float32
ldvr: IntNative // leading dimension of vr
): Int64
/* Complex64 version */
foreign func LAPACKE_zgeev(
matrix_layout: IntNative, // row or column major
jobvl: UInt8, // whether to compute left eigenvectors
jobvr: UInt8, // whether to compute right eigenvectors
n: IntNative, // size of A
a: CPointer<Unit>, // matrix A (over-written, Complex64)
lda: IntNative, // leading dimension of A
w: CPointer<Unit>, // eigenvalues, Complex64
vl: CPointer<Unit>, // left eigenvectors, Complex64
ldvl: IntNative, // leading dimension of vl
vr: CPointer<Unit>, // right eigenvectors, Complex64
ldvr: IntNative // leading dimension of vr
): Int64
/* Complex32 version */
foreign func LAPACKE_cgeev(
matrix_layout: IntNative, // row or column major
jobvl: UInt8, // whether to compute left eigenvectors
jobvr: UInt8, // whether to compute right eigenvectors
n: IntNative, // size of A
a: CPointer<Unit>, // matrix A (over-written, Complex32)
lda: IntNative, // leading dimension of A
w: CPointer<Unit>, // eigenvalues, Complex32
vl: CPointer<Unit>, // left eigenvectors, Complex32
ldvl: IntNative, // leading dimension of vl
vr: CPointer<Unit>, // right eigenvectors, Complex32
ldvr: IntNative // leading dimension of vr
): Int64
/* DSYEVR computes the eigenvalues and, optionally, the left and/or
right eigenvectors for SY matrices */
foreign func LAPACKE_dsyevr(
matrix_layout: IntNative,
jobz: UInt8,
range: UInt8,
uplo: UInt8,
n: IntNative,
a: CPointer<Unit>, // Float64
lda: IntNative,
vl: Float64,
vu: Float64,
il: IntNative,
iu: IntNative,
abstol: Float64,
m: CPointer<Unit>, // Int64
w: CPointer<Unit>, // Float64
z: CPointer<Unit>, // Float64
ldz: IntNative,
isuppz: CPointer<Unit> // Int64
): Int64
foreign func LAPACKE_dsbevd(
matrix_layout: IntNative,
jobz: UInt8,
uplo: UInt8,
n: IntNative,
kd: IntNative,
ab: CPointer<Unit>, // Float64
ldab: IntNative,
w: CPointer<Unit>, // Float64
z: CPointer<Unit>, // Float64
ldz: IntNative
): Int64
/* DSBEVX computes selected eigenvalues and, optionally, eigenvectors
of a real symmetric band matrix A. */
foreign func LAPACKE_dsbevx(
matrix_layout: IntNative,
jobz: UInt8,
range: UInt8,
uplo: UInt8,
n: IntNative,
kd: IntNative,
ab: CPointer<Unit>, // Float64
ldab: IntNative,
q: CPointer<Unit>, // Float64
ldq: IntNative,
vl: Float64,
vu: Float64,
il: IntNative,
iu: IntNative,
abstol: Float64,
m: CPointer<Unit>, // Int64
w: CPointer<Unit>, // Float64
z: CPointer<Unit>, // Float64
ldz: IntNative,
ifail: CPointer<Unit> // Int64
): Int64
/* DSTEMR computes selected eigenvalues and, optionally, eigenvectors
of a real symmetric tridiagonal matrix T. */
foreign func LAPACKE_dstemr(
matrix_layout: IntNative,
jobz: UInt8,
range: UInt8,
n: IntNative,
d: CPointer<Unit>, // Float64
e: CPointer<Unit>, // Float64
vl: Float64,
vu: Float64,
il: IntNative,
iu: IntNative,
m: CPointer<Unit>, // Int64
w: CPointer<Unit>, // Float64
z: CPointer<Unit>, // Float64
ldz: IntNative,
nzc: IntNative,
isuppz: CPointer<Unit>, // Int64
tryrac: CPointer<Unit> // Int64
): Int64
func testDstemr() {
var d = vector<Float64>(
[3.0, 3.0, 3.0, 3.0]
)
var e = vector<Float64>(
[-1.0, -1.0, -1.0, -1.0]
)
let N: Int64 = d.size()
var vl: Float64 = 0.0
var vu: Float64 = 0.0
var il: Int64 = 0
var iu: Int64 = 0
let ldz: Int64 = N
let nzc: Int64 = N
var w = zeros<Float64>(N)
var m = zeros<Int64>(N)
var z = zeros<Float64>(ldz, N)
var isuppz = zeros<Int64>(2 * N)
var tryrac = zeros<Int64>(N)
let info = unsafe { LAPACKE_dstemr(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_ALL_EIGENVALUES,
IntNative(N), d.ptr, e.ptr, vl, vu, IntNative(il), IntNative(iu), m.ptr,
w.ptr, z.ptr, IntNative(ldz), IntNative(nzc), isuppz.ptr, tryrac.ptr) }
var expected_w = vector<Float64>(
[1.381966, 2.381966, 3.618034, 4.618034]
)
var expected_z = matrix<Float64>(
[[0.371748, 0.601501, 0.601501, 0.371748],
[0.601501, 0.371748, -0.371748, -0.601501],
[0.601501, -0.371748, -0.371748, 0.601501],
[0.371748, -0.601501, 0.601501, -0.371748]]
)
assertApproxEqual(w, expected_w, atol:1e-6)
assertApproxEqual(z, expected_z, atol:1e-6)
}
public func eighTridiagonal(d: Vector<Float64>, e: Vector<Float64>, eigvals_only!:Bool = false):
(Vector<Float64>, Matrix<Float64>) {
let N: Int64 = d.size()
var vl: Float64 = 0.0
var vu: Float64 = 0.0
var il: Int64 = 0
var iu: Int64 = 0
let ldz: Int64 = N
let nzc: Int64 = N
var w = zeros<Float64>(N)
var m = zeros<Int64>(N)
var z = zeros<Float64>(ldz, N)
var isuppz = zeros<Int64>(2 * N)
var tryrac = zeros<Int64>(N)
if (eigvals_only) {
let info = unsafe { LAPACKE_dstemr(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_ALL_EIGENVALUES,
IntNative(N), d.ptr, e.ptr, vl, vu, IntNative(il), IntNative(iu), m.ptr,
w.ptr, z.ptr, IntNative(ldz), IntNative(nzc), isuppz.ptr, tryrac.ptr) }
} else {
let info = unsafe { LAPACKE_dstemr(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_ALL_EIGENVALUES,
IntNative(N), d.ptr, e.ptr, vl, vu, IntNative(il), IntNative(iu), m.ptr,
w.ptr, z.ptr, IntNative(ldz), IntNative(nzc), isuppz.ptr, tryrac.ptr) }
}
return (w, z)
}
public func eigvalshTridiagonal(d: Vector<Float64>, e: Vector<Float64>): Vector<Float64> {
let (w, z) = eighTridiagonal(d, e, eigvals_only:true)
return w
}
func testDsbevd() {
var ab = matrix<Float64>(
[[1.0, 2.0, 3.0, 4.0],
[5.0, 5.0, 5.0, 0.0],
[2.0, 2.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0]]
)
let N: Int64 = ab.getCols()
let kd: Int64 = ab.getRows() - 1
let ldab: Int64 = kd + 1
let ldz: Int64 = N
var w = zeros<Float64>(N)
var z = zeros<Float64>(ldz, N)
var ifail = zeros<Int64>(N)
let info = unsafe { LAPACKE_dsbevd(LAPACK_COL_MAJOR, LAPACK_VECTOR,
LAPACK_LOWER_TRIANGULAR, IntNative(N), IntNative(kd), ab.ptr, IntNative(ldab),
w.ptr,z.ptr, IntNative(ldz)) }
}
public func eigBanded(ab: Matrix<Float64>, eigvals_only!:Bool = false): (Vector<Float64>, Matrix<Float64>) {
let N: Int64 = ab.getCols()
let kd: Int64 = ab.getRows() - 1
let ldab: Int64 = kd + 1
let ldz: Int64 = N
var w = zeros<Float64>(N)
var z = zeros<Float64>(ldz, N)
var ifail = zeros<Int64>(N)
if (eigvals_only) {
let info = unsafe { LAPACKE_dsbevd(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR,
LAPACK_LOWER_TRIANGULAR, IntNative(N), IntNative(kd), ab.ptr, IntNative(ldab),
w.ptr,z.ptr, IntNative(ldz)) }
} else {
let info = unsafe { LAPACKE_dsbevd(LAPACK_COL_MAJOR, LAPACK_VECTOR,
LAPACK_LOWER_TRIANGULAR, IntNative(N), IntNative(kd), ab.ptr, IntNative(ldab),
w.ptr,z.ptr, IntNative(ldz)) }
}
return (w, z)
}
public func eigvalsBanded(ab: Matrix<Float64>): Vector<Float64> {
let (w, z) = eigBanded(ab, eigvals_only:true)
return w
}
func testDsbevx() {
var ab = matrix<Float64>(
[[1.0, 2.0, 3.0, 4.0],
[5.0, 5.0, 5.0, 0.0],
[2.0, 2.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0]]
)
let N: Int64 = ab.getCols()
let kd: Int64 = ab.getRows() - 1
let ldab: Int64 = kd + 1
let ldq: Int64 = N
var q = zeros<Float64>(ldq, N)
var vl: Float64 = 0.0
var vu: Float64 = 0.0
var il: Int64 = 0
var iu: Int64 = 0
let abstol: Float64 = 1e-6
let ldz: Int64 = N
var w = zeros<Float64>(N)
var m = zeros<Int64>(N)
var z = zeros<Float64>(ldz, N)
var ifail = zeros<Int64>(N)
let info = unsafe { LAPACKE_dsbevx(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_ALL_EIGENVALUES,
LAPACK_LOWER_TRIANGULAR, IntNative(N), IntNative(kd), ab.ptr, IntNative(ldab),
q.ptr, IntNative(ldq), vl, vu, IntNative(il), IntNative(iu), abstol, m.ptr,
w.ptr,z.ptr, IntNative(ldz), ifail.ptr) }
}
public func eigh(a: Matrix<Float64>, eigvals_only!:Bool = false): (Vector<Float64>, Matrix<Float64>) {
let N: Int64 = a.getCols()
var vl: Float64 = 0.0
var vu: Float64 = 0.0
var il: Int64 = 0
var iu: Int64 = 0
let abstol: Float64 = 1e-6
let ldz: Int64 = N
var w = zeros<Float64>(N)
var m = zeros<Int64>(N)
var z = zeros<Float64>(ldz, N)
var isuppz = zeros<Int64>(2 * N)
if (eigvals_only) {
let info = unsafe { LAPACKE_dsyevr(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_ALL_EIGENVALUES,
LAPACK_LOWER_TRIANGULAR, IntNative(N), a.ptr, IntNative(N),
vl, vu, IntNative(il), IntNative(iu), abstol, m.ptr, w.ptr,z.ptr,
IntNative(ldz), isuppz.ptr) }
} else {
let info = unsafe { LAPACKE_dsyevr(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_ALL_EIGENVALUES,
LAPACK_LOWER_TRIANGULAR, IntNative(N), a.ptr, IntNative(N),
vl, vu, IntNative(il), IntNative(iu), abstol, m.ptr, w.ptr,z.ptr,
IntNative(ldz), isuppz.ptr) }
}
return (w, z)
}
public func eigvalsh(a: Matrix<Float64>): Vector<Float64> {
let (w, z) = eigh(a, eigvals_only:true)
return w
}
func testDsyevr() {
var a = matrix<Float64>(
[[34.0, -4.0, -10.0, -7.0, 2.0],
[-4.0, 7.0, 2.0, 12.0, 0.0],
[-10.0, 2.0, 44.0, 2.0, -19.0],
[-7.0, 12.0, 2.0, 79.0, -34.0],
[2.0, 0.0, -19.0, -34.0, 29.0]]
)
let N: Int64 = a.getCols()
var vl: Float64 = 0.0
var vu: Float64 = 0.0
var il: Int64 = 0
var iu: Int64 = 0
let abstol: Float64 = 1e-6
let ldz: Int64 = N
var w = zeros<Float64>(N)
var m = zeros<Int64>(N)
var z = zeros<Float64>(ldz, N)
var isuppz = zeros<Int64>(2 * N)
let info = unsafe { LAPACKE_dsyevr(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_ALL_EIGENVALUES,
LAPACK_LOWER_TRIANGULAR, IntNative(N), a.ptr, IntNative(N),
vl, vu, IntNative(il), IntNative(iu), abstol, m.ptr, w.ptr,z.ptr,
IntNative(ldz), isuppz.ptr) }
var expected_w = vector<Float64>(
[0.000001, 9.119382, 31.587600, 51.433367, 100.859651]
)
var expected_z = matrix<Float64>(
[[-0.042499, -0.243066, 0.892827, -0.349513, 0.140762],
[ 0.662479, -0.721978, -0.158056, -0.029858, -0.118300],
[-0.281040, -0.315582, 0.263317, 0.840945, -0.211894],
[-0.353773, -0.155937, -0.070130, -0.367218, -0.843067],
[-0.595967, -0.543836, -0.321907, -0.186846, 0.458836]]
)
assertApproxEqual(w, expected_w, atol:1e-6)
assertApproxEqual(z, expected_z, atol:1e-6)
}
func testDgeev() {
var a = matrix<Float64>(
[[0.0, -1.0],
[1.0, 0.0]]
)
let N: Int64 = a.getCols()
var wr = zeros<Float64>(N)
var wi = zeros<Float64>(N)
var vl = zeros<Float64>(N, N)
var vr = zeros<Float64>(N, N)
let info = unsafe { LAPACKE_dgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_VECTOR,
IntNative(N), a.ptr, IntNative(N), wr.ptr, wi.ptr, vl.ptr,
IntNative(N), vr.ptr, IntNative(N)) }
var expected_wr = vector<Float64>(
[0.0, 0.0]
)
var expected_wi = vector<Float64>(
[1.0, -1.0]
)
assertApproxEqual(wr, expected_wr, atol:1e-6)
assertApproxEqual(wi, expected_wi, atol:1e-6)
}
public func eig(a: Matrix<Float64>, right!:Bool = true): (Vector<Complex64>, Matrix<Complex64>) {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var wr = empty<Float64>(N)
var wi = empty<Float64>(N)
var v_read = empty<Float64>(N, N)
var v_empty = empty<Float64>(N, N)
let a_copy = a.copy()
var info: Int64
if (right) {
info = unsafe { LAPACKE_dgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), wr.ptr, wi.ptr, v_empty.ptr,
IntNative(N), v_read.ptr, IntNative(N)) }
} else {
info = unsafe { LAPACKE_dgeev(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), wr.ptr, wi.ptr, v_read.ptr,
IntNative(N), v_empty.ptr, IntNative(N)) }
}
// Now read the eigenvalues and eigenvectors
let eigvals = empty<Complex64>(N)
let eigvecs = empty<Complex64>(N, N)
var i: Int64 = 0
while (i < N) {
eigvals[i] = Complex64(wr[i], wi[i])
if (wi[i] > 0.0) {
// Pair of complex eigenvalues.
eigvals[i+1] = Complex64(wr[i+1], wi[i+1])
for (j in 0..N) {
eigvecs[j,i] = Complex64(v_read[j,i], v_read[j,i+1])
eigvecs[j,i+1] = Complex64(v_read[j,i], -v_read[j,i+1])
}
i += 2
} else {
// Real eigenvalue.
for (j in 0..N) {
eigvecs[j,i] = Complex64(v_read[j,i], 0.0)
}
i += 1
}
}
return (eigvals, eigvecs)
}
public func eigvals(a: Matrix<Float64>): Vector<Complex64> {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var wr = empty<Float64>(N)
var wi = empty<Float64>(N)
var v_e1 = empty<Float64>(N, N)
var v_e2 = empty<Float64>(N, N)
let a_copy = a.copy()
var info = unsafe { LAPACKE_dgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), wr.ptr, wi.ptr, v_e1.ptr,
IntNative(N), v_e2.ptr, IntNative(N)) }
// Now read the eigenvalues
let eigvals = empty<Complex64>(N)
for (i in 0..N) {
eigvals[i] = Complex64(wr[i], wi[i])
}
return eigvals
}
public func eig(a: Matrix<Float32>, right!:Bool = true): (Vector<Complex32>, Matrix<Complex32>) {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var wr = empty<Float32>(N)
var wi = empty<Float32>(N)
var v_read = empty<Float32>(N, N)
var v_empty = empty<Float32>(N, N)
let a_copy = a.copy()
var info: Int64
if (right) {
info = unsafe { LAPACKE_sgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), wr.ptr, wi.ptr, v_empty.ptr,
IntNative(N), v_read.ptr, IntNative(N)) }
} else {
info = unsafe { LAPACKE_sgeev(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), wr.ptr, wi.ptr, v_read.ptr,
IntNative(N), v_empty.ptr, IntNative(N)) }
}
// Now read the eigenvalues and eigenvectors
let eigvals = empty<Complex32>(N)
let eigvecs = empty<Complex32>(N, N)
var i: Int64 = 0
while (i < N) {
eigvals[i] = Complex32(wr[i], wi[i])
if (wi[i] > Float32(0.0)) {
// Pair of complex eigenvalues.
eigvals[i+1] = Complex32(wr[i+1], wi[i+1])
for (j in 0..N) {
eigvecs[j,i] = Complex32(v_read[j,i], v_read[j,i+1])
eigvecs[j,i+1] = Complex32(v_read[j,i], -v_read[j,i+1])
}
i += 2
} else {
// Real eigenvalue.
for (j in 0..N) {
eigvecs[j,i] = Complex32(v_read[j,i], 0.0)
}
i += 1
}
}
return (eigvals, eigvecs)
}
public func eigvals(a: Matrix<Float32>): Vector<Complex32> {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var wr = empty<Float32>(N)
var wi = empty<Float32>(N)
var v_e1 = empty<Float32>(N, N)
var v_e2 = empty<Float32>(N, N)
let a_copy = a.copy()
var info = unsafe { LAPACKE_sgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), wr.ptr, wi.ptr, v_e1.ptr,
IntNative(N), v_e2.ptr, IntNative(N)) }
// Now read the eigenvalues
let eigvals = empty<Complex32>(N)
for (i in 0..N) {
eigvals[i] = Complex32(wr[i], wi[i])
}
return eigvals
}
public func eig(a: Matrix<Complex64>, right!:Bool = true): (Vector<Complex64>, Matrix<Complex64>) {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var w = empty<Complex64>(N)
var v_read = empty<Complex64>(N, N)
var v_empty = empty<Complex64>(N, N)
let a_copy = a.copy()
var info: Int64
if (right) {
info = unsafe { LAPACKE_zgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), w.ptr, v_empty.ptr,
IntNative(N), v_read.ptr, IntNative(N)) }
} else {
info = unsafe { LAPACKE_zgeev(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), w.ptr, v_read.ptr,
IntNative(N), v_empty.ptr, IntNative(N)) }
}
return (w, v_read)
}
public func eigvals(a: Matrix<Complex64>): Vector<Complex64> {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var w = empty<Complex64>(N)
var v_e1 = empty<Complex64>(N, N)
var v_e2 = empty<Complex64>(N, N)
let a_copy = a.copy()
var info = unsafe { LAPACKE_zgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), w.ptr, v_e1.ptr,
IntNative(N), v_e2.ptr, IntNative(N)) }
return w
}
public func eig(a: Matrix<Complex32>, right!:Bool = true): (Vector<Complex32>, Matrix<Complex32>) {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var w = empty<Complex32>(N)
var v_read = empty<Complex32>(N, N)
var v_empty = empty<Complex32>(N, N)
let a_copy = a.copy()
var info: Int64
if (right) {
info = unsafe { LAPACKE_cgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), w.ptr, v_empty.ptr,
IntNative(N), v_read.ptr, IntNative(N)) }
} else {
info = unsafe { LAPACKE_cgeev(LAPACK_COL_MAJOR, LAPACK_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), w.ptr, v_read.ptr,
IntNative(N), v_empty.ptr, IntNative(N)) }
}
return (w, v_read)
}
public func eigvals(a: Matrix<Complex32>): Vector<Complex32> {
if (!a.isSquare()) {
throw IllegalArgumentException("eig: input matrix must be square")
}
let N: Int64 = a.getCols()
var w = empty<Complex32>(N)
var v_e1 = empty<Complex32>(N, N)
var v_e2 = empty<Complex32>(N, N)
let a_copy = a.copy()
var info = unsafe { LAPACKE_cgeev(LAPACK_COL_MAJOR, LAPACK_NO_VECTOR, LAPACK_NO_VECTOR,
IntNative(N), a_copy.ptr, IntNative(N), w.ptr, v_e1.ptr,
IntNative(N), v_e2.ptr, IntNative(N)) }
return w
}
@Test
public class TestEig {
@TestCase
func testLapackDgeev(): Unit {
testDgeev()
}
@TestCase
func testEigFloat64(): Unit {
var a = matrix<Float64>(
[[0.0, -1.0], [2.0, 0.0]]
)
var (evals, evecs) = eig(a)
var expected_evals = vector<Complex64>([
Complex64(0.0, 1.414214), Complex64(0.0, -1.414214)
])
var expected_evecs = matrix<Complex64>([
[Complex64(0.0, 0.577350), Complex64(0.0, -0.577350)],
[Complex64(0.816497, 0.0), Complex64(0.816497, 0.0)]
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-6))
@Assert(approxEqualC(evecs, expected_evecs, atol:1e-6))
var (evals2, evecs2) = eig(a, right:false)
var expected_evecs2 = matrix<Complex64>([
[Complex64(0.816497, 0.0), Complex64(0.816497, 0.0)],
[Complex64(0.0, -0.577350), Complex64(0.0, 0.577350)]
])
@Assert(approxEqualC(evals2, expected_evals, atol:1e-6))
@Assert(approxEqualC(evecs2, expected_evecs2, atol:1e-6))
}
@TestCase
func testEigvalsFloat64(): Unit {
var a = matrix<Float64>(
[[1.0, 4.0, 2.0, 4.0],
[0.25, 1.0, 0.5, 1.0],
[0.5, 2.0, 1.0, 0.5],
[0.25, 1.0, 2.0, 1.0]]
)
var evals = eigvals(a)
var expected_evals = vector<Complex64>([
Complex64(4.249226, 0.0), Complex64(-0.124613, 1.021513),
Complex64(-0.124613, -1.021513), Complex64(0.0, 0.0)
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-6))
}
@TestCase
func testEigFloat32(): Unit {
var a = matrix<Float32>(
[[0.0, -1.0], [2.0, 0.0]]
)
var (evals, evecs) = eig(a)
var expected_evals = vector<Complex32>([
Complex32(0.0, 1.414214), Complex32(0.0, -1.414214)
])
var expected_evecs = matrix<Complex32>([
[Complex32(0.0, 0.577350), Complex32(0.0, -0.577350)],
[Complex32(0.816497, 0.0), Complex32(0.816497, 0.0)]
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-5))
@Assert(approxEqualC(evecs, expected_evecs, atol:1e-5))
var (evals2, evecs2) = eig(a, right:false)
var expected_evecs2 = matrix<Complex32>([
[Complex32(0.816497, 0.0), Complex32(0.816497, 0.0)],
[Complex32(0.0, -0.577350), Complex32(0.0, 0.577350)]
])
@Assert(approxEqualC(evals2, expected_evals, atol:1e-5))
@Assert(approxEqualC(evecs2, expected_evecs2, atol:1e-5))
}
@TestCase
func testEigvalsFloat32(): Unit {
var a = matrix<Float32>(
[[1.0, 4.0, 2.0, 4.0],
[0.25, 1.0, 0.5, 1.0],
[0.5, 2.0, 1.0, 0.5],
[0.25, 1.0, 2.0, 1.0]]
)
var evals = eigvals(a)
var expected_evals = vector<Complex32>([
Complex32(4.249226, 0.0), Complex32(-0.124613, 1.021513),
Complex32(-0.124613, -1.021513), Complex32(0.0, 0.0)
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-5))
}
@TestCase
func testEigComplex64(): Unit {
var a = matrix<Complex64>(
[[Complex64(0.0, 0.0), Complex64(-1.0, 1.0)],
[Complex64(2.0, 0.0), Complex64(0.0, 0.0)]]
)
var (evals, evecs) = eig(a)
var expected_evals = vector<Complex64>([
Complex64(0.643594, 1.553774), Complex64(-0.643594, -1.553774)
])
var expected_evecs = matrix<Complex64>([
[Complex64(0.246293, 0.594604), Complex64(-0.246293, -0.594604)],
[Complex64(0.765367, 0.0), Complex64(0.765367, 0.0)]
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-6))
@Assert(approxEqualC(evecs, expected_evecs, atol:1e-6))
var (evals2, evecs2) = eig(a, right:false)
var expected_evecs2 = matrix<Complex64>([
[Complex64(0.765367, 0.0), Complex64(0.765367, 0.0)],
[Complex64(0.246293, -0.594604), Complex64(-0.246293, 0.594604)]
])
@Assert(approxEqualC(evals2, expected_evals, atol:1e-6))
@Assert(approxEqualC(evecs2, expected_evecs2, atol:1e-6))
}
@TestCase
func testEigvalsComplex64(): Unit {
var a = matrix<Complex64>(
[[Complex64(0.0, 0.0), Complex64(-1.0, 1.0)],
[Complex64(2.0, 0.0), Complex64(0.0, 0.0)]]
)
var evals = eigvals(a)
var expected_evals = vector<Complex64>([
Complex64(0.643594, 1.553774), Complex64(-0.643594, -1.553774)
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-6))
}
@TestCase
func testEigComplex32(): Unit {
var a = matrix<Complex32>(
[[Complex32(0.0, 0.0), Complex32(-1.0, 1.0)],
[Complex32(2.0, 0.0), Complex32(0.0, 0.0)]]
)
var (evals, evecs) = eig(a)
var expected_evals = vector<Complex32>([
Complex32(0.643594, 1.553774), Complex32(-0.643594, -1.553774)
])
var expected_evecs = matrix<Complex32>([
[Complex32(0.246293, 0.594604), Complex32(-0.246293, -0.594604)],
[Complex32(0.765367, 0.0), Complex32(0.765367, 0.0)]
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-5))
@Assert(approxEqualC(evecs, expected_evecs, atol:1e-5))
var (evals2, evecs2) = eig(a, right:false)
var expected_evecs2 = matrix<Complex32>([
[Complex32(0.765367, 0.0), Complex32(0.765367, 0.0)],
[Complex32(0.246293, -0.594604), Complex32(-0.246293, 0.594604)]
])
@Assert(approxEqualC(evals2, expected_evals, atol:1e-5))
@Assert(approxEqualC(evecs2, expected_evecs2, atol:1e-5))
}
@TestCase
func testEigvalsComplex32(): Unit {
var a = matrix<Complex32>(
[[Complex32(0.0, 0.0), Complex32(-1.0, 1.0)],
[Complex32(2.0, 0.0), Complex32(0.0, 0.0)]]
)
var evals = eigvals(a)
var expected_evals = vector<Complex32>([
Complex32(0.643594, 1.553774), Complex32(-0.643594, -1.553774)
])
@Assert(approxEqualC(evals, expected_evals, atol:1e-5))
}
@TestCase
func testLapackDstemr(): Unit {
testDstemr()
}
@TestCase
func testEighTridiagonal(): Unit {
var d = vector<Float64>(
[3.0, 3.0, 3.0, 3.0]
)
var e = vector<Float64>(
[-1.0, -1.0, -1.0, -1.0]
)
var (w, z) = eighTridiagonal(d, e)
var expected_w = vector<Float64>(
[1.381966, 2.381966, 3.618034, 4.618034]
)
var expected_z = matrix<Float64>(
[[0.371748, 0.601501, 0.601501, 0.371748],
[0.601501, 0.371748, -0.371748, -0.601501],
[0.601501, -0.371748, -0.371748, 0.601501],
[0.371748, -0.601501, 0.601501, -0.371748]]
)
@Assert(approxEqual(w, expected_w, atol:1e-6))
@Assert(approxEqual(z, expected_z, atol:1e-6))
}
@TestCase
func testEigvalshTridiagonal(): Unit {
var d = vector<Float64>(
[3.0, 3.0, 3.0, 3.0]
)
var e = vector<Float64>(
[-1.0, -1.0, -1.0, -1.0]
)
var (w, z) = eighTridiagonal(d, e)
var expected_w = vector<Float64>(
[1.381966, 2.381966, 3.618034, 4.618034]
)
@Assert(approxEqual(w, expected_w, atol:1e-6))
}
@TestCase
func testLapackDsbevx(): Unit {
testDsbevx()
}
// @TestCase
// func testLapackDsyevr(): Unit {
// testDsyevr()
// }
@TestCase
func testEigvalsh(): Unit {
var a = matrix<Float64>(
[[6.0, 3.0, 1.0, 5.0],
[3.0, 0.0, 5.0, 1.0],
[1.0, 5.0, 6.0, 2.0],
[5.0, 1.0, 2.0, 2.0]]
)
var w = eigvalsh(a)
var expected_w = vector<Float64>(
[-3.74637491, -0.76263923, 6.08502336, 12.42399079]
)
@Assert(approxEqual(w, expected_w, atol:1e-6))
}
// @TestCase
// func testEigh(): Unit {
// var a = matrix<Float64>(
// [[34.0, -4.0, -10.0, -7.0, 2.0],
// [-4.0, 7.0, 2.0, 12.0, 0.0],
// [-10.0, 2.0, 44.0, 2.0, -19.0],
// [-7.0, 12.0, 2.0, 79.0, -34.0],
// [2.0, 0.0, -19.0, -34.0, 29.0]]
// )
// var (w, z) = eigh(a)
// var expected_w = vector<Float64>(
// [0.000001, 9.119382, 31.587600, 51.433367, 100.859651]
// )
// var expected_z = matrix<Float64>(
// [[-0.042499, -0.243066, 0.892827, -0.349513, 0.140762],
// [ 0.662479, -0.721978, -0.158056, -0.029858, -0.118300],
// [-0.281040, -0.315582, 0.263317, 0.840945, -0.211894],
// [-0.353773, -0.155937, -0.070130, -0.367218, -0.843067],
// [-0.595967, -0.543836, -0.321907, -0.186846, 0.458836]]
// )
// @Assert(approxEqual(w, expected_w, atol:1e-6))
// @Assert(approxEqual(z, expected_z, atol:1e-6))
// }
@TestCase
func testLapackDsbevd(): Unit {
testDsbevd()
}
// @TestCase
// func testEigBanded(): Unit {
// var ab = matrix<Float64>(
// [[1.0, 2.0, 3.0, 4.0],
// [5.0, 5.0, 5.0, 0.0],
// [2.0, 2.0, 0.0, 0.0],
// [0.0, 0.0, 0.0, 0.0]]
// )
// var (w, z) = eigBanded(ab)
// var expected_w = vector<Float64>(
// [-4.262005, -2.229872, 3.952223, 12.539654]
// )
// var expected_z = matrix<Float64>(
// [[ 0.545851, -0.490263, -0.589435, 0.338015],
// [-0.734039, 0.048276, -0.411239, 0.538274],
// [ 0.398961, 0.671053, 0.158026, 0.604604],
// [-0.063753, -0.554075, 0.677109, 0.480063]]
// )
// @Assert(approxEqual(w, expected_w, atol:1e-6))
// @Assert(approxEqual(z, expected_z, atol:1e-6))
// }
@TestCase
func testEigvalsBanded(): Unit {
var ab = matrix<Float64>(
[[1.0, 2.0, 3.0, 4.0],
[5.0, 5.0, 5.0, 0.0],
[2.0, 2.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0]]
)
var w = eigvalsBanded(ab)
var expected_w = vector<Float64>(
[-4.262005, -2.229872, 3.952223, 12.539654]
)
@Assert(approxEqual(w, expected_w, atol:1e-6))
}
}